Shaanxi Province, located in the heart of northwestern China (Figure 1), lies in the middle reaches of the Yellow River and covers an area of approximately 205,624 square kilometers (Li and Li, 2022). The region experiences an average annual temperature of 9–16 °C and annual precipitation ranging from 340 to 1,240 mm. Spanning three climatic zones, Shaanxi exhibits diverse natural landscapes, including mountains, hills, basins, and plains, with precipitation decreasing from south to north. The northern mountains and the Qinling Mountains divide Shaanxi into three distinct natural regions: the Loess Plateau (northern Shaanxi), the Guanzhong Plain (central Shaanxi), and the Qinba Mountains (southern Shaanxi). These regions demonstrate significant disparities in economic development and characteristic spatial heterogeneity in land-use patterns (Li et al., 2020). Shaanxi’s ecosystems host diverse flora and fauna, playing an irreplaceable role in HQ maintenance, SC, and WY functions (Song et al., 2017). To mitigate environmental risks associated with land use and LUCC, it is critical to explore optimization strategies for ESs under multiple scenarios. Such efforts not only provide scientific decision-making support for regional sustainable development but also promote a harmonious balance between economic growth and ecological conservation.

This study mainly combines natural condition datasets, socioeconomic development datasets, and land use data to achieve multi-scenario simulation of land use and ESs evaluation. The specific data sources are introduced in Table 1.

Land use data (2017 and 2022) are subject to uniform quality control, including geographic correction, visual interpretation and field survey, to ensure that all land use types at all levels meet the accuracy requirements (more than 80%) (Liu et al., 2024). The primary land types are divided into Cropland, Forest, Grassland, Water, Constructed Land, and Bare Land. The National Compilation of Cost and Income of Major Agricultural Products is provided by the National Bureau of Statistics of China¹ to record the production costs and income of major agricultural products. The Shaanxi Statistical Yearbook is published by the Shaanxi Bureau of Statistics² to obtain crop planting area data. To ensure calculation stability and spatial data accuracy consistency, the land cover data are resampled to a resolution of 250 m $\times$ 250 m using the bilinear interpolation method (Shi et al., 2023).

Table 2 summarizes the abbreviations used throughout this study. This study follows an integrated technical workflow built upon a four-stage progression: driver identification, spatial constraint construction, scenario-based land-use simulation, and ecosystem service response and mechanism analysis (Figure 2). The methodological framework is structured as follows:

Fifteen natural and socio-economic variables were evaluated with the LightGBM model to quantify their contributions to land-use transitions. This ensures that the selected drivers are scientifically sound and suitable for subsequent simulations.

Key ESs (WY, SC, CS, HQ) were assessed with InVEST to derive Ecological Protection Importance (EPI), while Urban Development Suitability (UDS) was constructed to capture development potential. Integrating EPI and UDS yielded the Ecological Security Pattern (ESP). The Grey Multi-Objective Programming (GMOP) model then optimized future land-use demand under joint ecological and economic objectives.

Incorporating the validated drivers, ESP constraints, and GMOP-optimized demands, the PLUS model simulated four 2037 scenarios: Natural Development (ND), Economic Priority (EP), Ecological Land Protection (ELP), and Coordinated Development (CD). These simulations reveal how alternative policies reshape spatial land-use patterns.

ESs under each scenario were quantified using InVEST. Redundancy Analysis (RDA) was then used to identify dominant natural and socio-economic drivers of ES dynamics, clarifying the mechanisms underlying WSCH changes.

In this study, the LightGBM model was employed to predict multiple land-use categories, and the Receiver Operating Characteristic (ROC) curves were generated using observed land-use data (Figure 3). Model performance was systematically evaluated through the Area Under the Curve (AUC) and the Root Mean Square Error (RMSE). The results indicate that all land-use classes achieved AUC values above 0.85, and the RMSE values remained generally low, demonstrating the model’s strong discriminatory power and robust stability (Guo et al., 2021). Notably, Water and Bare land exhibited exceptionally low RMSE values (0.047 and 0.030, respectively), suggesting minimal prediction error. Although the RMSE for Cropland and Grassland was relatively higher, their AUC values (0.89 and 0.96, respectively) still reflect reliable classification performance overall (Guo et al., 2021). Feature importance analysis further revealed significant differences in driving factors across various land use types. The distribution of Cropland and Forest is primarily driven by natural environmental factors, with key explanatory variables being mean annual precipitation, vegetation cover, and elevation. Grassland distribution is closely associated with precipitation and vegetation cover, while the relatively high contribution of GDP suggests a potential link between Grassland distribution and economic activities, such as livestock farming or tourism. Water distribution is influenced by the combined effects of slope, elevation, and population density. Constructed Land is highly correlated with human activities, with nighttime light, population density, and economic development level (GDP) serving as major driving factors. Bare land distribution is mainly affected by topography (slope and elevation) and spatial distance to water sources, reflecting the combined influence of natural factors and human activities. Overall, natural factors predominantly govern the spatial distribution of land use types, while social factors significantly impact only the distribution of Constructed Land. These findings provide scientific evidence for optimizing regional land resource allocation and formulating sustainable development plans.

To validate the reliability of the model, this study selected nine natural factors and six social factors, and combined them with 2017 land use data to simulate the land use patterns of 2022 using the coupled LightGBM-PLUS model. The simulation results were then compared with the actual 2022 land use data for accuracy validation, employing the Kappa coefficient and Overall Accuracy (OA) for evaluation (Nie et al., 2023). The validation results show that the Overall Accuracy (OA) achieved 0.93, indicating a classification accuracy of 93% across the overall sample, with strong classification performance. The Kappa coefficient reached 0.90, further demonstrating the model’s reliability and high precision. Diagnostic results of the model indicate that the coupled LightGBM-PLUS model exhibits excellent simulation accuracy and classification performance, making it a robust tool for accurately predicting LUCC.

Figure 4 presents the spatial distribution of UDS and EPI. High UDS values are primarily located in the Guanzhong Plain, including Baoji, Xi’an, Xianyang, Weinan, and northern Yulin, indicating significant development potential and adaptability in these areas. Conversely, low UDS values are found in southern Shaanxi and parts of Yan’an, predominantly covered by Grassland and Forest. High EPI values are concentrated in the Qinling-Daba Mountains, characterized by their unique ecological environment and biodiversity (Wang et al., 2024), thus warranting priority protection. In contrast, low EPI values are mainly distributed in the Guanzhong Plain and northern Shaanxi, where ecological functions are relatively weak due to the impacts of agriculture and urbanization (Zhu et al., 2022). Ecological source areas, derived based on UDS and EPI, cover an area of 13424.81 km², primarily located in the Qinling-Daba Mountains in southern Shaanxi. This region serves as a core area for ecological protection due to its unique ecological environment and critical ecological functions. A 600-meter buffer zone was subsequently established around these areas to maintain landscape diversity. Using the MCR model and circuit theory, 45 ecological corridors with a total length of 537.62 km were identified. Considering the role of ecological corridors in enhancing ecological flows, a 200-meter buffer zone was established on both sides of the corridors, based on existing studies. Finally, the buffered ecological source areas and ecological corridors were overlaid to construct the ESP. As shown in Figure 4c, the constructed ESP covers an area of 21,784.88 km², accounting for 10.61% of the total study area. The ESP developed in this study effectively preserves the core areas of ESs while considering urban development potential.

Figure 5 shows the spatial distribution of land use and transition patterns under four scenarios in Shaanxi Province. The overall spatial structure remains consistent, with Forest concentrated in the Qinling–Daba Mountains, Grassland in northern Shaanxi, and Cropland in the Guanzhong Plain. In the baseline year (2022), Forest, Grassland, and Cropland accounted for 46.25, 25.24, and 25.42% of the total area, respectively.

The ND and EP scenarios reveal contrasting development trajectories. Under ND, LUCC is minimal: Cropland remains stable, Forest shows strong persistence, and Grassland transitions mainly follow natural succession processes. In contrast, EP accelerates land development, with Constructed land expanding most rapidly (6777.94 km²), largely at the cost of Forest and Grassland. Although Cropland continues to increase under strengthened protection policies, ecological pressure intensifies significantly due to large-scale conversion of ecological land to built-up areas.

In the ELP and CD scenarios, ecological protection emerges as the dominant driver. ELP promotes substantial ecological restoration, with extensive Cropland converted to Forest and Grassland, and strict limits placed on urban expansion. CD represents a balanced pathway: Cropland remains stable to ensure food security, Forest increases moderately, Grassland decline is alleviated, and Constructed land expands in a controlled manner (5,588.94 km²). Overall, EP imposes the greatest ecological pressure, ELP maximizes ecological gains, ND maintains historical trends, and CD achieves coordinated development between economic and ecological goals.

From 2022 to 2037, ESs in Shaanxi exhibit pronounced spatial differentiation, consistently showing higher values in the south and mountainous regions than in northern arid and plain areas (Figure 6). WY displays a southwest-to-northeast gradient driven by climatic and hydrological conditions. Although WY remains relatively stable across scenarios, ND and EP show slight declines, with EP reaching the lowest value due to intensified land development. In contrast, the CD scenario achieves the highest WY (1983.27 × 10⁶ mm), indicating that balancing development with ecological constraints effectively safeguards key hydrological functions.

SC is consistently highest in the Qinling Mountains, where terrain, vegetation, and humid climate jointly reduce erosion. Among the scenarios, ELP achieves the greatest SC improvement, reflecting the strong effectiveness of ecological protection in restoring soil quality. EP also shows localized SC enhancement linked to Grassland expansion, while CD achieves moderate gains beyond ND, demonstrating that coordinated land use remains beneficial for erosion control even without strict ecological prioritization.

CS and HQ show strong responses to land development intensity. High-value CS and HQ areas are concentrated in forest- and grassland-dominated regions, particularly the Qinling Mountains and the southern Loess Plateau. EP leads to the sharpest declines in both CS and HQ, underscoring the ecological costs of rapid development. Conversely, ELP delivers the greatest improvements through extensive ecological land restoration, while CD performs nearly as well, achieving substantial gains in CS (1932.31 Tg) and noticeable enhancements in HQ. Overall, the comparative scenario analysis reveals that EP imposes the greatest ecological pressure, ELP maximizes ecological benefits, and CD offers an optimal balance, demonstrating that integrated land-use policies aligned with ESP constraints can effectively enhance ecosystem multifunctionality.

Figure 7 illustrates the trends in the total ecological value and total economic benefits from 2022 to 2037 across different scenarios, revealing a significant trade-off relationship between ESs and economic benefits. In the ND scenario, both economic benefits and ecological value show an increase; however, the WY decreases, and the improvement of other ESs is limited, unable to cope with the escalating ecological pressure in the region. This indicates that the traditional land-use development model needs optimization to meet the new ecological and economic demands. In the EP scenario, economic benefits reach their highest value (41,917.51 × 10⁸ CNY), but ecological value drops to its lowest level. While this scenario maximizes economic benefits, it has a significant negative impact on ESs, such as a decrease in WY, CS, and HQ. This suggests that in the EP scenario, economic growth comes at the cost of partial loss in ecosystem functions, highlighting the risks of intensified land development and resource utilization leading to ecological degradation.

The ELP scenario demonstrates a significant increase in ecological value (4708.54 × 10⁸ CNY), making it the scenario with the highest ecological value among all, but its economic benefits are relatively low, amounting to only 34167.93 × 10⁸ CNY. It is the only scenario where total economic value experiences a decline. This outcome results from stricter ecological protection measures that effectively enhance the functions of ESs, but also significantly inhibit economic growth. This illustrates the deep conflict between ecological protection and economic development objectives in policymaking. In contrast, the CD scenario strikes a more balanced relationship between ecological value and economic benefits, with ecological value and economic benefits reaching 4668.02 × 10⁸ CNY and 35158.72 × 10⁸ CNY, respectively. By optimizing land use structure and rationally planning resource allocation, this scenario achieves steady improvements in ESs while maintaining a high level of economic benefits (35158.72 × 10⁸ CNY). The CD scenario embodies a sustainable development pathway, preserving critical ecological function zones while meeting economic growth targets, offering valuable insights for formulating future development strategies.

This study explores the response of ESs to LUCC in Shaanxi from multiple perspectives. Through literature review, machine learning validation, and accuracy analysis, 15 selected factors were identified as effective in explaining the driving mechanisms of LUCC. Among these, temperature, precipitation, NDVI, and GDP exert a significant influence on each land type. The final Kappa coefficient and OA accuracy of the model both exceed 90%, indicating a high level of reliability in the simulation results.

Under the ND scenario, land use pattern changes are minimal, with Cropland and Forest remaining stable. However, Grassland exhibits a significant transition to Cropland and Forest, reflecting the influence of natural driving factors under historical evolution trends. In the EP scenario, economic development demands drive the rapid expansion of Constructed land, leading to a significant reduction in Grassland and Forest areas, thereby exerting greater pressure on the ecological environment. Under the ELP scenario, the ecological benefits of Forests are well preserved, while urban expansion is restricted. Due to the implementation of the Grain for Green policy, the Cropland area has decreased. In Shaanxi, regions with high vegetation coverage are primarily concentrated in the Qinba Mountains, particularly in the Qinling region, which boasts superior ecological conditions. Under this scenario, forest vegetation in the Qinling region of southern Shaanxi is effectively protected and restored.

The spatial distributions of UDS and EPI reveal distinct regional disparities in urbanization potential and ecological conservation priorities within Shaanxi Province. The Guanzhong Plain, characterized by well-developed infrastructure and robust economic growth, exhibits high development potential. Conversely, areas in southern Shaanxi and parts of Yan’an are constrained by complex topography and limited transportation access, resulting in lower UDS values. These contrasts necessitate the formulation of differentiated development strategies. The Qinba Mountain region, as a national key ecological function area (Zhang and Liang, 2020), possesses significant ecological value but is also a contiguous poverty-stricken area, where balancing ecological protection and economic development becomes a core challenge (Li R. et al., 2024). By integrating UDS and EPI to construct the ESP, reasonable constraints can be established in ecologically suitable areas, avoiding unnecessary restrictions on urban construction, while formulating differentiated protection strategies based on the ecological characteristics of different regions, thereby achieving coordinated development of urbanization and ecological protection. The CD scenario using ESP aims to balance economic development with ecological protection by stabilizing Cropland, moderately expanding Constructed land, and rationally protecting ecological land to achieve coordinated development. The implementation of ESP strengthens ecological protection in key areas while allowing for modest growth in Constructed land to meet urbanization demands within the context of economic growth. This scenario not only optimizes the ecological protection pattern but also provides a pathway for regional sustainable development, ensuring that urbanization and ecological security progress in parallel, promoting coordinated development and the construction of ecological civilization.

The LUCC of each scenario reflects the impact of different goal-oriented approaches. The ND scenario emphasizes the continuation of historical trends, while the EP scenario focuses on economic development, leading to greater ecological pressure. The ELP scenario prioritizes ecological protection, significantly expanding ecological land. In the CD scenario, land development policies demonstrate strong adaptability and comprehensiveness, supporting economic growth while effectively maintaining the functions of ESs. This scenario provides an important reference for future land use planning in Shaanxi, particularly in terms of reinforcing ecological protection policies for key areas such as the Qinba Mountains, to achieve a win-win situation for economic, ecological, and social benefits.

The PCC analysis of WSCH shows significant synergies between WY, SC, CS, and HQ (Figure 8a). A healthy soil structure and vegetation cover help to enhance water source conservation capacity, reduce soil erosion, and simultaneously improve CS (Abdallah et al., 2021). The stable supply of water sources promotes plant growth, which in turn enhances carbon storage capacity. Soil conservation not only protects soil organic matter but also improves the long-term carbon storage potential, while increased vegetation cover helps to improve soil structure and enhance SC (Francaviglia et al., 2023). As an indicator of biodiversity, HQ has the strongest correlation with CS. Water, soil, carbon, and habitat form an interconnected ecological network, collectively influencing the overall function and health of the ecosystem.

The study of the spatiotemporal dynamics of ESs is a scientific basis for formulating ecological protection policies. However, previous studies have often been limited to the analysis of single factors and have not systematically revealed the changes under the combined effect of multiple factors (Liu et al., 2019). RDA as a multivariate statistical method, can effectively explore the relationships between multiple explanatory variables and response variables. This study explores the potential driving factors of ES changes in Shaanxi from 2022 to 2037 based on 15 selected variables, including 9 natural factors and 6 socio-economic factors. As shown in Figure 8b, the redundancy analysis (RDA) explains over 70% of the variance across the entire study area, with the majority of the driving factors exhibiting statistically significant effects (p < 0.001), indicating that these factors have a pronounced influence on changes in ESs (Shi et al., 2023). Among them, slope (r = 0.23 ~ 0.22, p < 0.01) and precipitation (r = 0.22 ~ 0.21, p < 0.01) were identified as the two factors that contributed the most to the overall changes in the four ESs from 2022 to 2037. The study found that in both the 2022 and 2037 CD scenario, slope is a key factor influencing SC, while precipitation has the most significant effect on WY. Figure 8b further reveals the negative impacts of nighttime lights, GDP, and population density on HQ and CS, with HQ showing a highly significant negative correlation with nighttime lights. In the 2037 CD scenario, the impact of NDVI and altitude on ESs increased compared to 2022. The direction of the arrows in the figure reflects the positive or negative effects of the driving factors on ESs, while the arrow length indicates the intensity of the effect, i.e., the contribution rate of the variable’s explanation.

The factor contribution results reveal the key roles played by topographical conditions and climate in regulating ESs (Figure 8c). Slope determines the spatial distribution of SC and WY. Steep slope areas are more prone to soil erosion due to gravity, while gentle slope areas are more conducive to the stable retention of soil and water (Han et al., 2019). Changes in precipitation not only directly affect the water yield capacity but also indirectly influence CS and HQ by regulating vegetation growth and surface runoff (Jia et al., 2024). Although the changes in Shaanxi’s ESs are mainly driven by natural factors, the negative impacts of socio-economic factors should not be overlooked. Future ecological protection must comprehensively consider both natural and human activity influences.

In future development, the spatiotemporal sequence analysis of ESs reveals the dynamic evolution characteristics and spatial differentiation patterns of four ES functions under different scenarios. The study indicates that between 2022 and 2037, WY, SC, CS, and HQ are driven by climate, topography, and LUCC, showing significant regional differences with the south performing better than the north, and mountainous areas outperforming plains. In the various scenarios, the ND scenario shows some improvement in both economic benefits and ecological value, but the enhancement of ESs is very limited. The limitations of the traditional land use development model make it difficult to address the increasing ecological pressures in the future (Long et al., 2021). In contrast, the EP scenario seeks to maximize economic benefits (41917.51 × 10⁸ CNY), but ecological value declines to its lowest point, with key ecological functions such as WY and CS significantly degraded, reflecting the ecological costs brought by rapid economic growth.

In the ELP scenario, ESs’ overall level is significantly improved, especially in SC and HQ. Through strict protection measures, this scenario achieves the highest ecological value (4708.54 × 10⁸ CNY), but economic benefits fall to 34167.93 × 10⁸ CNY, highlighting the suppressive effect of the ecological priority policy on economic growth. The CD scenario performs well across multiple ESs. By optimizing land use structure and reasonably allocating resources, it effectively balances ecological value and economic benefit, demonstrating a sustainable development path.

Recent studies on watershed dynamics, soil conditions and pollution patterns underscore that balancing social development with ecological conservation requires integrated, evidence-based land management. Morphometric analyses indicate that watershed fragility and erosion risks must be incorporated into regional planning (Shoumik et al., 2025), while advances in plant growth-promoting microbial technologies offer viable pathways to sustain agricultural productivity with reduced chemical inputs (Shoumik et al., 2025). In parallel, research on antibiotic removal, heavy-metal pollution control and saline-alkali soil remediation highlights the critical role of low-cost treatment technologies and soil-quality restoration in maintaining ecosystem functioning (Asghar et al., 2023; Ilić et al., 2024). Moving forward, a deeper understanding of the spatiotemporal evolution of ESs, combined with systematic monitoring, pollution mitigation and adaptive management, will be essential for aligning ecological enhancement with socioeconomic needs and for providing a scientific basis for differentiated ecological-protection policies and optimized land-use strategies.

The 15 driving factors selected in this study cover highly correlated natural factors and socio-economic development factors. The diagnostic results of these driving factors indicate that the selection is reasonable. However, the interactions between these factors were not fully considered, which may introduce some uncertainty in the simulation results. Furthermore, although the four development scenarios explored in this study provide a reference for future predictions, they do not cover all possible outcomes.

Future research should enhance ES assessments by integrating additional models and algorithms, such as the SWAT model (Akoko et al., 2021), while explicitly accounting for climate change impacts on LUCC to capture more comprehensive ES responses. Strengthening collaboration with regional land-planning authorities is crucial to incorporate local policy constraints and develop LUCC strategies that reflect practical development needs, thereby providing a scientific basis for CD scenario design and informing regional land-use and ecological protection policies. Moreover, uncertainties arising from inconsistent spatial resolution among datasets and the omission of future climate variability limit the robustness of current ES estimates (Zamora-Maldonado et al., 2025); integrating SSP-RCP climate scenarios into PLUS and InVEST simulations would improve temporal reliability and strengthen scenario-based projections.