Ecosystem Services and Ecological Restoration in the Northern Shaanxi Loess Plateau, China, in Relation to Climate Fluctuation and Investments in Natural Capital

Accurately identifying the spatiotemporal variations and driving factors of ecosystem services (ES) in ecological restoration is important for ecosystem management and the sustainability of nature conservation strategies. As the Green for Grain project proceeds, food provision, water regulation and climate regulation services in the Northern Shaanxi Loess Plateau (NSLP) are changing and have caused broad attention. In this study, the dynamic pattern of the normalized differential vegetation index (NDVI) and the main drivers of grain production (GP), water yield (WY) and net primary production (NPP) in the NSLP from 2000–2013 are identified by incorporating multiple data and methods, in order to provide a better understanding of how and why ES change during ecological restoration. WY was simulated by hydrological modeling, and NPP was estimated with the Carnegie Ames Stanford Approach (CASA) model. The results show that vegetation restoration continued from 2000–2013, but fluctuated because of the comprehensive influence of climate and human activity. GP and NPP both exhibited significantly increasing trends, while changes in WY occurred in two stages: decline (2000–2006) and growth (2007–2013). Spatially, significantly increasing trends in NPP and WY were detected in 52.73% and 24.76% of the region, respectively, in areas that correspond with the Green for Grain project and high precipitation growth. Correlation and partial correlation analyses show that there were different dominant factors (i.e., natural vs. anthropogenic) driving ES change in the NSLP from 2000–2013. The change in WY was mainly driven by precipitation, while the improvements in GP and NPP can be attributed to investments in natural capital (i.e., chemical fertilizer, agricultural machinery power and afforestation). We also found that vegetation restoration can produce positive effects on NPP, but negative effects on WY by using response analyses of WY or NPP change to NDVI change, demonstrating that additional research on the role of water in vegetation restoration is needed. results provide support for ES management and the sustainable development of ecological restoration in the NSLP.


Introduction
China has been the world's second largest economy by the pursuit of rapid economic development, but the costs of this success are reflected in high levels of ecological degradation [1].In China, natural disasters (e.g., flood) caused by environmental degradation occur more frequently and have severely impacted sustainable socioeconomic development in recent decades [2].The Chinese government has implemented a series of ecological projects for improving the eco-environment.For example, Green for Grain (GFG), launched in 1999, is the largest eco-environmental construction project and important nature conservation program in China.GFG involves the effort to return cultivated land on steep slopes to perennial vegetation and afforest barren hills and wasteland using a government payment scheme, which engages millions of rural households as core agents of project implementation [2,3].The Northern Shaanxi Loess Plateau (NSLP) is the core region of GFG, and it is also an ecologically-vulnerable area sensitive to climatic change [4,5].Vegetation coverage of the NSLP has gradually improved since the implementation of GFG [6][7][8].This improvement in vegetation cover has realized huge ecological effects, including the potential in sequestrating greenhouse gases [9,10] or regulating climate, which is a hot issue in the background of climate change.However, because some regions of GFG (e.g., NSLP) belong to the semi-arid area, some research has raised concerns about the problems of food security and water scarcity caused by GFG [11,12].Climate regulation, water supply (yield) and food provision are all ecosystem services (ES), which are defined as the benefits people derive from ecosystems [13].As a new round of GFG projects begin, we need to evaluate the benefits and impacts of previous GFG projects, which have been of great significance to the sustainability of GFG.ES certainly will change under ecological restoration, and this fact has received significant attention in recent years [1,[14][15][16].ES can build a "bridge" between ecological systems and human social systems and can be used to evaluate the benefits and impacts of nature conservation strategies [17].In recent years, the Ecosystem Services Partnership, the Intergovernmental Science Policy Platform on Biodiversity and Ecosystem Services (IBPES), ecoSERVICES (a core project of Future Earth investigating the impact of biodiversity change on ecosystem functioning and services, and human well-being) and the Group on Earth Observations Biodiversity Observation Network (GEO-BON) have been established.One objective of these scientific programs is to assess ES comprehensively and promote the management and application of ES [18], and this requires knowledge of their dynamic features and driving factors.
As co-productions of natural and human systems [19,20], ES are affected by both natural factors (e.g., climatic change) and anthropogenic factors.According to the Millennium Ecosystem Assessment, the impact of climatic change on ES is growing [13].Because it influences ecosystem structure and function, climatic change can significantly alter the capacity of ecosystems to provide services [21,22], but the relationships between climatic change and various ES can be either positive or negative, and they require careful analysis [21,23,24].Humans can change ES in direct ways (e.g., land cover change, overexploitation, pollution) and indirect ways (e.g., population growth, economic activity, technical progress) [25,26].One of the most direct anthropogenic factors, land cover change, has been an important contributor to changes in ES over the past 50 years [13].Extensive croplands have been converted into woodlands, and natural vegetation has been restored in the NSLP, which has had a dramatic impact on ES.Current research shows that ES significantly changed in the Loess Plateau after GFG [3,27,28].However, the processes and main drivers of different ES changes during the whole GFG period in the NSLP still need to be further elucidated.Quantitative relationships between vegetation restoration and different ES are not yet fully understood.Here, we focus on food provision, water regulation and climate regulation in the NSLP, which are all highly related to the GFG project and interact with each other.By incorporating multi-source data and statistical methods, we investigated the dynamic patterns of vegetation and the main drivers of different ES changes in the NSLP from 2000-2013.In particular, we: (i) assessed the processes and the trends in vegetation restoration using the normalized difference vegetation index (NDVI); (ii) assessed the processes and trends in the difference of ES changes and their responses to vegetation restoration; and (iii) quantified the relationships between natural or anthropogenic factors and ES and identified the main drivers of different ES changes.

Study Area
Located in northern Shaanxi Province, China (35 • 21 -39 • 34 N and 107 • 41 -110 • 31 E) (Figure 1b), the NSLP covers an area of approximately 8.03 × 10 4 km 2 and consists of 25 counties (Figure 1a) belonging to two prefectures, Yan'an and Yulin (Figure 1c).High in the west and becoming lower in the east, the NSLP is a hilly, gullied area of the Loess Plateau [29,30].The Mu Us desert in the northwest of Yulin has a primarily psammophyte vegetation type.Huanglong Mountain and Ziwu Mountain in southern Yan'an are the main forest-covered areas in the NSLP, and vegetation coverage there is relatively high.The eastern section of the NSLP is flat, and land use categories there mainly consist of cropland and grassland.Crops are mainly corn, soybeans and wheat, and the grassland is dominated by Miscanthus grass and clover [29].This area has a semi-arid continental monsoon climate [30] and is characterized by warm and rainy summers and cold and dry winters [31].
Sustainability 2017, 9,199 3 of 20 Here, we focus on food provision, water regulation and climate regulation in the NSLP, which are all highly related to the GFG project and interact with each other.By incorporating multi-source data and statistical methods, we investigated the dynamic patterns of vegetation and the main drivers of different ES changes in the NSLP from 2000-2013.In particular, we: (i) assessed the processes and the trends in vegetation restoration using the normalized difference vegetation index (NDVI); (ii) assessed the processes and trends in the difference of ES changes and their responses to vegetation restoration; and (iii) quantified the relationships between natural or anthropogenic factors and ES and identified the main drivers of different ES changes.

Study Area
Located in northern Shaanxi Province, China (35°21'-39°34' N and 107°41'-110°31' E) (Figure 1b), the NSLP covers an area of approximately 8.03 × 10 4 km 2 and consists of 25 counties (Figure 1a) belonging to two prefectures, Yan'an and Yulin (Figure 1c).High in the west and becoming lower in the east, the NSLP is a hilly, gullied area of the Loess Plateau [29,30].The Mu Us desert in the northwest of Yulin has a primarily psammophyte vegetation type.Huanglong Mountain and Ziwu Mountain in southern Yan'an are the main forest-covered areas in the NSLP, and vegetation coverage there is relatively high.The eastern section of the NSLP is flat, and land use categories there mainly consist of cropland and grassland.Crops are mainly corn, soybeans and wheat, and the grassland is dominated by Miscanthus grass and clover [29].This area has a semi-arid continental monsoon climate [30] and is characterized by warm and rainy summers and cold and dry winters [31].From the maps of land use in 2000 (Figure S1a) and 2010 (Figure S1b), croplands, forestlands and grasslands showed the main land use types in our site.From a 2000-2010 land use transition matrix for the NSLP (Table 1), it can be seen that 694 km 2 and 556 km 2 of croplands in 2000 have been replaced by woodlands (including scrublands) and grasslands, respectively, due to the GFG project.Woodland coverage increased by 10.13%, while cropland coverage decreased by 4.19% during the 11-year period.At the same time, the grasslands increased by only 0.13%, because many grasslands were converted into woodlands.Dramatic changes in land use in the NSLP caused a series of ecological effects [2,32,33].From the maps of land use in 2000 (Figure S1a) and 2010 (Figure S1b), croplands, forestlands and grasslands showed the main land use types in our site.From a 2000-2010 land use transition matrix for the NSLP (Table 1), it can be seen that 694 km 2 and 556 km 2 of croplands in 2000 have been replaced by woodlands (including scrublands) and grasslands, respectively, due to the GFG project.Woodland coverage increased by 10.13%, while cropland coverage decreased by 4.19% during the 11-year period.At the same time, the grasslands increased by only 0.13%, because many grasslands were converted into woodlands.Dramatic changes in land use in the NSLP caused a series of ecological effects [2,32,33].

Vegetation Restoration Quantification
Vegetation indices are developed for evaluating vegetative covers using spectral measurements in the applications of remote sensing.Over forty vegetation indices have been developed in order to enhance vegetation response and minimize the effects of soil brightness, soil moisture, shadow and environmental effects [35].NDVI is the normalized difference vegetation index and calculated by reflectance in the near-infrared and red wavelengths [36,37].The NDVI is sensitive to green leaf area and surface vegetation coverage [4,36] and has been widely used in ecological research [37][38][39].A few research works [4,[6][7][8][36][37][38][39] have demonstrated the potential of using NDVI to detect vegetation dynamics.After the GFG project was implemented for more than ten years, the low-yield slope croplands were transformed into grasslands and woodlands, and many barren hills and wastelands turned "green".High-resolution images of remote sensing can detect these changes.MODIS (Moderate Resolution Imaging Spectroradiometer) NDVI with a 250-m spatial resolution can be adapted for the higher resolution detection.To eliminate the influence of cloud cover and atmosphere, the annual maximum value of the NDVI was chosen for use in this study by the maximum value composites methods [40].Through that, the continuity of NDVI time series was enhanced, and the temporal-spatial change of vegetation could be identified more accurately.

Ecosystem Services Quantification
According to the Millennium Ecosystem Assessment [13], more than twenty kinds of ES were listed and could be classified into provisioning services, regulating services, supporting services and cultural services.We only chose three typical and representative services, including food provision (provisioning service), hydrology regulation (regulating service) and climate regulation (regulating service) according to the availability of data, the feasibility of the method and the importance of ES.The explanation was as follows.Food provision service was directly related to human well-being.Moreover, food provision would be affected because of the reduction of cropland in our site.Hydrology regulation was a significant and key service because our site was located in a semi-arid area.Moreover, the increase of forest area might consume more water because of more intense transpiration.In our site, ecosystems were very sensitive to climate changes [4], so climate regulation also was an important service.Moreover, in recent years, carbon emission growth caused by human activities was 13.96% per year in this area [41], and the government undertook the huge responsibility of carbon emission reduction in our site.

Grain Production
Grain production (GP) was selected as a proxy for the capacity of the food provision services.Data on wheat, corn and rice production were obtained at the provincial level from the Bureaus of Statistics [42] in our site (Yulin and Yan'an).

Water Yield
Water yield (WY) was selected as a proxy for hydrology regulation services.In arid and semi-arid areas, water budgets mainly depend on precipitation (PPT) and evapotranspiration (ET) [43,44].WY is the result of the balancing of PPT and ET [29] and can be calculated as: PPT was derived from the National Climatic Bureau and interpolated with the kriging method [45].ET datasets were obtained from MOD16 ET products and estimated using the improved ET algorithm [46].
The interpolated precipitation data and the ET data (MOD16A3) came from different sources.We reduced the uncertainty through the following ways.The interpolated PPT data were set as the same spatial and temporal resolutions as satellite-derived ET data.The PPT and ET layers had the same coordinated system, and the calculation was processed in ArcGIS software with a stable environment.A few extreme values in the results were removed, and they were replaced by adjacent values.

Net Primary Production
Net primary production (NPP) is the key driver of the terrestrial carbon cycling, balancing CO 2 and O 2 and regulating global temperature [33].Feng et al. [3] have found that the increase of carbon sequestration was mainly contributed by enhanced NPP on the Loess Plateau during the period of GFG implementation.The increment of NPP would enhance climate regulation services [47] in the NSLP, which is sensitive to climate changes.In this research, NPP was selected as a proxy for climate regulation services.The Carnegie Ames Stanford Approach (CASA) was employed to estimate per-pixel NPP at monthly intervals.CASA is a typical model of light use efficiency [48], and it needs relatively few parameters.Moreover, some parameters do not require actual measurement and can be obtained from remote sensing data.The total formula estimating NPP with CASA can be represented as [49]: where SOL is the total solar radiation (MJ/m 2 ), which can be calculated from Equations ( 4) and ( 5 Solar radiation can be calculated with the Angstrom formula that relates solar radiation to extraterrestrial radiation and sunshine percentage [53].The local intensity of extraterrestrial radiation is determined by the angle between the direction of the Sun's rays and the normal to the surface of the atmosphere [53].The key equations for the solar radiation calculation are as follows: where ESOL is the extraterrestrial radiation (MJ/m 2 ); S c is the solar constant (1360 W/m 2 or 49.20 MJ/m 2 •h 2 ); d r , w s and β are the inverse relative Earth-Sun distance, sunset hour angle and solar decimation, respectively, calculated following Allen et al. [53]; α is the latitude of the meteorological stations (rad); a s and b s are regression constants and variable in different latitudes, which were set according to the simulation results of Wen [54] in north China as follows: 0.344 for a s and 0.390 for b s ; n is the actual duration of sunshine; N is daylight hours; and n/N is the relative sunshine duration.

Change Trends
Linear regression analysis was used to simulate NDVI/WY/NPP change trends at the pixel level or over the whole regional level [4,55,56].The equation is as follows: where a positive or negative Slope indicates an increasing or decreasing trend for NDVI/WY/NPP; n is the number of years; and PA t is the corresponding NDVI/WY/NPP value.The F-test method was used to analyze the statistical significance of the NDVI/WY/NPP change.

Response Relationships
The linear regression analysis was also used to explore the response relationships of ES change to vegetation restoration at the pixel scale.This paper is mainly focused on the response of regulating services to vegetation restoration, which has not been fully understood in current research [25][26][27][28]57].The responses of WY change and NPP change to NDVI change at the pixel level were calculated using linear regression analysis.An increased or decreased result indicates a positive or negative response.

Correlation and Partial Correlation Analysis
Correlation and partial correlation analysis were used to analyze the relationships between ES and different driving factors.Both natural and anthropogenic factors have an important impact on ES [13].In this paper, the climate fluctuations were chosen to be the main natural driving factors.The anthropogenic driving factors in our site were different for provisioning service (GP) and regulating services (WY and NPP).Temperature (T), precipitation (PPT), agricultural machinery power consumption, chemical fertilizer consumption and labor force input were chosen to be the driving factors of GP change, which was based on the study of Zhang et al. [58].T, PPT and cumulative afforested area were chosen to be the driving factors of WY and NPP change, which agreed with the research of Xie et al. [59].The cumulative afforested area could be understood as an indicator of GFG policy.
Partial correlation analysis is used when the two variables are related to the other variables at the same time; it will exclude the impact of the other variables and only analyze the correlation of the two variables [4,60].For example, if there are three variables, the equation is as follows: where R xy-z is the partial correlation coefficient between x and y excluding the impact of z; R xy is the correlation coefficient between x and y; R xz is the correlation coefficient between x and z; R yz is the correlation coefficient between y and z.The calculation of the correlation coefficient R xy (or R xz or R xz ) is detailed in Xu [60].

Data Preparation
Meteorological data from 32 weather stations (Table S1) within and near the NSLP were collected from the National Climatic Bureau [61], including monthly T, monthly PPT and sunshine duration during the study period.
The tile numbers of the remote sensing data in the NSLP were h26v05 and h27v05.A series of Moderate Resolution Imaging Spectroradiometer (MODIS) data for the CASA model during the 14-year period was derived from the MODIS on NASA's Terra satellite [62], including surface albedo (MOD09A1), vegetation index (MOD13A2) and land cover type (MCD12Q1).The annual ET product (MOD16A3) was used to estimate WY, and the annual NPP product (MOD17A3) was used for comparison with the simulated NPP results of the CASA model.MOD16A3 and MOD17A3 were derived from the Numerical Terradynamic Simulation Group [63].
Data on agricultural machinery power consumption, chemical fertilizer consumption, labor force input in GP and cumulative afforested area during the 14-year period were obtained from provincial-level Bureaus of Statistics [42] in the NSLP and used for quantifying the relationships between ES and different driving factors.

Climate Change Background
According to the original meteorological station data within our site, the area exhibited a warming and drying trend during the period from 1955-2013.We divided the complete interval into two periods, 1955-1999 and 2000-2013.From 1955-1999, we found that the annual mean T was 8.97 • C with an increase of 0.019 • C/year, and the annual mean PPT was 469.84 mm with a decrease of 2.61 mm/year (Figure 2).From 2000-2013, the annual mean T was 9.78 • C with a decrease of 0.024 • C/year, and the annual mean PPT was 433.48 mm with an increase of 8.46 mm/year (Figure 2).The annual mean T was higher and PPT lower from 2000-2013, while the changes of T and PPT were opposite to those of 1955-1999.At the pixel level, T and PPT both had higher values in the south than the north (Figure S2).In the last 14-year period, the T of 59.32% of the area was decreasing, while the PPT of the whole area was increasing (Figure S2).The global surface T was increasing over the past few decades, while the situation seemed to be different in our site.GFG policy might play a part of the role, because the most T decreasing regions were distributed in the GFG regions (Figure S2).

Data Preparation
Meteorological data from 32 weather stations (Table S1) within and near the NSLP were collected from the National Climatic Bureau [61], including monthly T, monthly PPT and sunshine duration during the study period.
The tile numbers of the remote sensing data in the NSLP were h26v05 and h27v05.A series of Moderate Resolution Imaging Spectroradiometer (MODIS) data for the CASA model during the 14-year period was derived from the MODIS on NASA's Terra satellite [62], including surface albedo (MOD09A1), vegetation index (MOD13A2) and land cover type (MCD12Q1).The annual ET product (MOD16A3) was used to estimate WY, and the annual NPP product (MOD17A3) was used for comparison with the simulated NPP results of the CASA model.MOD16A3 and MOD17A3 were derived from the Numerical Terradynamic Simulation Group [63].
Data on agricultural machinery power consumption, chemical fertilizer consumption, labor force input in GP and cumulative afforested area during the 14-year period were obtained from provincial-level Bureaus of Statistics [42] in the NSLP and used for quantifying the relationships between ES and different driving factors.

Climate Change Background
According to the original meteorological station data within our site, the area exhibited a warming and drying trend during the period from 1955-2013.We divided the complete interval into two periods, 1955-1999 and 2000-2013.From 1955-1999, we found that the annual mean T was 8.97 °C with an increase of 0.019 °C/year, and the annual mean PPT was 469.84 mm with a decrease of 2.61 mm/year (Figure 2).From 2000-2013, the annual mean T was 9.78 °C with a decrease of 0.024 °C/year, and the annual mean PPT was 433.48 mm with an increase of 8.46 mm/year (Figure 2).The annual mean T was higher and PPT lower from 2000-2013, while the changes of T and PPT were opposite to those of 1955-1999.At the pixel level, T and PPT both had higher values in the south than the north (Figure S2).In the last 14-year period, the T of 59.32% of the area was decreasing, while the PPT of the whole area was increasing (Figure S2).The global surface T was increasing over the past few decades, while the situation seemed to be different in our site.GFG policy might play a part of the role, because the most T decreasing regions were distributed in the GFG regions (Figure S2).

Vegetation Change
From 2000-2013, the NDVI exhibited a significantly increasing trend of 0.01/year (R 2 = 0.858, p < 0.05), and thus, it improved 33.02% from 2000-2013 (Figure 3).From 2000-2002, the NDVI was in a period of rapid growth with an increase of 0.02/year, because the scale of GFG was larger in the

Vegetation Change
From 2000-2013, the NDVI exhibited a significantly increasing trend of 0.01/year (R 2 = 0.858, p < 0.05), and thus, it improved 33.02% from 2000-2013 (Figure 3).From 2000-2002, the NDVI was in a period of rapid growth with an increase of 0.02/year, because the scale of GFG was larger in the early days of the project [4].From 2003-2007, the NDVI was in a fluctuating period and kept decreasing from 2005-2007.Some vegetation (e.g., Hippophae rhamnoides L.) was planted, but it was not suitable for the semi-arid climate in the NSLP and needed more water when growing [64].When PPT was low, these plants started to wither.More suitable vegetation was subsequently planted (e.g., Pinus sylvestris), and the NDVI started increasing from 2008.
PPT was low, these plants started to wither.More suitable vegetation was subsequently planted (e.g., Pinus sylvestris), and the NDVI started increasing from 2008.
In spatial terms, the NDVI in 85.62% of the area showed an increase (Figure 4b), and the NDVI in 52.96% of the area showed a significantly increasing trend (p < 0.05) (Figure 4c).Those areas were in the mid-eastern NSLP, which originally had low vegetation coverage (Figure 4a) and was therefore heavily involved in the GFG program from 2000-2013.At the same time, the NDVI in 14.38% of the area showed a slight decrease (Figure 4b), and these decreases were mostly not significant.The decreased areas were located in the southern NSLP (i.e., Ziwu Mountain and Huanglong Mountain) with high forest coverage.Some other researchers have found similar results either on larger scales [4,6] or at the same scale [64,65].

GP Change
During 2000-2013, GP increased by 1.63-times and showed a significantly increasing trend of 24.87 t/km 2 •year (R 2 = 0.824, p < 0.01) (Figure 5a).The croplands area decreased because of the GFG program, so the gross productivity increased only 1.3-times and showed a growth of 16.84 × 10 4 t/year (Figure 5a).In 2005 and 2006, GP decreased (Figure 5a) because of drought stress (Figure 2b).A significant improvement in GP and gross productivity also was detected by other researchers [2,66,67].In spatial terms, the NDVI in 85.62% of the area showed an increase (Figure 4b), and the NDVI in 52.96% of the area showed a significantly increasing trend (p < 0.05) (Figure 4c).Those areas were in the mid-eastern NSLP, which originally had low vegetation coverage (Figure 4a) and was therefore heavily involved in the GFG program from 2000-2013.At the same time, the NDVI in 14.38% of the area showed a slight decrease (Figure 4b), and these decreases were mostly not significant.The decreased areas were located in the southern NSLP (i.e., Ziwu Mountain and Huanglong Mountain) with high forest coverage.Some other researchers have found similar results either on larger scales [4,6] or at the same scale [64,65].early days of the project [4].From 2003-2007, the NDVI was in a fluctuating period and kept decreasing from 2005-2007.Some vegetation (e.g., Hippophae rhamnoides L.) was planted, but it was not suitable for the semi-arid climate in the NSLP and needed more water when growing [64].When PPT was low, these plants started to wither.More suitable vegetation was subsequently planted (e.g., Pinus sylvestris), and the NDVI started increasing from 2008.
In spatial terms, the NDVI in 85.62% of the area showed an increase (Figure 4b), and the NDVI in 52.96% of the area showed a significantly increasing trend (p < 0.05) (Figure 4c).Those areas were in the mid-eastern NSLP, which originally had low vegetation coverage (Figure 4a) and was therefore heavily involved in the GFG program from 2000-2013.At the same time, the NDVI in 14.38% of the area showed a slight decrease (Figure 4b), and these decreases were mostly not significant.The decreased areas were located in the southern NSLP (i.e., Ziwu Mountain and Huanglong Mountain) with high forest coverage.Some other researchers have found similar results either on larger scales [4,6] or at the same scale [64,65].

GP Change
During 2000-2013, GP increased by 1.63-times and showed a significantly increasing trend of 24.87 t/km 2 •year (R 2 = 0.824, p < 0.01) (Figure 5a).The croplands area decreased because of the GFG program, so the gross productivity increased only 1.3-times and showed a growth of 16.84 × 10 4 t/year (Figure 5a).In 2005 and 2006, GP decreased (Figure 5a) because of drought stress (Figure 2b).A significant improvement in GP and gross productivity also was detected by other researchers [2,66,67].5a).The croplands area decreased because of the GFG program, so the gross productivity increased only 1.3-times and showed a growth of 16.84 × 10 4 t/year (Figure 5a).In 2005 and 2006, GP decreased (Figure 5a) because of drought stress (Figure 2b).A significant improvement in GP and gross productivity also was detected by other researchers [2,66,67].5b).Overall, WY in our site showed only a slight increase from 2000-2013, which was consistent with the PPT (Figure 2b).Other research [2] showed WY decreased after GFG in the whole Loess Plateau from 2000-2008.The results in this study show that WY decreased 24.81 mm during the same period.
Substantial regional variation in WY was detected under the influence of climate and vegetation coverage (Figure 6a).From the spatial change patterns, we see that 90% of the area experienced an increase in WY (Figure 6b).A decrease in WY was identified in the northwest NSLP, including Dingbian County, Jingbian County and Wuqi County.The annual growth rate of PPT was low in those areas.The change in 75.24% of the area was not significant, and 24.76% of the area experienced significant growth (p < 0.05), as seen in northern Yulin, including Shenmu County and Fugu County (Figure 6c), where there was a high annual growth rate of PPT (>8 mm/year) (Figure S2).5b).Overall, WY in our site showed only a slight increase from 2000-2013, which was consistent with the PPT (Figure 2b).Other research [2] showed WY decreased after GFG in the whole Loess Plateau from 2000-2008.The results in this study show that WY decreased 24.81 mm during the same period.
Substantial regional variation in WY was detected under the influence of climate and vegetation coverage (Figure 6a).From the spatial change patterns, we see that 90% of the area experienced an increase in WY (Figure 6b).A decrease in WY was identified in the northwest NSLP, including Dingbian County, Jingbian County and Wuqi County.The annual growth rate of PPT was low in those areas.The change in 75.24% of the area was not significant, and 24.76% of the area experienced significant growth (p < 0.05), as seen in northern Yulin, including Shenmu County and Fugu County (Figure 6c), where there was a high annual growth rate of PPT (>8 mm/year) (Figure S2).5b).Overall, WY in our site showed only a slight increase from 2000-2013, which was consistent with the PPT (Figure 2b).Other research [2] showed WY decreased after GFG in the whole Loess Plateau from 2000-2008.The results in this study show that WY decreased 24.81 mm during the same period.
Substantial regional variation in WY was detected under the influence of climate and vegetation coverage (Figure 6a).From the spatial change patterns, we see that 90% of the area experienced an increase in WY (Figure 6b).A decrease in WY was identified in the northwest NSLP, including Dingbian County, Jingbian County and Wuqi County.The annual growth rate of PPT was low in those areas.The change in 75.24% of the area was not significant, and 24.76% of the area experienced significant growth (p < 0.05), as seen in northern Yulin, including Shenmu County and Fugu County (Figure 6c), where there was a high annual growth rate of PPT (>8 mm/year) (Figure S2).5c).The decline in NPP may be attributed to the death of vegetation, as we have previously described [64].Overall, NPP in our site showed a significant increasing trend during the 14-year period (R 2 = 0.542, p < 0.01), which indicated that NPP in local ecosystems had improved.
NPP showed a similar spatial pattern to the NDVI (Figure 4a,b and Figure 6d,e).An increase was observed in 88.79% of the area, principally the mid-eastern NSLP that had undergone ecological restoration through the GFG program.NPP in Ziwu Mountain and Huanglong Mountain areas had a decrease, and this may be attributed to environmental factors and natural succession [68].The higher significance level indicates a more credible NPP variation trend.The increasing NPP trend was significant in 52.73% of the area, while the decreasing NPP was generally not significant (Figure 6f).

Response of ES Change to Vegetation Change
Response of regulating services to vegetation change is an important topic with the GFG project proceeding, but it has not been fully understood in current research [25][26][27][28]57].At the pixel scale in our site, response relationships between regulating services (WY and NPP) and NDVI were analyzed and identified for the period from 2000-2013 (Figure 7).The response of WY change to NDVI change was negative (Figure 7a) (R 2 = 0.012, R = −0.111,p < 0.01), which showed a negative relationship between WY and vegetation restoration.Vegetation restoration may have increased evapotranspiration (ET) and reduced regional WY, as has been found by other researchers [12,69,70].The response of NPP change to NDVI change was positive (Figure 7b), which showed a positive relationship between NPP and vegetation restoration and that was consistent with previous research [2,59,68].Moreover, Figure 7 showed that the response levels of regulating services to vegetation change were NPP > WY.  5c).The decline in NPP may be attributed to the death of vegetation, as we have previously described [64].Overall, NPP in our site showed a significant increasing trend during the 14-year period (R 2 = 0.542, p < 0.01), which indicated that NPP in local ecosystems had improved.
NPP showed a similar spatial pattern to the NDVI (Figures 4a,b and 6d,e).An increase was observed in 88.79% of the area, principally the mid-eastern NSLP that had undergone ecological restoration through the GFG program.NPP in Ziwu Mountain and Huanglong Mountain areas had a decrease, and this may be attributed to environmental factors and natural succession [68].The higher significance level indicates a more credible NPP variation trend.The increasing NPP trend was significant in 52.73% of the area, while the decreasing NPP was generally not significant (Figure 6f).

Response of ES Change to Vegetation Change
Response of regulating services to vegetation change is an important topic with the GFG project proceeding, but it has not been fully understood in current research [25][26][27][28]57].At the pixel scale in our site, response relationships between regulating services (WY and NPP) and NDVI were analyzed and identified for the period from 2000-2013 (Figure 7).The response of WY change to NDVI change was negative (Figure 7a) (R 2 = 0.012, R = −0.111,p < 0.01), which showed a negative relationship between WY and vegetation restoration.Vegetation restoration may have increased evapotranspiration (ET) and reduced regional WY, as has been found by other researchers [12,69,70].The response of NPP change to NDVI change was positive (Figure 7b), which showed a positive relationship between NPP and vegetation restoration and that was consistent with previous research [2,59,68].Moreover, Figure 7 showed that the response levels of regulating services to vegetation change were NPP > WY.

Relationships between ES and Different Driving Factors
From the results of correlation coefficients and partial correlation coefficients, both chemical fertilizer and agricultural machinery power were significant positive with GP from 2000-2013 in our site (Table 2).The GP correlation level of chemical fertilizer was higher than that of agricultural machinery power.The correlation between natural factors (T and PPT) and GP both were not significant (Table 2).A negative correlation was identified between labor force and GP from 2000-2013.During the study period, the labor force input in GP kept decreasing.However, other investments in GP including chemical fertilizer and agricultural machinery power were increasing, so the GP still kept increasing in our site.We divided the whole period into two stages: a

Relationships between ES and Different Driving Factors
From the results of correlation coefficients and partial correlation coefficients, both chemical fertilizer and agricultural machinery power were significant positive with GP from 2000-2013 in our site (Table 2).The GP correlation level of chemical fertilizer was higher than that of agricultural machinery power.The correlation between natural factors (T and PPT) and GP both were not significant (Table 2).A negative correlation was identified between labor force and GP from 2000-2013.During the study period, the labor force input in GP kept decreasing.However, other investments in GP including chemical fertilizer and agricultural machinery power were increasing, so the GP still kept increasing in our site.We divided the whole period into two stages: a GFG-extensive period (2000)(2001)(2002)(2003)(2004)(2005)(2006) and a GFG-weakened period (2007)(2008)(2009)(2010)(2011)(2012)(2013).It seemed that the correlation between different factors (natural or anthropogenic factors) and GP was stronger in the GFG-weakened period than in the GFG-extensive period (Table 2).Over the whole regional scale, neither T nor cumulative afforested area were significantly correlated with WY from 2000-2013, but a significant positive correlation was identified between WY and PPT (Table 3).The strong response of WY change to PPT change at the pixel scale also showed that WY was significantly correlated with PPT (Figure 8a).The correlation between natural factors and WY was stronger in the GFG-weakened period than in the GFG-extensive period (Table 3).Through the correlation coefficients between WY and cumulative afforested area in different periods, it seemed that we could not judge that WY was increased or decreased with the increasing of cumulative afforested area.However, the partial coefficients in different periods were consistent and indicated that WY was decreased with the increasing of cumulative afforested area.GFG-extensive period (2000)(2001)(2002)(2003)(2004)(2005)(2006) and a GFG-weakened period (2007)(2008)(2009)(2010)(2011)(2012)(2013).It seemed that the correlation between different factors (natural or anthropogenic factors) and GP was stronger in the GFG-weakened period than in the GFG-extensive period (Table 2).Over the whole regional scale, neither T nor cumulative afforested area were significantly correlated with WY from 2000-2013, but a significant positive correlation was identified between WY and PPT (Table 3).The strong response of WY change to PPT change at the pixel scale also showed that WY was significantly correlated with PPT (Figure 8a).The correlation between natural factors and WY was stronger in the GFG-weakened period than in the GFG-extensive period (Table 3).Through the correlation coefficients between WY and cumulative afforested area in different periods, it seemed that we could not judge that WY was increased or decreased with the increasing of cumulative afforested area.However, the partial coefficients in different periods were consistent and indicated that WY was decreased with the increasing of cumulative afforested area.Table 3. Correlation and partial correlation between WY (water yield) and T/PPT/CAA (cumulative afforested area) during the study period over the whole regional scale.Over the whole regional scale, neither T nor PPT were significantly correlated with NPP (Table 4).NPP was decreased with the increasing in T, and the response of NPP change to T change at the pixel scale could indicate that (Figure 8b).The NSLP is a semi-arid area, so higher T would intensify soil evaporation and restrain plant photosynthesis [4].The response of NPP change to T change at the pixel scale indicated that NPP was increased with the increasing in PPT (Figure 8c).Vegetation growth was restrained by moisture in the NSLP, so higher PPT enhanced plant activity and vegetation growth.During the 14-year period, the correlation between NPP and cumulative afforested area was significant and positive (Table 4).A stronger correlation was identified between NPP and cumulative afforested area during the GFG-extensive period than the GFG-weakened period (Table 4).

Credibility of the Results
NPP and ET are crucial to the process of assessing ES change.NPP was estimated by the CASA model using multi-source data.ET was derived from the MODIS product with strong variability and it was a key and uncertain parameter in simulating WY.NPP and ET are affected by vegetation physiological factors and complex environment factors, so it is necessary to evaluate the credibility of the NPP and ET values.
The several simulated NPP values in the Loess Plateau have differences, and the main factors responsible include the simulation method, the study period and the scope of the Loess Plateau (Table 5) [71].With respect to the scope of the Loess Plateau, annual NPP values simulated by Biome-BGC (Biogeochemical Cycles) were 10.08% lower for our site than for the whole area of the Loess Plateau.The simulated NPP value from the CASA model was lower by 11.69%-12.26%for our site than for the whole area of the Loess Plateau.This comparison showed that the CASA model used for our site performed well.Overall, the NPP values simulated by Biome-BGC were lower than the CASA values, partly because CASA simulates NPP at monthly intervals, and its estimates of the fraction of photosynthetic active radiation may therefore be larger [59].The average NPP estimated by the CASA model (272 g C/m 2 •year) was higher than the MOD17-NPP (214 g C/m 2 •year) (Table 5).In the 25 counties of our site, simulated NPP values ranged from 141 g C/m 2 •year-695 g C/m 2 •year, and MOD17-NPP values ranged from 148 g C/m 2 •year-357 g C/m 2 •year.The correlation between simulated NPP and MOD17-NPP is significant in all 25 counties (Figure S3a).However, the simulated NPP was much higher than the MOD17-NPP in Fu County, Yichuan County, and especially in Huanglong County and Huangling County (Figure 9 and Figure S3).The above four counties are primarily located near Ziwu Mountain and Huanglong Mountain, which have high vegetation coverage (Figure 4a), so the NPP value should be higher there than in other counties, and some research corroborates this [3,68].MOD17-NPP values were in middle range in the four counties, which shows that NPP simulated by the CASA model is more reasonable in those areas.Figure 9 shows that MOD17-NPP had small fluctuations in the 25 counties, which was mainly due to the simulation of MOD17-NPP at a global scale, while the NPP as estimated by CASA exhibited reasonable variation.During the 14-year period over the whole area, the annual simulated NPP and MOD17-NPP also had a significant correlation (R 2 = 0.61, p < 0.01) (Figure S3b).Overall, the NPP simulated by CASA in this study is reliable.
Sustainability 2017, 9, 199 13 of 20 148 g C/m 2 •year-357 g C/m 2 •year.The correlation between simulated NPP and MOD17-NPP is significant in all 25 counties (Figure S3a).However, the simulated NPP was much higher than the MOD17-NPP in Fu County, Yichuan County, and especially in Huanglong County and Huangling County (Figure 9 and Figure S3).The above four counties are primarily located near Ziwu Mountain and Huanglong Mountain, which have high vegetation coverage (Figure 4a), so the NPP value should be higher there than in other counties, and some research corroborates this [3,68].MOD17-NPP values were in middle range in the four counties, which shows that NPP simulated by the CASA model is more reasonable in those areas.Figure 9 shows that MOD17-NPP had small fluctuations in the 25 counties, which was mainly due to the simulation of MOD17-NPP at a global scale, while the NPP as estimated by CASA exhibited reasonable variation.During the 14-year period over the whole area, the annual simulated NPP and MOD17-NPP also had a significant correlation (R 2 = 0.61, p < 0.01) (Figure S3b).Overall, the NPP simulated by CASA in this study is reliable.In order to demonstrate the reliability of the ET dataset in our site, we offer the following points.First, the validity of MOD16-ET has been verified on different scales, as has been shown for Asia [73], Korea [74], China [75], Shaanxi Provence [76] and the Wei River Basin [77].Second, the synergy between NPP and ET, which has been found in our site by Jia et al. [33], means that ET can be enhanced with the increases in NPP.At the pixel scale, the correlation between annual NPP and ET from 2000-2013 was analyzed and found to be significantly positive (R 2 > 0.8, p < 0.01) (Figure 10a).The response of the ET change trend to the NPP change trend was also analyzed, and it too was significantly positive (R 2 = 0.282, p < 0.01) (Figure 10b).As a result of this positive correlation between NPP and ET, the reliability of MOD16-ET in the NSLP is further enhanced.In order to demonstrate the reliability of the ET dataset in our site, we offer the following points.First, the validity of MOD16-ET has been verified on different scales, as has been shown for Asia [73], Korea [74], China [75], Shaanxi Provence [76] and the Wei River Basin [77].Second, the synergy between NPP and ET, which has been found in our site by Jia et al. [33], means that ET can be enhanced with the increases in NPP.At the pixel scale, the correlation between annual NPP and ET from 2000-2013 was analyzed and found to be significantly positive (R 2 > 0.8, p < 0.01) (Figure 10a).The response of the ET change trend to the NPP change trend was also analyzed, and it too was significantly positive (R 2 = 0.282, p < 0.01) (Figure 10b).As a result of this positive correlation between NPP and ET, the reliability of MOD16-ET in the NSLP is further enhanced.
Sustainability 2017, 9,199 13 of 20 148 g C/m 2 •year-357 g C/m 2 •year.The correlation between simulated NPP and MOD17-NPP is significant in all 25 counties (Figure S3a).However, the simulated NPP was much higher than the MOD17-NPP in Fu County, Yichuan County, and especially in Huanglong County and Huangling County (Figure 9 and Figure S3).The above four counties are primarily located near Ziwu Mountain and Huanglong Mountain, which have high vegetation coverage (Figure 4a), so the NPP value should be higher there than in other counties, and some research corroborates this [3,68].MOD17-NPP values were in middle range in the four counties, which shows that NPP simulated by the CASA model is more reasonable in those areas.Figure 9 shows that MOD17-NPP had small fluctuations in the 25 counties, which was mainly due to the simulation of MOD17-NPP at a global scale, while the NPP as estimated by CASA exhibited reasonable variation.During the 14-year period over the whole area, the annual simulated NPP and MOD17-NPP also had a significant correlation (R 2 = 0.61, p < 0.01) (Figure S3b).Overall, the NPP simulated by CASA in this study is reliable.In order to demonstrate the reliability of the ET dataset in our site, we offer the following points.First, the validity of MOD16-ET has been verified on different scales, as has been shown for Asia [73], Korea [74], China [75], Shaanxi Provence [76] and the Wei River Basin [77].Second, the synergy between NPP and ET, which has been found in our site by Jia et al. [33], means that ET can be enhanced with the increases in NPP.At the pixel scale, the correlation between annual NPP and ET from 2000-2013 was analyzed and found to be significantly positive (R 2 > 0.8, p < 0.01) (Figure 10a).The response of the ET change trend to the NPP change trend was also analyzed, and it too was significantly positive (R 2 = 0.282, p < 0.01) (Figure 10b).As a result of this positive correlation between NPP and ET, the reliability of MOD16-ET in the NSLP is further enhanced.

Main Drivers of ES Changes
According to the relationships between GP and different factors (Table 2), neither T nor PPT were significantly correlated with GP during the study period, but the correlations between GP and the investments in agriculture (i.e., chemical fertilizer and agricultural machinery power) were significant (Table 2).Fu et al. [78] considered GP to be closely related to agricultural production conditions, artificial inputs and technical progress in the Loess Plateau.Our study showed that the increasing GP was mainly driven by the increasing agricultural investments in our site.Chemical fertilizer investment ranged from 155.71 kg/ha in 2000 to 330.29 kg/ha in 2013 (Figure 11a); it thus grew 1.12-times from 2000-2013 with a significantly increasing trend of 15.39 kg/ha•year (R 2 = 0.928, p < 0.01).Agricultural machinery power grew 2.07-times and ranged from 19.80 kW/ha in 2000 to 60.77 kW/ha in 2013, with a significantly increasing trend (R 2 = 0.996, p < 0.01) (Figure 11b).Moreover, our study also showed that chemical fertilizer input was more effective than agricultural machinery power input (Table 2).Other minor drivers of GP improvement in our site may include farm management and the construction of high quality basic croplands [2,65,66].

Main Drivers of ES Changes
According to the relationships between GP and different factors (Table 2), neither T nor PPT were significantly correlated with GP during the study period, but the correlations between GP and the investments in agriculture (i.e., chemical fertilizer and agricultural machinery power) were significant (Table 2).Fu et al. [78] considered GP to be closely related to agricultural production conditions, artificial inputs and technical progress in the Loess Plateau.Our study showed that the increasing GP was mainly driven by the increasing agricultural investments in our site.Chemical fertilizer investment ranged from 155.71 kg/ha in 2000 to 330.29 kg/ha in 2013 (Figure 11a); it thus grew 1.12-times from 2000-2013 with a significantly increasing trend of 15.39 kg/ha•year (R 2 = 0.928, p < 0.01).Agricultural machinery power grew 2.07-times and ranged from 19.80 kW/ha in 2000 to 60.77 kW/ha in 2013, with a significantly increasing trend (R 2 = 0.996, p < 0.01) (Figure 11b).Moreover, our study also showed that chemical fertilizer input was more effective than agricultural machinery power input (Table 2).Other minor drivers of GP improvement in our site may include farm management and the construction of high quality basic croplands [2,65,66].According to the relationships between WY and different factors (Table 3), both T and cumulative afforested area had a weak connection with WY change from 2000-2013 (Table 3).However, the negative effect of afforestation was significant in the GFG-extensive period (Table 3).The response of WY change to NDVI change showed there was a negative relationship between vegetation restoration and WY (Figure 7a).In this study, WY was calculated by subtracting ET from PPT.During the study period, ET increased with the growth of woodlands (Figure 11c), which can partly explain the negative relationship between afforestation or vegetation restoration and WY.However, the increase in ET was slight, and the annual growth in PPT (8.46 mm/year) was much larger than the growth in ET (1.90 mm/year), which suggests that the slight improvement of WY was mainly driven by PPT during the study period.The strong and positive response of WY change to PPT change (Figure 8a), significant correlation (Table 3) and linear fitting (Figure 12a) between WY and PPT can enhance the credibility of that conclusion.However, the afforestation producing the negative effect on WY still cannot be ignored.According to the relationships between WY and different factors (Table 3), both T and cumulative afforested area had a weak connection with WY change from 2000-2013 (Table 3).However, the negative effect of afforestation was significant in the GFG-extensive period (Table 3).The response of WY change to NDVI change showed there was a negative relationship between vegetation restoration and WY (Figure 7a).In this study, WY was calculated by subtracting ET from PPT.During the study period, ET increased with the growth of woodlands (Figure 11c), which can partly explain the negative relationship between afforestation or vegetation restoration and WY.However, the increase in ET was slight, and the annual growth in PPT (8.46 mm/year) was much larger than the growth in ET (1.90 mm/year), which suggests that the slight improvement of WY was mainly driven by PPT during the study period.The strong and positive response of WY change to PPT change (Figure 8a), significant correlation (Table 3) and linear fitting (Figure 12a) between WY and PPT can enhance the credibility of that conclusion.However, the afforestation producing the negative effect on WY still cannot be ignored.

Main Drivers of ES Changes
According to the relationships between GP and different factors (Table 2), neither T nor PPT were significantly correlated with GP during the study period, but the correlations between GP and the investments in agriculture (i.e., chemical fertilizer and agricultural machinery power) were significant (Table 2).Fu et al. [78] considered GP to be closely related to agricultural production conditions, artificial inputs and technical progress in the Loess Plateau.Our study showed that the increasing GP was mainly driven by the increasing agricultural investments in our site.Chemical fertilizer investment ranged from 155.71 kg/ha in 2000 to 330.29 kg/ha in 2013 (Figure 11a); it thus grew 1.12-times from 2000-2013 with a significantly increasing trend of 15.39 kg/ha•year (R 2 = 0.928, p < 0.01).Agricultural machinery power grew 2.07-times and ranged from 19.80 kW/ha in 2000 to 60.77 kW/ha in 2013, with a significantly increasing trend (R 2 = 0.996, p < 0.01) (Figure 11b).Moreover, our study also showed that chemical fertilizer input was more effective than agricultural machinery power input (Table 2).Other minor drivers of GP improvement in our site may include farm management and the construction of high quality basic croplands [2,65,66].According to the relationships between WY and different factors (Table 3), both T and cumulative afforested area had a weak connection with WY change from 2000-2013 (Table 3).However, the negative effect of afforestation was significant in the GFG-extensive period (Table 3).The response of WY change to NDVI change showed there was a negative relationship between vegetation restoration and WY (Figure 7a).In this study, WY was calculated by subtracting ET from PPT.During the study period, ET increased with the growth of woodlands (Figure 11c), which can partly explain the negative relationship between afforestation or vegetation restoration and WY.However, the increase in ET was slight, and the annual growth in PPT (8.46 mm/year) was much larger than the growth in ET (1.90 mm/year), which suggests that the slight improvement of WY was mainly driven by PPT during the study period.The strong and positive response of WY change to PPT change (Figure 8a), significant correlation (Table 3) and linear fitting (Figure 12a) between WY and PPT can enhance the credibility of that conclusion.However, the afforestation producing the negative effect on WY still cannot be ignored.According the relationships between NPP and different factors (Table 4), neither T nor PPT were significantly correlated with NPP during the study period, but the correlation between NPP and cumulative afforested area was significant.Xie et al. [59] consider land use change caused by GFG to be the main factor in increasing NPP in the Loess Plateau.Driven by GFG policy (the input of labor and funding), the decreasing croplands and increasing woodlands (including scrublands) were the main types of land use change in our site between 2000 and 2013.That change brought about the recovery of natural vegetation systems, which played a major role in the increase in NPP.Additionally, NPP was in fact significantly increasing in the mid-eastern NSLP, including Jia County, Hengshan County, Mizhi County, Wubao County, Zizhou County, Suide County, Zichang County, Qingjian County, Yanchuan County and Yan'an County, which correspond with the main GFG regions.The significant correlation (Table 4) and linear fitting (Figure 12b) between NPP and cumulative afforested area showed that the sustained growth of NPP was mainly caused by afforestation (or GFG policy).
The drivers involved in ES changes in our site from 2000-2013 are complex.For different ES, there are different dominant factors, some of them natural and other anthropogenic.The change in WY is attributed to the PPT belonging to the climate factor.The improvement in GP can be explained by the increasing in chemical fertilizer and agricultural machinery power, and the increasing in NPP can be explained by afforestation.Actually, these drivers (i.e., chemical fertilizer, agricultural machinery power and afforestation) can be attributed to the investments in natural capital.Ecosystems (e.g., forest ecosystem, farmland ecosystem) have been seen as natural capital assets, with the potential to generate a stream of vital life-support services [79].The stable and persistent investments in natural capital are crucial to improve some ES.Recent research at a large scale [1] have found that the investments in natural capital starting in 2000 have significantly improved the key ES in China.On a relatively small scale, our analysis showed that the investments in the farmland system improved food provision and that the afforestation increased climate regulation service.

Sustainability of Ecological Restoration and ES
GFG policies increased GP and NPP in the NSLP, and those policies provide a natural foundation for maintaining ES.However, the ecosystem and socioeconomic system are interactive, and during this interaction, farmers' attitudes and strategies are key to sustainable developments in ecological restoration [80].Close cooperation between government and local stakeholders can maximize the effects of GFG projects [2].The benefits produced by ecological restoration happen on a large scale (e.g., regional and national), while the costs of ecological restoration are undertaken by the local farmers.When the incomes of local farmers are temporarily reduced by ecological restoration projects, they typically re-plant later in the restored regions [81].The conflicts caused by the differing scales of the costs and benefits have an impact on the sustainability of ecological restoration.After a dozen years of GFG in the NSLP, ecological restoration has shifted focus from ecological benefit to reconciling ecological and social benefits.However, according to the data derived from the National Poverty Alleviation Office in 2015, 44% of the counties in the GFG region are still poverty-stricken (i.e., Qingjian County, Zizhou County, Suide County, Mizhi County, Jia County, Wubao County, Hengshan County, Zichang County, Ansai County and Wuqi County).Achieving a win-win situation involving both ES and poverty alleviation is the primary objective for the sustainable development of ecological restoration.
Ecological restoration improves regulating services, but it may at the same time cause a drop in provision services and, thus, have a negative impact on local farmers.A few studies have found that there is a trade-off between provision services and regulating services [82][83][84].However, our case study in the NSLP shows that the GP does not decrease, which mainly benefiting from sustaining investments in agriculture.This suggests that configuring effective compensation mechanisms (e.g., technological and economic) would contribute to the coordinated growth of both provision services and regulating services.Moreover, following the laws of ecological succession by choosing native tree species and developing more ecologically-informed GFG programs are both conducive to achieving stable and sustainable ecological restoration.

Conclusions
As a link between ecosystems and human well-being, ES are receiving more and more attention.As the GFG project proceeds, food provision, water regulation and climate regulation services in the NSLP are changing and have raised broad attention.To better understand the change trends and driving factors of the ES in the NSLP from 2000-2013, this study quantified the changes in NDVI/GP/WY/NPP and analyzed the relationships between GP/WY/NPP and different factors during the study period.The study found that, from 2000-2013, vegetation restoration was ongoing, but with fluctuation under the comprehensive influence of climate and human activity.GP and NPP both had a significantly increasing trend during the study period at the whole region scale in our site.Spatially, the regions with significant NPP growth corresponded strongly to the GFG counties.The changes in WY occurred in two stages: decline (2000)(2001)(2002)(2003)(2004)(2005)(2006) and growth (2007-2013).There were different dominant factors (i.e., natural factors vs. anthropogenic factors) driving GP/WY/NPP change in our site from 2000-2013.The changes in WY were influenced mainly by PPT.The improvement in GP and NPP can be attributed to investments in natural capital (i.e., chemical fertilizer, agricultural machinery power and afforestation).Meanwhile, we also found that vegetation restoration produced positive effects on NPP, but negative effects on WY.According to our research, we offer the following recommendations to other landscapes experiencing similar policies.In the process of vegetation restoration, we should follow the laws of ecological succession and choose native tree species to avoid climate inadaptation.Configuring effective compensation mechanisms (e.g., technological and economic) in the nature conservation programs could relieve the losses of local people and would be conducive to achieving stable and sustainable ecological restoration.S1: Information of weather stations within and near the Northern Shaanxi Loess Plateau (NSLP).

Figure 1 .
Figure 1.(a) The 25 counties of the Northern Shaanxi Loess Plateau (NSLP); (b) the location of NSLP; and (c) the main streams and key features in the study area.(DEM, digital elevation model.)

Figure 1 .
Figure 1.(a) The 25 counties of the Northern Shaanxi Loess Plateau (NSLP); (b) the location of NSLP; and (c) the main streams and key features in the study area.(DEM, digital elevation model.) y and z.The calculation of the correlation coefficient Rxy (or Rxz or Rxz) is detailed in Xu[60].

Figure 3 .
Figure 3. Changes in the annual NDVI from 2000-2013 in our site.

Figure 4 .
Figure 4.The (a) spatial pattern, (b) change rate and (c) significance of NDVI changes in our site from 2000-2013.

Figure 3 .
Figure 3. Changes in the annual NDVI from 2000-2013 in our site.

Figure 3 .
Figure 3. Changes in the annual NDVI from 2000-2013 in our site.

Figure 4 .
Figure 4.The (a) spatial pattern, (b) change rate and (c) significance of NDVI changes in our site from 2000-2013.

Figure 4 .
Figure 4.The (a) spatial pattern, (b) change rate and (c) significance of NDVI changes in our site from 2000-2013.

Figure 5 .
Figure 5. Changes in annual (a) grain production (GP) and gross productivity, (b) total water yield (WY) and (c) total net primary production (NPP) in our site from 2000-2013.

Figure 6 .
Figure 6.Spatial pattern of (a) WY and (d) NPP, change rate of (b) WY and (e) NPP, significance of (c) WY change and (f) NPP change in our site from 2000-2013.

Figure 5 .
Figure 5. Changes in annual (a) grain production (GP) and gross productivity, (b) total water yield (WY) and (c) total net primary production (NPP) in our site from 2000-2013.

Sustainability 2017, 9 , 199 9 of 20 Figure 5 .
Figure 5. Changes in annual (a) grain production (GP) and gross productivity, (b) total water yield (WY) and (c) total net primary production (NPP) in our site from 2000-2013.4.3.2.WY ChangeTotal WY changed substantially during the 14-year period, ranging from 1723.61 m 3 /year in 2004 to 20,896.08 m 3 /year in 2013.The changes in WY can be divided into two stages: decline(2000)(2001)(2002)(2003)(2004)(2005)(2006) and growth (2007-2013) (Figure5b).Overall, WY in our site showed only a slight increase from 2000-2013, which was consistent with the PPT (Figure2b).Other research[2] showed WY decreased after GFG in the whole Loess Plateau from 2000-2008.The results in this study show that WY decreased 24.81 mm during the same period.Substantial regional variation in WY was detected under the influence of climate and vegetation coverage (Figure6a).From the spatial change patterns, we see that 90% of the area experienced an increase in WY (Figure6b).A decrease in WY was identified in the northwest NSLP, including Dingbian County, Jingbian County and Wuqi County.The annual growth rate of PPT was low in those areas.The change in 75.24% of the area was not significant, and 24.76% of the area experienced significant growth (p < 0.05), as seen in northern Yulin, including Shenmu County and Fugu County (Figure6c), where there was a high annual growth rate of PPT (>8 mm/year) (FigureS2).

Figure 6 .
Figure 6.Spatial pattern of (a) WY and (d) NPP, change rate of (b) WY and (e) NPP, significance of (c) WY change and (f) NPP change in our site from 2000-2013.

Figure 6 .
Figure 6.Spatial pattern of (a) WY and (d) NPP, change rate of (b) WY and (e) NPP, significance of (c) WY change and (f) NPP change in our site from 2000-2013.

Figure 7 .
Figure 7. Response of (a) WY change to NDVI change and (b) NPP change to NDVI change at the pixel scale in our site from 2000-2013.

Figure 7 .
Figure 7. Response of (a) WY change to NDVI change and (b) NPP change to NDVI change at the pixel scale in our site from 2000-2013.
correlation coefficient; PR: partial coefficient; * correlation or partial correlation is significant at the 0.01 level; *** correlation or partial correlation is significant at the 0.1 level.

Figure 8 .
Figure 8. Response of (a) WY change to PPT change, (b) NPP change to T change and (c) NPP change to PPT change at the pixel scale in our site from 2000-2013.

Figure 8 .
Figure 8. Response of (a) WY change to PPT change, (b) NPP change to T change and (c) NPP change to PPT change at the pixel scale in our site from 2000-2013.

Figure 9 .
Figure 9.Comparison between MOD17-NPP and simulated NPP in each county during the study period.

Figure 10 .
Figure 10.The response of ET to NPP in terms of (a) spatial distribution and (b) the change trend at the pixel scale in our site from 2000-2013.

Figure 9 .
Figure 9.Comparison between MOD17-NPP and simulated NPP in each county during the study period.

Figure 9 .
Figure 9.Comparison between MOD17-NPP and simulated NPP in each county during the study period.

Figure 10 .
Figure 10.The response of ET to NPP in terms of (a) spatial distribution and (b) the change trend at the pixel scale in our site from 2000-2013.

Figure 10 .
Figure 10.The response of ET to NPP in terms of (a) spatial distribution and (b) the change trend at the pixel scale in our site from 2000-2013.

Figure 11 .
Figure 11.Changes in annual (a,b) agricultural investment and (c) ET (evapotranspiration) in our site from 2000-2013.

Figure 12 .
Figure 12.The linear fitting (a) between WY and PPT and (b) between NPP and cumulative afforested area in our site from 2000-2013.

Figure 11 .
Figure 11.Changes in annual (a,b) agricultural investment and (c) ET (evapotranspiration) in our site from 2000-2013.

Figure 11 .
Figure 11.Changes in annual (a,b) agricultural investment and (c) ET (evapotranspiration) in our site from 2000-2013.

Figure 12 .
Figure 12.The linear fitting (a) between WY and PPT and (b) between NPP and cumulative afforested area in our site from 2000-2013.

Figure 12 .
Figure 12.The linear fitting (a) between WY and PPT and (b) between NPP and cumulative afforested area in our site from 2000-2013.

Supplementary Materials:
The following are available online at www.mdpi.com/2071-1050/9/2/199/s1, Figure S1: The land use types in the NSLP in 2000 (a) and 2010 (b), Figure S2: Spatial distribution and change trends in temperature and precipitation from 2000 to 2013 in the NSLP, Figure S3: Comparisons of simulated NPP and MOD17-NPP values in our site using (a) the average annual value from the 25 counties, and (b) the annual value from the 14-year period as a whole, Table

Table 1 .
[34] use transition matrix for the Northern Shaanxi Loess Plateau (NSLP) during 2000-2010 (area × 10 2 km 2 ).The land use dataset was interpreted by Landsat TM and was provided by the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC)[34].

Table 2 .
Correlation and partial correlation between GP and T (temperature)/PPT (precipitation)/CF (chemical fertilizer)/AMP (agricultural machinery power)/LF (labor force) during the study period in our site.
R: correlation coefficient; PR: partial coefficient; * correlation or partial correlation is significant at the 0.01 level; ** correlation or partial correlation is significant at the 0.05 level.

Table 3 .
Correlation and partial correlation between WY (water yield) and T/PPT/CAA (cumulative afforested area) during the study period over the whole regional scale.
R: correlation coefficient; PR: partial coefficient; * correlation or partial correlation is significant at the 0.01 level; *** correlation or partial correlation is significant at the 0.1 level.

Table 2 .
Correlation and partial correlation between GP and T (temperature)/PPT (precipitation)/CF (chemical fertilizer)/AMP (agricultural machinery power)/LF (labor force) during the study period in our site.
R: correlation coefficient; PR: partial coefficient; * correlation or partial correlation is significant at the 0.01 level; ** correlation or partial correlation is significant at the 0.05 level.

Table 4 .
Correlation and partial correlation between NPP (net primary production) and T/PPT/CAA (cumulative afforested area) during the study period over the whole regional scale.
R: correlation coefficient; PR: partial coefficient; * correlation or partial correlation is significant at the 0.01 level; ** correlation or partial correlation is significant at the 0.05 level.