We analyzed the impact of the Grain-to-Green Program (GTGP) on ecosystem service (ES) and interactions of ES pairs at diverse spatial-temporal scales.
•
The spatial pattern of ES changes, synergies, trade-offs, and bundles have spatial heterogeneity, with varied changes in different services across GTGP phases.
•
The GTGP improves soil erosion control and net ecosystem production, but also introduces new risks resulting from increased trade-offs.
•
We integrate the information of ES changes at multiple spatial-temporal scales into ES bundles and propose an environmental management strategy.
Abstract
Understanding the relationships between ecosystem services (ES) and the factors driving their changes over long periods and multiple scales is key for landscape managers in decision-making. However, the widespread implementation of restoration programs has led to significant ES changes, with trade-offs across space and time that have been little explored empirically, making it challenging to provide effective experience for managers. We quantified changes and interactions among five ES across various stages of the Grain-to-Green Program in the eastern Loess Plateau, examining these dynamics at threefold spatial scales. We observed notable increases in soil retention and Net Ecosystem Production but declines in habitat quality and Landscape aesthetics under afforestation. Over time, and with more integrated restoration strategies, synergies between ES pairs weakened, and non-correlations (even trade-offs) increased. To avoid unnecessary trade-offs, we recommend incorporating socio-ecological factors driving ES changes and ES bundles, informed by empirical experience, into proactive spatial planning and environmental management strategies for multi-ES objectives. The temporal lags and spatial trade-offs highlighted by this study offer crucial insights for large-scale restoration programs worldwide.
ES exhibit diverse types and spatial heterogeneity (Cord et al., 2017). Human activities, like grazing (Li et al., 2023), deforestation (Portela and Rademacher, 2001), urban agglomeration (Li et al., 2018; Zhang et al., 2017), continuously alter ES through land use and landscape changes. Meanwhile, ES interact through synergies or trade-offs (Bennett et al., 2009). A trade-off arises when the supply of one ES decreases due to the increase of another, while synergies emerge when multiple ES are enhanced simultaneously (Bennett et al., 2009; Rodriguez et al., 2006). Recent research focuses on incorporating ES assessments into urban and conservation planning by analyzing ES flows, trade-offs, and responses to climate and human activities within socio-ecological contexts. Numerous studies indicate that understanding the complex changes and relationships of ES, especially in ecologically fragile areas facing unique challenges, is essential for achieving ES enhancement goals through spatial planning (Dade et al., 2019; Davies et al., 2015; Wang et al., 2022).
In response to the extreme societal consequences of ecosystem degradation and climate change, the implementation of ecological restoration projects has increased globally, significantly affecting ES. Examples include the European Union's agri-environment measures (Gimona et al., 2023), China's land consolidation projects (Li et al., 2019), and the Three-North Shelter Forest Program (Zhai et al., 2023). Among these, the most notable is the Grain-to-Green Program (GTGP) (Fu et al., 2005; Xu et al., 2004). These long-term, large-scale ecological restorations have further complicated the already complex mechanisms of ES change. The GTGP, launched in 1999 in the ecologically fragile Loess Plateau, completed its first phase by 2014 and has been in its second ongoing phase since then (National Forestry and Grassland Administration, 2020). For over two decades, GTGP has effectively expanded China's vegetative cover (Zhao et al., 2019), contributing to erosion reduction (Jin et al., 2021) and forest carbon sequestration (Wang et al., 2017; Zhai et al., 2023). Nevertheless, its effectiveness remains debated (Xu et al., 2004). The enhancement of a single ES might lead to the weakening of others (Chen et al., 2010; Jia et al., 2017; Kong et al., 2023). The methodology of GTGP has evolved between its two phases. The first phase primarily focused on large-scale afforestation as an ecological restoration strategy to address the region's severe soil erosion. The second phase integrated a more diverse set of environmental management strategies to promote environmental welfare and sustainable development (Li et al., 2021), resulting in increased vegetation coverage and carbon storage (Wang et al., 2012). However, previous studies often did not distinguish between the phases of GTGP or analyzed only for specific years (Pu et al., 2023). Given the profound impact of GTGP as a large-scale restoration project, a more precise assessment of its phase-specific effects is necessary.
Bundle constitution, referring to the repeated spatiotemporal manifestation of trade-offs and synergies among ES, is critical for understanding ES dynamics. Methods like k-means (Raudsepp-Hearne et al., 2010) and self-organizing maps (SOM) (Dittrich et al., 2017) help elucidate the complex trade-off and synergistic relationships among ES and identify ES bundles. These researches reveal the internal connections among various ES and transcends static perspectives. Considering the dynamic interactions among multiple ES, shaped by changing ecological processes and escalating human activities, ES bundles capture both contemporary and historical data (Gou et al., 2021). This approach provides significant benefits by encompassing temporal and spatial dimensions. Consequently, due to their capacity to identify interrelations among multiple ES and encompass temporal evolution and spatial distribution, ES bundles are identified as effective management units. They are particularly suitable for constructing spatial planning and ecosystem management within a sustainable development framework incorporating multi-objective programming.
To bridge the gap between existing research and its application in environmental management, this study aims to uncover the long-term interactions of ES and their socio-ecological drivers. This is particularly relevant for ecologically fragile areas in China that require sustained ecological restoration. The study seeks to provide scientific evidence for ecological restoration projects and environmental management in such regions. The following scientific questions are posed: during different stages of the GTGP (1) where and how did changes and trade-offs/synergies of ES occur, (2) and what are the differences at the grid, sub-watershed, and city scales? (3) Did clusters of ES changes form units amenable to precise management? (4) How can the mechanisms of ES change inform spatial planning and management strategies?
2. Materials and methods
2.1. Study area
The study area, Shanxi Province, covering 156,700 km2 within the Yellow River watershed (110°14′–114°33′ E, 30°34′–40°44′ N), is located in the eastern Loess Plateau, China (Fig. 1a). It includes 11 prefecture-level cities and 12 sub-watersheds. Characterized by its semiarid climate and limited annual rainfall averaging 468.3 mm, it possesses a fragile ecosystem (Shanxi Climate Center, 2022). As of 2020, the area's population is approximately 34.91 million (Statistics Bureau of Shanxi Province, 2022). This region, with a long history of human activities, has experienced persistent soil erosion, exacerbated by climate change and intensified agricultural and urban development (Fu et al., 2017; Xiong et al., 2023). The Loess Plateau has been the focus of the GTGP since its inception in 1999, with full coverage achieved by 2002 (National Forestry and Grassland Administration, 2020). Initially aimed at mitigating severe soil erosion by converting sloped farmland into forest or grassland, these efforts significantly improved vegetation cover and ecosystem quality. Expanding the scope, the Chinese government, from 2014, introduced a new GTGP phase with diverse measures for the eastern Loess Plateau, encompassing mining restoration, pollution control, integrating fruit tree forests, and developing advanced farmland.
2.2. Data sources
This study utilized a combination of natural and socio-economic data. Natural data included surface temperature, soil texture, precipitation, evapotranspiration, land use, elevation, and third-level sub-watersheds. Socio-economic data encompassed population density data, yearbooks, and administrative boundary data, with specific details available in supporting information. Prior to analysis, all data were standardized to a spatial resolution of 1 km using bilinear interpolation in the UTM 49 N coordinate system. Computational tasks were conducted on the Google Earth Engine (GEE) platform using the “geemap” library (Wu, 2020). Raster operations and surface analysis, such as slope and aspect, were performed using the “Xarray” (Hoyer and Hamman, 2017) and “Xarray-spatial” (makepath, 2024) packages. Vector operations were carried out using “GeoPandas” (Jordahl et al., 2020).
2.3. Quantification of ES at threefold spatial scales
This study conducted a quantitative analysis of ES, focusing on four categories: provisioning, regulating, cultural, and supporting services, as classified by the Millennium Ecosystem Assessment (2005) of the United Nations. Based on previous studies (Cord et al., 2017; Cui et al., 2019; Hu et al., 2023; Wang et al., 2022) and insights from the study area, we specifically quantified five ES: water conservation (WC), soil retention (SR), habitat quality (HQ), net ecosystem production (NEP), and landscape aesthetics (LA). Detailed information regarding the selection process and computational methodologies can be found in the supporting information.
Research indicates that spatial variability is fundamental to ecosystem management (Xia et al., 2023). Consequently, this study employs three scales—grid, sub-watershed, and city—for analysis. We calculated the aforementioned five key ES and compared their states at three critical time points: 2001, 2014, and 2022. Average ES were computed at the grid, sub-watershed, and city scales to quantify their spatiotemporal variability. To facilitate comparisons across different years, ES were normalized on a scale from 0 to 1 (low to high) for all three years.
2.4. Assessment of relationships between ES
2.4.1. Correlation analysis between ES pairs
We conducted temporal and spatial correlation analyses to reveal the synergies and trade-offs among ES. The grid scale is crucial for reflecting spatial variability, and by calculating different ES within 6 km grids (Gong et al., 2023), we obtained the average ES at this scale. Based on annual grid-scale ES data, we performed pairwise Spearman correlation analyses for each grid, considering the data distribution characteristics across two stages (Lyu and Zhang, 2019). This process yielded a matrix of ES trade-offs and synergies with spatial distribution characteristics at the grid scale. At the sub-watershed scale, to mitigate the impact of size differences among sub-watersheds on ES change assessment, we quantified synergies by measuring the increase/decrease of pairwise ES. In the city scale analysis, after normality testing, we used Pearson correlation to analyze the relationship of ES changes, assessing the impact of different urban development strategies on ES. Spearman and Pearson correlations were computed using the “xskillscore” (Bell et al., 2021) and “SciPy” (Virtanen et al., 2020) packages, respectively.
2.4.2. Identify bundle
After elucidating the relationships of ES pairs across multiple scales, it's important to determine if these changes form units more suitable for unified management. Identifying clusters of various ES changes can help maximize management effectiveness. Bundle constitution, involving the repetitive spatiotemporal appearance of trade-offs and synergies of ES (Cord et al., 2017; Raudsepp-Hearne et al., 2010), reveals the intrinsic connections among different ES and enhances understanding of their supply, demand, and transport mechanisms. We employed SOM to identify ES bundles at the grid scale. SOM, an unsupervised learning neural network method, classified each grid into ES bundles based on the similarities of ES co-occurrence in space (Shen et al., 2021). The “MiniSom” (Vettigli, 2024) package was used to categorize ES at three stages. Transition matrices were then used to calculate the conversion and spatiotemporal distribution of each ES bundle.
2.5. Identification of the social-ecological drivers of ES
The formation, spatial distribution, and temporal evolution of ES bundles are affected by various drivers (Raudsepp-Hearne et al., 2010). Leveraging social-ecological drivers identified in previous research, we selected four categories encompassing a total of thirteen independent variables (Table S2) to explore the determinants of ES changes (Gong et al., 2023). Natural factors included daytime temperature (TD), nighttime temperature (TN), precipitation (PREC), and Normalized Difference Vegetation Index (NDVI). Landscape patterns comprised built-up ratio (BR), crop ratio (CR), forest ratio (FR), grassland ratio (GR), largest patch index (LPI), contagion (CON), and number of patches (NP). Human activities and economic development factors included population (POP) and gross domestic product (GDP). Additionally, at the grid scale, geographical conditions were also considered, specifically sub-watershed (SW) and city (CI). The GeoDetector method is a widely used approach for exploring spatial driving factors (Chen et al., 2020; Lyu et al., 2023). In this study, we used factor detection to determine how variations in different social-ecological factors influence changes in ES. The correlations were analyzed by constructing a statistic, which ranges from 0 to 1. The value reflects the extent to which the independent variables account for the dependent variables. To enhance the robustness of our model, we employed optimal discretization for stratifying independent variables. Both optimal discretization and factor detection were conducted using the “GD” package (Song et al., 2018) in R. The statistic was defined as follows:where denotes the levels of independent variables; and represent the number of grid cells in the level and the entire area, respectively; and represent the variance of a specific ES change in the level and the entire area, respectively. The statistic quantifies the driving forces of independent variable on ES. The larger the value, the stronger the influence of independent variable on ES, with implying that independent variable has no effect on ES and indicating that independent variable entirely controls ES.
3. Results
3.1. Spatiotemporal characteristics of ES in the east of Loess Plateau
3.1.1. Spatiotemporal pattern of ES
The analysis revealed that while the distribution of ES in the eastern Loess Plateau is relatively stable, it exhibits spatial heterogeneity (Fig. 2). Overall, WC, SR, and NEP have strengthened, whereas HQ and LA have weakened. These trends display similar spatial distribution patterns at both the grid (Fig. 2a) and sub-watershed scales (Fig. 2b). Enhancements in WC and SR initially appeared in the southern regions (sub-watersheds I, II, III) and expanded to the central and western areas (sub-watersheds XI, XII) in the second phase. NEP showed the most significant increase among the five ES. The northwestern Yellow River coastal area (sub-watersheds IV, XI, XII, X) had weaker NEP compared to other regions, but it improved almost universally in both phases. HQ and LA both exhibited spatial changes from northwest to southeast, with noticeable differences between sub-watersheds IV, XI, XII, X and others. HQ strengthened in the western and northern areas but was weaker in the central and southern lower-elevation regions. LA was weaker in the western and northern areas, improving in the central and southern regions. The five ES also demonstrated variations in cities (Fig. 2c). For instance, LA consistently showed a distribution from weak in the north to strong in the south across 2001, 2014, and 2022.
3.1.2. Characteristics of ES changes
We quantified the specific changes in different ES across two phases at threefold spatial scales. At the grid scale (Fig. 3a, Table S3), four ES, excluding HQ, were promoted in most regions during the first stage, while only two ES showed improvement in most regions during the second stage. Notably, NEP was enhanced across the entire study area in the first stage. In the second stage, only SR saw larger-scale enhancement, covering 94.66% of the entire area, compared to the first stage. The areas where LA and WC weakened were more than 20% greater than in the first stage.
At the sub-watershed scale (Fig. 3b, Table S4), although NEP was enhanced in all areas, 10 sub-watersheds showed less enhancement in the second stage compared to the first. Almost all sub-watersheds experienced a weakening of HQ, but the extent of this reduction was smaller in 10 sub-watersheds during the second phase compared to the first. SR was strengthened in almost all sub-watersheds, with 8 sub-watersheds exhibiting greater enhancement in the second phase than the first. Fig. 3b illustrates the rich spatiotemporal variations of different ES across the two phases in specific sub-watersheds.
Due to the smaller analysis units compared to sub-watersheds, the city scale reveals richer variability (Fig. 3c, Table S5). NEP was the most noticeably improved ES in most cities, while HQ showed more declining trends than at the sub-watershed scale. We selected three representative cities from north to south based on latitude: Datong (northern), Taiyuan (central), and Yuncheng (southern). A key observation is that NEP in Datong and Taiyuan was enhanced in both stages. In Datong, WC saw improvements in more areas, while LA experienced weakening in more regions. In Taiyuan, all ES had more enhanced than weakened areas in the second stage. Areas of enhancement for all ES, except LA, expanded in the second stage compared to the first. In Yuncheng, SR, NEP, and WC improved throughout the area in the first stage. However, in the second stage, SR and WC weakened in Yuncheng's more areas. Like the other cities, LA's weakening occurred more extensively in the second stage. More details can be found in Fig. 3c.
Combining insights from Fig. 2, we observe that NEP has almost universally improved, but the extent and magnitude of this enhancement in the second phase were less pronounced than in the first. SR has strengthened across the entire study area, with more extensive optimization in the second phase. However, specific local regions (e.g., Yuncheng, sub-watershed I) still face risks. WC experienced significant temporal and spatial variations. It enhanced in the southern region during the first phase, but weakened in many central and southern areas and strengthened in some northwestern regions in the second phase (e.g., Linfen and Yuncheng, sub-watersheds I, II, III). HQ exhibited weakening in both phases, particularly in the central and southern regions (e.g., Taiyuan, Jinzhong and Yuncheng), with the degree of weakening being lesser in the second phase. LA saw more extensive weakening in the second phase compared to the first.
3.2. Trade-offs and synergies between ES pairs at different scales
At the grid scale, pairwise correlation analysis of ES over consecutive years highlights changes in their relationships during two stages of GTGP (Fig. 4a, Table S6). Initially, three ES pairs showed trade-offs, while two exhibited synergies. In the later stage, only three pairs displayed either synergistic or trade-off relationships. A notable finding was the consistent synergy between SR and NEP across both stages, suggesting concurrent growth, as inferred from Fig. 2. Conversely, LA and NEP consistently demonstrated a trade-off. Many ES relationships evolved between the stages. In the second phase, the SR-WC synergy expanded, especially in the southeastern region. However, the SR-NEP synergy decreased in areas comprising 12% of the total. The correlation and trade-off coverage between HQ-NEP and SR-HQ reduced, and the LA-NEP and LA-SR synergy weakened, with an intensified trade-off, shifting from synergy to trade-off in some areas. Consequently, most ES pairs tended towards becoming non-correlated or more trade-off oriented in the second stage.
At the sub-watershed scale (Fig. 4b), the most prominent relationship among all ES pairs is the trade-off between NEP and HQ, observed in over 80% of the sub-watersheds. This relationship holds for both stages, where an enhancement in NEP corresponds with a reduction in HQ. Furthermore, the relationships among ES vary based on the geographical locations of the sub-watersheds, with pairs like WC-NEP, WC-LA, WC-SR, and NEP-SR showing distinct north-south differences. In the first stage, WC-HQ tended towards a synergistic decrease in northern sub-watersheds like X, while displaying a predominantly trade-off relationship in southern watersheds such as I. In the second stage, WC-SR indicated a synergistic increase in northern watersheds like VI and XII, but a synergistic decrease in southern watersheds like I and II. No sub-watershed exhibited consistent ES relationships across both stages, indicating complex and diverse changes. Taking sub-watershed IV as an example, the first stage showed trade-offs between WC-RS and NEP-RS and a synergistic increase for WC-NEP. The second stage, however, revealed trade-offs among WC-NEP, NEP-HQ, NEP-LA, and a synergistic decrease for WC-LA. Therefore, the variation in ES pair relationships at the sub-watershed scale is more diverse and complex than at the grid scale.
Among the cities (Fig. 4c), the changes of trade-offs and synergies between ES pairs were quite varied across the two stages of GTGP. However, larger cities showed a more stable situation. Taiyuan, the largest city in the eastern part of the Loess Plateau, exhibited the smallest variation in ES correlations across the two stages, with six pairs of ES consistently correlated. WC-SR and SR-HQ showed strong synergies in both phases indicating a relatively stable ES state in Taiyuan. Additionally, different cities located in the same or adjacent watershed show some similarities due to the similarity of geographical environment. For example, SR and WC showed a significant positive correlation in both Linfen and Lüliang, similar to their synergy observed in sub-watershed XI.
Although there was a general reduction in the correlation or an expansion of trade-offs between NEP and other ES, this trend was not uniform. In northern sub-watersheds, a decrease in WC correlated with an increase in NEP in the first stage, but in the second stage, they transitioned to a synergistic increase. In southern sub-watersheds, an increase in NEP coincided with a decrease in SR in the first stage, but in the second stage, they shifted to a synergistic increase (Fig. 4). Taiyuan and Yuncheng showed trends of synergistic weakening and increasing trade-offs between NEP and other ES, aligning with findings at the grid scale. However, in Datong, the synergy between NEP and LA increased in the second stage compared to the first.
3.3. Spatial-temporal patterns of ES bundles
We identified five ES bundles (B1–B5) at the grid scale and examined their distribution in 2001, 2014, and 2022. The main characteristics of these ES bundles are presented in Table 1 and Fig. 5.
Table 1. The area and mean altitude of different ES bundles.
Bundle
Name
Distribution
Altitude (m)
B1
Integrated ecological bundle
•
Across the eastern Loess Plateau
•
With a concentrated presence in areas of higher elevation, except for the northern bank of the Yellow River
1383.1
B2
Ecological Transition Bundle
•
In central and northern regions
1134.6
B3
NEP-LA Synergy Bundle
•
Forming north-south oriented strips
•
predominantly in the central and eastern
•
Mainly floodplains and basins
904.9
B4
HQ-WC Synergy Bundle
•
Mountainous region
1264.0
B5
LA Bundle
•
In the river valleys and basins
•
Near urban outskirts
723.2
B1 not only covers the largest land area but also provides the most extensive range of ES (Fig. 5, Fig. 6). The spatial distribution of B1 remained relatively stable in 2001, 2014, and 2022, with minor changes in its covered area (Fig. 5b). Prior to 2014, B1 experienced only a 2.4% reduction overall. Its most significant change was an 8.8% conversion to B3, a pattern consistent with Fig. 2a and b, reflecting the impact of the first phase of GTGP on enhancing NEP. Post-2014, B1 expanded by 4.7%, primarily due to assimilation from B4 in the northeastern part of the study area (accounting for 6.0% of B1's total area). This shift further accentuated the characteristic concentrated distribution of B1.
B2 though smaller in coverage, provides relatively balanced support for all five ES (Fig. 5). It is almost absent in the southern areas (sub-watersheds I, II, Fig. 6). Despite expanding by 12.6% from 2001 in the first phase, B2 shrank by 33.2% in the second phase. As a transitional zone between B1 and B3, B2 offers more HQ but less of the other ES compared to B3 (Fig. 5a).
B3 provides strong support for LA and NEP. Its coverage area is second only to B1. Throughout both phases of GTGP implementation, B3 has continuously expanded (Fig. 5), growing by 1.4 times by 2022. This expansion mainly resulted from areas transitioning to B1.
B4 is predominantly located in the high-altitude areas between the Yellow River and the Lüliang Mountains, with a few occurrences in the Taihang Mountains. After the first phase of the GTGP, B4 significantly shrank, particularly in the mountainous areas on both sides of the central Fen River valley (e.g., the eastern slope of the southern section of the Lüliang Mountains and the Zhongtiao Mountains in the southern part of the study area). The areas reduced in B4 largely transformed into B1 and B2, reducing the fragmentation of the bundle distribution and promoting a more contiguous pattern for B1. This aligns with the GTGP's phase one focus on converting sloped farmland into forest, reflecting the large-scale afforestation's improvement of ES like SR. In the second phase, B4 experienced a slight expansion, yet its coverage in 2022 only reached 75.8% of its extent in 2001.
B5 provides more LA than any other cluster. It saw little change over the two stages, but a slight expansion occurred in the central Fen River valley. B5 likely provides a higher level of LA due to its rich landscape types and complex patterns, but it offers the lowest support for HQ and SR among the five ES.
3.4. Social-ecological drivers of ES across time and scales
In the GeoDetector result at grid scale (Fig. S1), natural factors, especially PREC, were found to have a greater impact on WC and SR than socio-economic factors. PREC's influence on WC was notably stronger in the second phase of the GTGP (q = 0.86). Despite the initial phase's focus on forest restoration, the FR showed limited explanatory power for WC and SR changes. Geographical conditions like SW and CI also significantly affected these changes, though their impact diminished in the second phase. At the sub-watershed scale, PREC remained a primary factor, but the effects of landscape patterns, human activities, and economic development varied, requiring attention. For instance, in northern sub-watersheds during the second phase, POP and CR significantly influenced WC and SR changes in specific sub-watersheds. At the grid scale, CI and SW changes were primarily influenced by factors like the NDVI and TD, with NDVI being more explanatory for NEP changes in the first phase and TD in the second. POP's explanatory power for NEP changes increased significantly in the second phase, emphasizing the impact of human activities, a trend also observed at the sub-watershed scale, particularly in sub-watersheds III and V. Landscape compositions, such as GR and FR, varied in their impact on NEP changes across different watersheds. For HQ, landscape composition and pattern were key at the grid scale, with CR (q = 0.91) and GR (q = 0.57) being the most explanatory factors in the first phase. Their impact decreased in the second phase but remained significant. This suggests that afforestation may have limited effects on improving HQ compared to agricultural-forest landscapes. At the sub-watershed scale, CR was the most explanatory factor in most sub-watersheds, while GR was more important in watershed X. Landscape pattern also significantly influenced HQ changes in watershed XII. Regarding LA, while the explanatory power of various factors was generally weak at the grid scale, landscape pattern, indicated by LPI and CON, was important in specific sub-watersheds. CR had a stronger explanatory power for HQ changes across all sub-watersheds. Additionally, socio-economic factors in sub-watersheds IX and X were highlighted as important considerations for managers.
At the city scale, different factors had more specific and varied effects on ES changes. Notably, CR significantly impacted HQ (Fig. S3), a trend observed in most cities. Specifically, in Datong, PREC was the dominant factor affecting ES changes, especially for NEP, SR, and WC. GR consistently remained the key factor influencing HQ changes. In Taiyuan, the critical factors for SR and WC shifted from PREC to FR and POP, and for NEP, from TD to NDVI. However, CR was consistently the most significant factor affecting HQ changes, a trend also observed in Yuncheng. Moreover, in Yuncheng, WC changes shifted from being primarily influenced by PREC to more by POP, highlighting the importance of socio-economic activities on WC. Contrary to Taiyuan, Yuncheng's SR and NEP changes shifted from being mainly influenced by FR and NDVI to TD. These results suggest that urban ES changes are related to the timing and scale of afforestation under GTGP, indicating that human-led ecological restoration strategies indeed lead to ES changes. More details can be found in supplementary material (Fig. S3).
We also observed a general decrease in the explanatory power of various factors for all ES in the second phase across all threefold spatial scales. At the grid scale, this trend affected 69% of all ES-factor pairs, including the impact of SW on WC and SR, and NDVI on NEP. A similar pattern was noted at the sub-watershed scale. For instance, in the largest sub-watershed V, 60% of factors showed reduced explanatory power for different ES in the second phase. Some cities also exhibited the similar trend. In Yuncheng, the explanatory power of factors for all ES except WC notably decreased in the second phase. This pattern might reflect a diminished long-term impact of GTGP's first phase, which primarily focused on large-scale afforestation, on ES changes. The more integrated restoration strategies of the second phase led to more complex ES changes, but the degree of change caused to individual ES weakened. This suggests that as ecological restoration programs evolve and diversify, their impacts on specific ES become less direct and more nuanced, necessitating a broader and more adaptive management approach.
In summary, our study revealed that the dominant drivers influencing changes in various ES varied over time and across grid, sub-watershed, and city scales. While some ES had consistent dominant drivers across both stages (e.g., HQ, WC), others experienced shift in their primary drivers through different stages (e.g., NEP, LA). Notably, no single ES exhibited the same main influencing factors across all threefold spatial scales. Consequently, knowledge of multi-spatial scale attributions of ES changes should be integrated into landscape planning and environmental management strategies at various levels.
4. Discussion
4.1. Changes and risks in ES resulting from the long-term environmental management
Our study demonstrates that the GTGP has significantly influenced ES in the eastern Loess Plateau over the past two decades. ES in our study area exhibit distinct spatial heterogeneity and temporal changes, consistent with trends observed in previous research (Hu et al., 2021). Specifically, SR and NEP have notably increased, triggering diverse changes in other ES across different spatial scales. The region, known for its loess gullies and heavy rainstorms, has historically faced severe erosion, with escalating soil loss prior to 1999 (Borrelli et al., 2017). Soil erosion is primary cause for terrain fragmentation and ecological environment degradation (Fu et al., 2017). Therefore, GTGP initially aimed to control erosion and improve SR. Previous studies showed that average soil retention rate (the ratio of soil retention to potential soil loss in percentage) was up to 63.3% during 2000–2008 (Fu et al., 2011). Our study verifies that SR has improved in the eastern Loess Plateau through a different method. Indeed, SR improvements were observed at both grid and sub-watershed scales (Fig. 2). At the sub-watershed scale, the results (Fig. 2b) also show the similar trend and spatial distribution to other studies quantifying long-term SR changes (Hu et al., 2021). In eleven cities, across both phases, there was a greater increase in areas with enhanced SR compared to those with decreased SR (Fig. 3c). This indicates that large-scale afforestation, a key method of ecological restoration, positively impacted the initial target ES. However, these improvements displayed spatial heterogeneity, highlighting the need for environmental management authorities to focus more on regions with reduced SR, such as Yuncheng.
Despite carbon sequestration not being the primary objective of local ecological restoration and environmental management, NEP significantly improved across both stages and at nearly all spatial scales. This aligns with findings from numerous studies, which underscore the positive contribution of ecological restoration strategies to carbon sequestration (Feng et al., 2017; Lü et al., 2012). Forests, with their high carbon stock attributable to significant biomass (Feng et al., 2020), facilitated the overall enhancement of NEP due to extensive forest expansion in the first stage (Fig. 2a; 4a). In the eastern part of the Loess Plateau, 462,700 hm2 of land was restored or added to forest in the first stages of GTGP under the requirements of ecological restoration policy (Zheng, 2019). However, as ecological restoration goals and methods diversified in the second stage, NEP declined in 10.2% of the study area.
While effective in improving soil erosion and carbon sequestration, our findings indicate that LA and HQ are weakened under the GTGP. This could be a consequence of landscape homogenization and loss of distinct features due to large-scale afforestation and urbanization, such as rapid urbanization and road network growth (Zhang et al., 2018). The reduction in LA might lead to risks including diminished living experiences, inequity in infrastructure, and impacts on residents' health. A decrease in HQ directly affects biodiversity conservation, potentially escalating agricultural risks like reduced natural pollination and increased pests and diseases. The results showed that crop ratio was more important than forest ratio in Landscape patterns affecting HQ (Fig. S2, S3). Therefore, we argue that large-scale afforestation is not a long-term solution for resolving ecosystem issues and achieving sustainable development goals. Although initially effective for soil erosion mitigation, it may negatively impact other ES. Understanding the rich variations in ES under current strategies and the relationships between different ES is essential for making informed trade-offs and decisions.
4.2. Characteristics and spatial scale effects of ES interactions
Interaction analysis is crucial when integrating ES for spatial planning, environmental management, and decision-making. However, restoration programs, acting as experiments, often result in unintended consequences like spatial and temporal trade-offs. Our grid-scale study confirmed that large-scale ecological restoration programs in the eastern Loess Plateau primarily enhanced the synergy between SR and WC. However, substantial non-correlations and some trade-offs have emerged between NEP and other ES during the second phase. Previous research suggests that in areas with limited water resources and high biodiversity value, carbon sequestration is more likely to display potential trade-offs with other ES (Chisholm, 2010). Historically, the eastern Loess Plateau has faced issues of both resource-based and quality-based water scarcity. The per capita water resource availability here is only 281 m³, about one-sixth of the mainland China average (Fan, 2022). Moreover, nearly two decades of large-scale afforestation might have increased local biodiversity value but also intensified water resource consumption. We hypothesize that these factors led to the evolving trade-off between NEP and other ES, most notably the deepening trade-off between NEP and LA in the second phase (Fig. 4a).
Furthermore, the implementation of large-scale ecological restoration programs has induced multiscale consequences. Notably, in some eastern sub-watersheds, the cost of enhanced NEP was a significant reduction in WC, a trade-off that became apparent only in the second phase of GTGP. For instance, the NEP-WC relationship shifted from synergy in the first phase to trade-off in the second (Fig. 4b). Previous studies indicate that the expansion and growth of artificial forests are among the primary causes of reduced surface runoff and increased evapotranspiration (Li et al., 2021), particularly in dry and semi-arid climates. In the second stages of GTGP, these artificial forests reached their half-mature and mature stages, characterized by peak water consumption (Schwärzel et al., 2020). The maturing artificial forests in the revegetated subbasins enhanced NEP while reducing WC, thus expanding the trade-off between them. Our study confirms that FR and NDVI are among the main factors influencing NEP changes in sub-watershed III (Fig. S2), aligning with prior findings (Hu et al., 2021). In the most cities, similar to sub-watershed III, the NEP-WC trade-off relationship showed no significant change. The most noticeable relationship shift in Yuncheng was the transition from synergy to trade-off between NEP and SR (Fig. 4c). Additionally, the FR, which significantly impacted SR changes in the first phase, almost lost its influence in the second phase. This suggests that the impact of large-scale afforestation on SR is not long-lasting, and increasing NEP may pose greater risks of soil erosion.
In summary, the pathways influencing changes in ES are complex, heavily swayed by local factors particularly at sub-watershed and city scales. This underscores the need for a detailed spatiotemporal and cross-scale analysis to thoroughly understand ES dynamics. In the context of the Chinese government's strong emphasis on carbon neutrality goals, local managers at various administrative levels must carefully weigh the trade-offs between NEP and other ES. This cautious approach is crucial to prevent an excessive focus on enhancing carbon sequestration from causing significant declines in other vital ES.
4.3. Implications for targeted spatial planning and ecosystem management based on ES bundles and social-ecological drivers
Undoubtedly, accurately identifying drivers of ES is crucial for shaping corresponding spatial planning and ecosystem management strategies. These strategies aim to manage multiple ES concurrently, promoting sustainable development (Dittrich et al., 2017; Wang et al., 2024). We observed that the spatial distribution of ES bundles (Fig. 6) is similar to that of sub-watersheds (Fig. 1a) in eastern Loess Plateau. Therefore, it is efficient and appropriate to treat them both as units of environmental management. Based on the multi-scale analysis of the GTGP, we synthesized the dominant socio-ecological drivers of prominent ES within each ES bundle in 2022. Additionally, we analyzed the spatial-temporal changes in ES interactions and proposed spatial planning and management strategies at both grid and sub-watershed scales (Fig. 7).
Firstly, our analysis at the grid scale informed macro-level environmental management objectives based on the changes in ES and ES bundles. Continuously controlling soil erosion remains a formidable task in the eastern Loess Plateau, especially in areas B2, B3, and B5. Moreover, tailored strategies for different bundles are vital. Enhancing carbon sequestration in B2 and B4, and improving habitat quality in B3 and B5, are critical for achieving multi-ES management objectives. As spatiotemporal analysis revealed the degradation of most synergistic relationships, we proposed effective measures in these bundles to enhance ecosystem resilience and prevent new risks from escalating trade-offs (Fig. 4a). For instance, B1 should avoid trade-offs between NEP and HQ that could lead to ES decline and area reduction; B2 should promote synergy between HQ and LA, preventing trade-offs between WC and HQ, and aim to transform more areas into B1 instead of B3; B3 should avoid exacerbating trade-offs between SR and other ES; B4 should prevent intensifying trade-offs between SR and LA; B5 should promote synergies between NEP and LA, expanding the possibilities for achieving multiple ES objectives.
The sub-watershed scale, as a physical geographical unit, is adept at reflecting biophysical characteristics and offers a unique advantage in addressing the scale mismatch between ecological processes and human management (Feng et al., 2020). However, research into the application of this unit has been limited (Zhao et al., 2018). Our study addresses this gap, highlighting the spatial and temporal heterogeneity in ES relationships and their socio-ecological drivers at the sub-watershed scale. This provides decision-makers with more specific and detailed information for effective planning and management. The GTGP, as one of China's primary ecological restoration initiatives, has been integrated into Chinese spatial planning at multiple levels, from national to township (Liu and Zhou, 2021). This reflects the Chinese government's strong commitment to ecological protection and fulfilling its responsibilities to enhance ES for sustainable development. Although the integration of ecological restoration with spatial planning is comprehensive and advanced, its multi-tiered system still struggles to break through administrative boundaries to achieve multi-ES management. Our research confirms the necessity and feasibility of incorporating sub-watershed scale spatial planning and ecosystem management into the national territorial spatial planning framework and ecological restoration program. In China, investments in restoration, including GTGP, have exceeded USD 378.5 billion over the past decade (Li et al., 2021). Cooperative planning across administrative boundaries can optimize the benefits from these substantial management costs, enhancing the effectiveness of ecosystem conservation efforts (Guo et al., 2022).
5. Conclusion
We explored the spatiotemporal patterns of ES interactions and socio-ecological drivers across grid, sub-watershed, and city scales. Under the context of long-term implementation of large-scale ecological restoration projects, our findings support future planning and environmental management strategies across spatial scales. Firstly, the spatial pattern of ES changes was heterogeneous, with varied changes in different ES across GTGP phases. The current ecological restoration measures have achieved the first-phase goal of reducing soil erosion but have yet to fully enhance ecological quality. Secondly, synergies mainly occurred between SR and NEP, largely influenced by GTGP's initial large-scale afforestation activities, while trade-offs primarily existed between NEP and other ES, notably LA. Thirdly, ES management should focus not only on declining ES but also on pairs with reduced synergies (enhanced trade-offs) to avoid new threats from solely pursuing the improvement of a single ES. Understanding the relationships among multiple ES is crucial for future ecological restoration projects. Additionally, ES bundles and socio-ecological drivers of ES varied across spatial scales, indicating different management priorities and methods at various scales. If continuous human intervention in ecologically fragile areas cannot be avoided, we underscore the importance of considering spatial trade-offs and temporal accumulations in the implementation of large-scale ecological restoration projects. By integrating bundles and synergies, we delineated more optimal spatial management units. We propose spatial planning and management strategies based on our results to guide the sustainable development of regional ecosystems and enhance human well-being.
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgements
This research was funded by Natural Science Foundation of Shanxi Province, China, grant number 202303021222221. We would like to express our gratitude to Mr. Junpeng Liu and Mr. Baoshu Jia for their valuable suggestions on visualizing the research results.
S. Baumgärtner, C. Becker, K. Frank, B. Müller, M. Quaas
Relating the philosophy and practice of ecological economics: the role of concepts, models, and case studies in inter- and transdisciplinary sustainability research
P. Borrelli, D.A. Robinson, L.R. Fleischer, E. Lugato, C. Ballabio, C. Alewell, K. Meusburger, S. Modugno, B. Schütt, V. Ferro, V. Bagarello, K.V. Oost, L. Montanarella, P. Panagos
An assessment of the global impact of 21st century land use change on soil erosion
A.F. Cord, B. Bartkowski, M. Beckmann, A. Dittrich, K. Hermans-Neumann, A. Kaim, N. Lienhoop, K. Locher-Krause, J. Priess, C. Schröter-Schlaack, N. Schwarz, R. Seppelt, M. Strauch, T. Václavík, M. Volk
Towards systematic analyses of ecosystem service trade-offs and synergies: main concepts, methods and the road ahead
R. Costanza, R. d'Arge, R. de Groot, S. Farber, M. Grasso, B. Hannon, K. Limburg, S. Naeem, R.V. O'Neill, J. Paruelo, R.G. Raskin, P. Sutton, M. van den Belt
The value of the world's ecosystem services and natural capital
J. Förster, J. Barkmann, R. Fricke, S. Hotes, M. Kleyer, S. Kobbe, D. Kübler, C. Rumbaur, M. Siegmund-Schultze, R. Seppelt, J. Settele, J. Spangenberg, V. Tekken, T. Václavík, H. Wittmer
Assessing ecosystem services for informing land-use decisions: a problem-oriented approach
Assessing the soil erosion control service of ecosystems change in the Loess Plateau of China. Ecological complexity, special section: complexity of Coupled
J. Guo, C. Jiang, Y. Wang, J. Yang, W. Huang, Q. Gong, Y. Zhao, Z. Yang, W. Chen, H. Ren
Exploring ecosystem responses to coastal exploitation and identifying their spatial determinants: Re-orienting ecosystem conservation strategies for landscape management
K. Jordahl, J.V.D. Bossche, M. Fleischmann, J. Wasserman, J. McBride, J. Gerard, J. Tratner, M. Perry, A.G. Badaracco, C. Farmer, G.A. Hjelle, A.D. Snow, M. Cochran, S. Gillies, L. Culbertson, M. Bartos, N.Maxalbert Eubank, A. Bilogur, S. Rey, C. Ren, D. Arribas-Bel, L. Wasser, L.J. Wolf, M. Journois, J. Wilson, A. Greenhall, C. Holdgraf, Leblanc F. Filipe
Possibilities and requirements for introducing agri-environment measures in land consolidation projects in China, evidence from ecosystem services and farmers' attitudes
P. Virtanen, R. Gommers, T.E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S.J. van der Walt, M. Brett, J. Wilson, K.J. Millman, N. Mayorov, A.R.J. Nelson, E. Jones, R. Kern, E. Larson, C.J. Carey, İ. Polat, Y. Feng, E.W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E.A. Quintero, C.R. Harris, A.M. Archibald, A.H. Ribeiro, F. Pedregosa, P. van Mulbregt
SciPy 1.0: fundamental algorithms for scientific computing in Python
Understanding the effects of socio-ecological factors on trade-offs and synergies among ecosystem services to support urban sustainable management: a case study of Beijing, China
Impacts of urban expansion on ecosystem services in the Beijing-Tianjin-Hebei urban agglomeration, China: a scenario analysis based on the Shared Socioeconomic Pathways
With complex topography and geomorphology, mountainous cities possess abundant natural resources. They are constrained by ecological environment and topographic conditions, leading to a prominent contradiction between urbanization development and ecological protection. As a result, ecosystem services (ESs) are under greater regulatory pressure. The identification of ecosystem services bundles (ESBs) can be the foundation for developing zonal ecological protection planning policies. We took Chongqing as a case study, investigated the impact mechanisms of socio-ecological factors on the level of ES supply in each ESB. The findings reveal that: (1) The quantitative assessment of ESs for 2000, 2010, and 2020 showed that ESs were temporally stable and spatially heterogeneous. Areas with high supplies of food production (FP) and water yield (WY) were predominantly found in the northwestern cropland and urban built-up regions, whereas high supply areas for the other four ESs were primarily located in the northeastern Dabashan Mountains and the southern Wuling Mountains. (2) The quantification of trade-offs and synergies between ESs showed that FP had a trade-off effect with all five other ESs, while most other ES pairs exhibited synergistic effects. It was found that the interrelationships produced changes over time. (3) Then, three types of ESBs were identified. After examining the influence mechanisms of socio-ecological factors across the three ESBs, individual ESs were found to have essentially the same types of main impact factors in three ESBs, but varies in impact. (4) Finally, with reference to changes in ES levels and interrelationships and the driving mechanisms of socio-ecological factors in each zone, this study proposed zonal strategies for managing ecosystem services and optimizing territorial space based on the geographic characteristics and socio-economic development in different ESBs, with the goal of attaining sustainable urban development and improving human welfare.