The Grain for Green Program Intensiﬁes Trade-Offs between Ecosystem Services in Midwestern Shanxi, China

: Ecological engineering is a widely used strategy to address environmental degradation and enhance human well-being. A quantitative assessment of the impacts of ecological engineering on ecosystem services (ESs) is a prerequisite for designing inclusive and sustainable engineering pro-grams. In order to strengthen national ecological security, the Chinese government has implemented the world’s largest ecological project since 1999, the Grain for Green Program (GFGP). We used a professional model to evaluate the key ESs in Lvliang City. Scenario analysis was used to quantify the contribution of the GFGP to changes in ESs and the impacts of trade-offs/synergy. We used spatial regression to identify the main drivers of ES trade-offs. We found that: (1) From 2000 to 2018, the contribution rates of the GFGP to changes in carbon storage (CS), habitat quality (HQ), water yield (WY), and soil conservation (SC) were 140.92%, 155.59%, − 454.48%, and 92.96%, respectively. GFGP compensated for the negative impacts of external environmental pressure on CS and HQ, and signiﬁcantly improved CS, HQ, and SC, but at the expense of WY. (2) The GFGP promotes the synergistic development of CS, HQ, and SC, and also intensiﬁes the trade-off relationships between WY and CS, WY and HQ, and WY and SC. (3) Land use change and urbanization are signiﬁcantly positively correlated with the WY–CS, WY–HQ, and WY–SC trade-offs, while increases in NDVI helped alleviate these trade-offs. (4) Geographically weighted regression explained 90.8%, 94.2%, and 88.2% of the WY–CS, WY–HQ, and WY–SC trade-offs, respectively. We suggest that the ESs’ beneﬁts from the GFGP can be maximized by controlling the intensity of land use change, optimizing the development of urbanization, and improving the effectiveness of afforestation. This general method of quantifying the impact of ecological engineering on ESs can act as a reference for future ecological restoration plans and decision-making in China and across the world.


Introduction
Ecosystem services (ESs) refer to all the benefits that human beings obtain directly or indirectly from the natural ecosystem to meet and maintain their living needs [1,2]. The Millennium Ecosystem Assessment (MEA) divides ESs into four basic types, including regulating services (e.g., soil conservation and carbon storage), provisioning services (e.g., water and timber), supporting services (e.g., biodiversity conservation and nutrient cycling), and cultural services (e.g., forest recreation) [2,3]. Human development patterns over the last few centuries have detrimentally affected the health and resilience of natural ecosystems [3,4]. Declines in ESs have been observed at global and regional scales, and these declines pose a significant threat to human well-being [2,5]. Ecological engineering is a widely adopted countermeasure that attempts to mitigate the contradiction between human and eastern regions is relatively high, and human activities and industrial development are mainly concentrated in the southeastern plains; the western loess hilly regions have broken terrain, barren soil, and sparse vegetation [34,35]. Due to the low coverage rate of surface vegetation coupled with the landform type of prevalent ravines, the area has serious soil erosion and is typically an ecologically fragile area [34]. In recent years, because of the GFGP, the vegetation coverage rate in this area has increased significantly, the functions of various ecosystems such as climate regulation and soil conservation have improved significantly, and ESs have, accordingly, changed significantly [24,35].
Remote Sens. 2021, 13, x FOR PEER REVIEW 3 of 18 heat, and sufficient sunlight. The average annual temperature is between 0.4 °C and 12.2 °C, and the average annual precipitation is between 438 and 588 mm. The elevation of the study area ranges from 561 to 2806 m a.s.l., with high terrain in the middle of the study area and lower terrain on the edges (Figure 1). Vegetation cover in the mountains of the central and eastern regions is relatively high, and human activities and industrial development are mainly concentrated in the southeastern plains; the western loess hilly regions have broken terrain, barren soil, and sparse vegetation [34,35]. Due to the low coverage rate of surface vegetation coupled with the landform type of prevalent ravines, the area has serious soil erosion and is typically an ecologically fragile area [34]. In recent years, because of the GFGP, the vegetation coverage rate in this area has increased significantly, the functions of various ecosystems such as climate regulation and soil conservation have improved significantly, and ESs have, accordingly, changed significantly [24,35].

Data Sources and Descriptions
In this study, we used multi-source data products, such as land use, meteorology, soil, and digital elevation models, to evaluate ESs. Detailed descriptions and data sources are shown in Table 1. In ArcGIS 10.2, all data are converted to the same projected coordinate system (WGS_1984_UTM_Zone_49N), and the "Resample" tool is used to unify the raster data resolution to 30 m.

Data Sources and Descriptions
In this study, we used multi-source data products, such as land use, meteorology, soil, and digital elevation models, to evaluate ESs. Detailed descriptions and data sources are shown in Table 1. In ArcGIS 10.2, all data are converted to the same projected coordinate system (WGS_1984_UTM_Zone_49N), and the "Resample" tool is used to unify the raster data resolution to 30 m.

Quantifying Ecosystem Services
The InVEST model is used to quantify four key ESs: water yield (WY), soil conservation (SC), habitat quality (HQ), and carbon storage (CS). WY is calculated based on the difference between annual average precipitation and actual annual evapotranspiration [29,36]. SC refers to the erosion control ability of the ecosystem to prevent soil loss and the ability to store and maintain sediment [36,37]. The sediment delivery ratio module calculates soil conservation services based on the difference between potential (under extremely degraded conditions without vegetation cover) and actual (under current land cover and management conditions) soil loss [19,36,37]. HQ refers to the ability to provide resources and environmental conditions for the survival and development of species or populations, which depends on the abundance of natural resources [36,38]. The habitat quality module calculates the HQ according to the habitat suitability of each land use type, the impact distance and weight of threat factors, and the sensitivity of each land use type to threat factors [36,38]. Through previous studies [38][39][40][41], we determined the impact distance and weight of threat factors, the habitat suitability of each land use type, and the sensitivity parameters to each threat factor (Tables S1 and S2). The carbon module quantifies CS using previous local research on the carbon density of different land use types [42][43][44] (Table S3). To avoid the influence of abnormal climate fluctuations in a single year, we selected the average rainfall and temperature from 2000 to 2018 as the general results from the study area [45,46]. Table 2 provide greater detail on the process of assessment of each ES.
annual water yield for each grid cell; AET x : annual actual evapotranspiration for pixel x; Px: annual precipitation on pixel x; Biophysical coefficients of model input are shown in Table S3.

SC
InVEST model sediment delivery ratio module SC: soil conservation; R: rainfall erosion factor; K: soil erosion factor; LS: slope length and gradient factor; C: vegetation cover factor; P: support practice factor. R and K are calculated to refer to the method of Yang et al. [32] and Zhang et al. [37]. We assigned C and P values according to existing literature [17,36,39] (Table S3).

HQ
InVEST model habitat quality module HQ: habitat quality; H j : habitat suitability for habitat type j; D xj : degree of habitat degradation in pixel x that is in habitat type j; K: half-saturation constant; Z: default parameter of the normalized constant model.

CS
InVEST model carbon module CS = C a + C b + C s + C d CS: carbon storage; C a , C b , C s , and C d are carbon densities in aboveground biomass, belowground biomass, soil, and dead matter, respectively, for each land use type.

Calculation of Trade-Offs Between Ecosystem Services
Correlation analysis is an effective tool to identify relationships between pairs of ESs, with significant negative correlations representing trade-offs and positive correlations representing synergies [47]. The size of the Pearson correlation coefficient indicates the strength of the trade-off and synergy relationships [47]. Obviously, this method ignores the Remote Sens. 2021, 13, 3966 5 of 16 difference in the geographical space of the change rate of the ES trade-offs. The root mean squared error (RMSE) quantifies the average difference between the standard deviation of a single ES and the average ES's standard deviation [47,48]. The dispersion degree of the standard deviation of distance of average ESs is described, and reflects the difference in the geographical space of the change rate of the ES's trade-offs [49,50]. Therefore, this study uses the RMSE to quantify the trade-offs between ESs. To eliminate the influence of ES unit differences, we first standardize the value of each ES.
where ES i is the standardized value; ES i,obs is the raw value; and ES i,min and ES i,max are the minimum and maximum values of the i ESs, respectively. RMSE is calculated as follows: where ES is the expected value of n kinds of ESs. In two dimensions, RMSE represents the distance from the coordinate point to the diagonal, and the relative position of the coordinate point represents the relative benefit of a certain ecosystem service [49]. Lu et al. [48] and Luo et al. [50] provide detailed instructions and procedures for the calculation of such trade-offs.

Actual Land Use Changes and Scenarios
The local administrative department of the GFGP provided vector data for the implementation area of the GFGP in Lvliang City as of the end of 2018. We set up a scenario where the GFGP was not implemented and quantified the impact of the GFGP on regional ESs by comparing this alternative scenario with the actual scenario.
(1) Actual scenario: we evaluated ESs before (2000) and after (2018) the implementation of the GFGP based on actual land use. By comparing ESs in 2000 and 2018, we can understand actual changes of ESs under the implementation of the GFGP.
(2) Alternative scenarios where the GFGP was not implemented (2018S): this is a simulated scenario. We assume that during the period 2000-2018, the actual GFGP implementation area did not implement the GFGP; that is, the land use types remained, unchanged, at their state in 2000, while the land use types in other regions were consistent with actual changes. By comparing ESs in 2018S and 2018, we were able to quantify the impact of the GFGP on ESs.
Terrain fragmentation due to soil erosion is the main cause of ecological degradation in the Loess Plateau [51]. The design and implementation of the GFGP on the sub-watershed scale to carry out comprehensive control of soil erosion has achieved good results [19,52]. In addition, as a physical geographical unit, the sub-watershed scale can more accurately reflect biophysical characteristics [19]. Therefore, we obtained the average value of each ES at the sub-watershed scale through the Zonal Statistics tool in ArcGIS 10.2. At the sub-watershed scale, the impact of GFGP on ESs was quantified and the trade-offs among ESs and their influencing factors were analyzed.

Geographically Weighted Regression Model
Previous studies have confirmed that there are obvious geographical differences in ESs [53,54]. It is difficult for the classic global regression to reflect the differences in the relationship between ES trade-offs and influencing factors in geographic space, and not fully reflect actual local processes [55]. Geographically weighted regression (GWR) obtains local coefficients by minimizing residuals, taking into account differences in the spatial variation in the relationship between ES trade-offs and influencing factors, which improves the reliability of the model [56].
ESs and their trade-off relationships are affected by factors such as land use [57], climate [58], vegetation [19], geomorphology [59], and urbanization [60]. We selected eight influencing factors to include in our model: the dynamic degree of comprehensive land use change (LUD; Refer to Li et al. [61] method for calculation), annual average temperature (TEM), NDVI, annual average precipitation (PRE), percentage of construction land (CON), elevation (DEM), slope (SLO), and potential evapotranspiration (PET) (Table S4 provides detailed calculation methods or sources). To avoid the influence of multicollinearity, all factors were tested for multicollinearity in SPSS 21, and the factors with VIF greater than five were eliminated (Table S5). These preliminary analyses left us with LUD, NDVI, PRE, DEM, and CON as independent variables and the trade-off relationships between ESs as dependent variables for our GWR model. The lower the AIC c value of the model output, the more concise the model and the more reliable the regression estimation. The higher adjusted R 2 indicates a higher explanatory power and a better fit [56]. The mathematical expression of the model is as follows: where y is the dependent variable; (u i , v i ) is the spatial location of the i-th sample; β 0 (u i , v i ) is the intercept; p is the number of influencing factors; x ij represents the independent variables; β j (u i , v i ) is the estimated coefficient of the i-th sample for the j-th driving factors; and ε i is the error term.  Table 3). The area of farmland and forest under the 2018S scenario decreased by 5.21% and 10.92%, respectively, and the area of grassland and shrub land increased by 4.20% and 0.87%, respectively (Table 3). Our results show that the GFGP led to a decrease of 23.69% in the area of farmland, and an increase of 24.62%, 19.78%, and 3.64% in the area of forest, grassland, and shrub land, respectively (Table 3), and thus was the main driving force for the significant increase in regional vegetation cover.

Trade-Offs Between ESs
The correlation between changes in ESs from 2000 to 2018 and 2018S (the scenarios if the GFGP were not implemented) was analyzed at the sub-watershed scale. CS, HQ, and SC have a significant synergistic relationship, and there is a significant trade-off between these ESs and WY (Table 4). In addition, the correlation coefficients (including positive and negative correlations) between paired ESs in the actual scenario are larger than those in the alternative scenario if the GFGP were not implemented (Table 4). This indicates that the GFGP has intensified the trade-offs and synergies between ESs.
We visualized the WY-CS, WY-HQ, and WY-SC trade-offs using root mean squared error (RMSE), and our results show that the west and southeast are the high value areas of the trade-offs (

Trade-Offs Between ESs
The correlation between changes in ESs from 2000 to 2018 and 2018S (the scenarios if the GFGP were not implemented) was analyzed at the sub-watershed scale. CS, HQ, and SC have a significant synergistic relationship, and there is a significant trade-off between these ESs and WY (Table 4). In addition, the correlation coefficients (including positive and negative correlations) between paired ESs in the actual scenario are larger than those in the alternative scenario if the GFGP were not implemented (Table 4). This indicates that the GFGP has intensified the trade-offs and synergies between ESs.    We visualized the WY-CS, WY-HQ, and WY-SC trade-offs using root mean squared error (RMSE), and our results show that the west and southeast are the high value areas of the trade-offs (Figure 4). Average tradeoff values of WY-CS, WY-HQ, and WY-SC are 0.051, 0.050, and 0.016, respectively, in the actual scenario, and the average tradeoff values of WY-CS, WY-HQ, and WY-SC are 0.028, 0.030, and 0.014, respectively, in the alternative scenario if the GFGP were not implemented (Figure 4). This indicates that the implementation of the GFGP strengthens the trade-offs between WY-CS, WY-HQ, and WY-SC.

FACTORS Influencing ESs Trade-Offs
We built a GWR model to explore the geospatial relationship between ES trade-offs and the factors that influence them. Compared with OLS, the adjusted R 2 of the GWR model was greater, and the AICc value decreased significantly ( Table 5), indicating that the GWR results have higher explanatory power and can more accurately reflect the processes at play.

FACTORS Influencing ESs Trade-Offs
We built a GWR model to explore the geospatial relationship between ES trade-offs and the factors that influence them. Compared with OLS, the adjusted R 2 of the GWR model was greater, and the AIC c value decreased significantly ( Table 5), indicating that the GWR results have higher explanatory power and can more accurately reflect the processes at play. The correlation coefficient of the GWR model reflects the spatial non-stationary response of ES trade-offs to influencing factors ( Figure 5). LUD is significantly positively correlated with WY-CS, WY-HQ, and WY-SC trade-offs, and the correlation coefficient is relatively high in the northeast ( Figure 5 (a1-a3), Table 6). NDVI is significantly negatively correlated with WY-CS, WY-HQ, and WY-SC trade-offs, and the correlation coefficient is relatively high in the southwest ( Figure 5 (b1-b3), Table 6). PRE is positively correlated with WY-CS and WY-HQ trade-offs, but negatively correlated with the WY-SC trade-off, and the correlation coefficient is high in the southeast ( Figure 5 (c1-c3), Table 6). DEM is negatively correlated with the WY-HQ trade-off, but positively correlated with WY-CS and WY-SC trade-offs, and the correlation coefficient is larger in the north ( Figure 5 (d1-d3), Table 6). CON is negatively correlated with the WY-CS trade-off, but significantly positively correlated with WY-HQ and WY-SC trade-offs, and the correlation coefficient is larger in the south (Figure 5 (e1-e3), Table 6).
negatively correlated with the WY-HQ trade-off, but positively correlated with WY-CS and WY-SC trade-offs, and the correlation coefficient is larger in the north ( Figure 5 (d1-d3), Table 6). CON is negatively correlated with the WY-CS trade-off, but significantly positively correlated with WY-HQ and WY-SC trade-offs, and the correlation coefficient is larger in the south ( Figure 5 (e1-e3), Table 6).   Local R 2 maps describe spatial differences in model goodness of fit, ranging between 0 and 1. Our results show that the selected five influencing factors are closely related to ES trade-offs, explaining 90.8%, 94.2%, and 88.2% of the WY-CS, WY-HQ, and WY-SC trade-offs, respectively ( Figure 6, Table 5). In general, LUD, CON, and NDVI are the most important driving factors of ES trade-offs, and they are significantly positively correlated with LUD and CON while being negatively correlated with NDVI ( Figure 6, Table 6). This means that increasing vegetation cover, controlling the intensity of land use change, and optimizing the development of urbanization are effective ways to alleviate the trade-offs between ESs and realize the synergistic promotion of multiple ESs.
ES trade-offs, explaining 90.8%, 94.2%, and 88.2% of the WY-CS, WY-HQ, and WY-SC trade-offs, respectively ( Figure 6, Table 5). In general, LUD, CON, and NDVI are the most important driving factors of ES trade-offs, and they are significantly positively correlated with LUD and CON while being negatively correlated with NDVI ( Figure 6, Table 6). This means that increasing vegetation cover, controlling the intensity of land use change, and optimizing the development of urbanization are effective ways to alleviate the trade-offs between ESs and realize the synergistic promotion of multiple ESs.

Effects of the GFGP on ESs
The GFGP is a successful program for coping with environmental degradation and increasing the supply of ESs [31]. Although the quality of the regional ecological environment has been greatly improved [17,62], realizing the coordinated development of multiple ESs is still a key consideration in optimizing GFGP policies. Our results show that the implementation of the GFGP significantly increased forest, grassland, and shrub land area (Table 3), and vegetation cover increased significantly [24], leading to significant increases in CS, HQ, and SC. This indicates that the GFGP promoted the synergistic relationship among CS, HQ, and SC. However, CS and HQ, under an alternative scenario where the GFGP was not implemented, significantly decreased (Figure 3), indicating that the GFGP effectively compensated for the negative impacts of external environmental pressures on CS and HQ. By comparing the actual scenario in 2018 with the alternative scenario, we found that the GFGP is the most important driving force for the increase in SC, with a contribution rate of 92.96% (Figure 3). In addition, our results show that the GFGP has had a significant negative impact on WY, with a contribution rate of −454.48% (Figure 3).

Effects of the GFGP on ESs
The GFGP is a successful program for coping with environmental degradation and increasing the supply of ESs [31]. Although the quality of the regional ecological environment has been greatly improved [17,62], realizing the coordinated development of multiple ESs is still a key consideration in optimizing GFGP policies. Our results show that the implementation of the GFGP significantly increased forest, grassland, and shrub land area (Table 3), and vegetation cover increased significantly [24], leading to significant increases in CS, HQ, and SC. This indicates that the GFGP promoted the synergistic relationship among CS, HQ, and SC. However, CS and HQ, under an alternative scenario where the GFGP was not implemented, significantly decreased (Figure 3), indicating that the GFGP effectively compensated for the negative impacts of external environmental pressures on CS and HQ. By comparing the actual scenario in 2018 with the alternative scenario, we found that the GFGP is the most important driving force for the increase in SC, with a contribution rate of 92.96% (Figure 3). In addition, our results show that the GFGP has had a significant negative impact on WY, with a contribution rate of −454.48% (Figure 3). This is mainly due to the large-scale planting of non-native vegetation, which leads to a significant increase in water consumption and evapotranspiration [27,63]. Our research confirms that the GFGP produces significant ecological benefits while also exacerbating regional water resource conflicts, which is consistent with previous studies [26,27,51,64]. However, unlike previous studies, we quantified the contribution rate of the GFGP to ES changes and the impact on the trade-offs/synergies between ESs, providing a more direct reference for alleviating regional water resource conflicts and realizing the synergistic promotion of multiple ESs.

Suggestions on the Inclusiveness and Sustainable Development of the GFGP
Identifying the dominant factors influencing the trade-offs between ESs is critical to formulating an inclusive and sustainable plan for the GFGP. There are obvious spatial differences in the relationships between ES trade-offs and their influencing factors [65,66], and classical global regression models did not fully reflect the relationships between the two in geographic space [59]. The local coefficients obtained by the GWR model by minimizing the residuals reflect the spatial non-stationary relationships between them [56], effectively overcoming the problems with classic regression models. We used the GWR model to explore the spatially non-stationary relationship between ES trade-offs and their influencing factors, and our results show that LUD, CON, and NDVI are the most important driving factors for ES trade-offs ( Figure 5, Table 6). The WY-CS, WY-HQ, and WY-SC trade-offs were significantly positively correlated with LUD and CON, but negatively correlated with NDVI ( Figure 5 (a1-a3, b1-b3, e1-e3), Table 6).
Land use change and urbanization are the main drivers of declines in CS, HQ, and SC [57,67], and also have negative impacts on the water conservation capacity of ecosystems [68]. However, land use change and urbanization also reduced the evapotranspiration of surface vegetation to a certain degree [69], and their impacts on precipitation at smaller timescales are also limited [70]. Therefore, increases in LUD and CON intensify the tradeoffs between WY-CS, WY-HQ, and WY-SC. NDVI is the most direct manifestation of the effectiveness of afforestation [71]. The GFGP is the main driver of the increase in regional NDVI [24], which not only improves CS, HQ, and SC but also improves the water conservation capacity of the ecosystem [72,73]. Therefore, increasing NDVI helps to alleviate the trade-offs between WY-CS, WY-HQ, and WY-SC.
The correlation coefficient of the GWR model reflects the spatial non-stationary response of ES trade-offs to their influencing factors ( Figure 5). In the northeast and south of the study area, urbanization developed rapidly and the intensity of human activity was high (Figure 2), so the correlation coefficient between LUD, CON, and ES trade-offs was relatively large ( Figure 5 (a1-a3, e1-e3)). In the southwestern region, the terrain is rugged and vegetation is relatively scarce [34], so the correlation coefficient between NDVI and ESs trade-offs is relatively high ( Figure 5 (b1-b3)). Therefore, controlling LUD and CON in the northeast and south, and increasing vegetation cover in the southwest, is essential to alleviate the WY-CS, WY-HQ, and WY-SC trade-offs.
In summary, we propose that future engineering projects should take into account the geospatial relationships between ES trade-offs and their influencing factors. By controlling the intensity of land use change, optimizing the development of urbanization, and improving the effectiveness of afforestation, the inclusive and sustainable development of the regional GFGP can be realized.

Uncertainties and Limitations
Our research provides a direct and flexible method to quantify the impacts of the GFGP on ESs, but it still has certain limitations. First, changing ESs is a complex process driven by factors such as nature, human activities, and climate change [17,74]. It is very difficult to completely quantify the impact of the GFGP on ESs. Our study used the average climate parameters from 2000 to 2018. Although this method is widely adopted [45,46], climate change during the research period will certainly have had an impact on ESs. Second, the input parameters of the model evaluation are taken from previous studies, but due to the limitations of our data sources, quality, and availability, we did not verify the results of the ES evaluations. Third, because of the limited availability of data, we had to ignore some details of the GFGP, such as tree species selection and configuration, vegetation management methods, etc., although these practices certainly could have a strong impact on ESs [71]. These problems may introduce some uncertainty into our model results. Therefore, it is necessary to obtain long-term positioning observation data and conduct more detailed research on the impacts of the GFGP.

Conclusions
Based on scenario analysis, we quantified the impacts of the GFGP on changes in ESs in Lvliang City, a typical ecologically fragile area, and analyzed the main forces driving ES trade-offs through spatial regression. Our research shows that the GFGP compensated for the negative impacts of external environmental pressures on CS and HQ, and significantly improved CS, HQ, and soil conservation (SC), but this improvement came at the expense of water yield (WY). While the GFGP promotes the synergistic development of CS, HQ, and SC, it also intensifies the trade-off relationships between these services and WY. Land use change and urbanization are significantly positively correlated with the trade-offs between WY-CS, WY-HQ, and WY-SC, while increasing NDVI helped to alleviate these trade-offs. Therefore, controlling the intensity of land use change, optimizing the development of urbanization, and improving the effectiveness of afforestation are effective ways to realize the inclusive and sustainable development of the GFGP. The general methods used in this study to quantify the impacts of ecological engineering on ESs can provide a reference for future ecological restoration plans and decision-making in China and around the world.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/rs13193966/s1, Table S1: Threats and their maximum distance of influence and weights. Table S2: The sensitivity of habitat types to each threat. Table S3: Table of biophysical coefficients  for InVEST. Table S4: Description of variables selected in this study.   Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Publicly available datasets were analyzed in this study. Data sources and access links are detailed in Table 1 and Table S4.