Spatiotemporal patterns and driving mechanisms of interactions among cropland ecosystem services in the Loess Hilly–Gully Region: a case study of Yanchuan County, Shaanxi Province

The Loess Hilly-Gully Region has undergone profound land-use transitions under the coupled implementation of the Grain for Green Program and gully land consolidation. However, how these interventions have reshaped interactions among multiple cropland ecosystem services (ESs) and their driving mechanisms remains inadequately understood. Taking Yanchuan County as a case study, we quantified water yield (WY), soil conservation (SC), and food production (FP) for 2000, 2010, and 2020 using the InVEST model and statistical downscaling. Spearman correlation analysis and geographically weighted regression (GWR) were used to characterize pairwise ES relationships and their spatial heterogeneity, while a GeoShapley-based interpretable machine learning model was applied to identify the drivers of individual ESs and their pairwise interactions. All three services showed net increases, with the largest gains in the western valley systems where gully land consolidation was intensively implemented. The three service pairs exhibited distinct patterns: WY–SC had a stable moderate positive association across all periods, SC–FP remained weakly negative, and WY–FP showed no stable spatial association. GWR revealed pronounced spatial heterogeneity, with stronger positive WY–SC associations in the western valley systems. For individual services, temperature was the dominant negative driver of WY, slope was the leading positive driver of SC, and GDP was positively associated with WY and SC; attribution for FP was limited. For ES interaction intensities, the grouped geographic component dominated the GeoShapley magnitude across all pairs. Among substantive predictors, GDP ranked first for WY–SC, precipitation for WY–FP, and nighttime light intensity for SC–FP; the WY–SC attribution was interpreted cautiously because of limited spatial transferability. These findings indicate that the spatial reorganization of cropland under coordinated ecological restoration and gully land consolidation is associated with differentiated changes in ES interactions—stable positive associations for some service pairs and weak negative or negligible associations for others—providing scientific support for spatially targeted multifunctional cropland management in semi-arid loess landscapes.

The Loess Plateau is one of the most severely eroded regions in China and globally, characterized by long-standing ecological fragility and intense human–land conflicts (Hou et al., 2017; Wu et al., 2020). The large-scale Grain for Green program, initiated at the end of the last century, substantially restored vegetation cover and markedly reduced soil erosion, bringing about a historic improvement in the ecological condition of the Loess Hilly-Gully Region (Jin et al., 2021; Liu Y. et al., 2024; Xie et al., 2025). However, the extensive retirement of sloping cropland also reduced the agricultural production space and squeezed the traditional livelihood base of farming households (Wu et al., 2021; Ye et al., 2023). As a result, the tension between food security and ecological conservation has become increasingly salient (Liu Y. et al., 2024; Yang and Guo, 2024; Yang et al., 2025). To address this dilemma, the gully land consolidation project was implemented in representative areas such as Yan’an. Following the model of “returning farmland to forest on the ridges and consolidating gullies to create cropland on the valley floor,” this project converted wasteland, inefficient water bodies, and unused gully land into high-yield check-dam land and terraces. It compensated for the cropland lost to afforestation, optimized the spatial allocation of cropland resources, and substantially improved rural livelihoods (Li Y. et al., 2019; Qu et al., 2022; Guo et al., 2024). This drastic land-use transition has profoundly reshaped the structure and function of cropland ecosystems. Yet, the newly created gully-floor cropland and the remaining slope cropland together form a heterogeneous agricultural landscape. How do multiple ecosystem services interact across this landscape? What factors drive these interactions? These questions remain inadequately understood.

Cropland ecosystem services form a critical nexus connecting natural processes and human wellbeing (Ding et al., 2023; Gong et al., 2024). In the Loess Hilly-Gully Region provides both provisioning services—food and specialty agricultural products—and vital regulating and supporting services, including soil conservation, carbon sequestration, water retention, and habitat maintenance (Liu Y. et al., 2024). These services do not exist independently but exhibit complex trade-off and synergy relationships. For instance, maximizing grain production may impair soil conservation, whereas enhancing soil organic carbon stocks can synergistically improve soil fertility and water-holding capacity (Yang et al., 2026). Since the implementation of the Grain for Green program, many studies have assessed the values of individual services delivered by regional forests, grasslands, and croplands (Li G. et al., 2019; Song et al., 2024; Zhao et al., 2024) and examined the impacts of land-use change on ecosystem services (Giri et al., 2019; Chen et al., 2025; Yang et al., 2025). However, the new landscape pattern of “forest on the ridges, productive cropland on the valley floor” created by gully land consolidation may have fundamentally altered the interaction patterns, spatial differentiation, and temporal dynamics of cropland ESs compared with traditional slope-cropland systems (Qu et al., 2022; Huang et al., 2023). Current studies predominantly identify ES correlations from a global perspective, overlooking the spatial non-stationarity of driving factors, and thus fail to reveal how natural endowments, engineering interventions, and household management behaviors differentially modulate service interactions across different landscape positions.

Yanchuan County is located in the hinterland of the Loess Hilly-Gully Region in northern Shaanxi and serves as a pioneer and typical demonstration area for the gully land consolidation (Wang et al., 2019; Yurui et al., 2019; Zhang et al., 2020). Since 2000, the Grain for Green program and gully land consolidation have jointly driven drastic changes in cropland area, spatial configuration, and internal composition—slope cropland, check-dam land, and terraces—substantially reorganizing the agricultural production space (Qu et al., 2022; Yang et al., 2025). These spatial shifts provide an ideal natural laboratory for analyzing the spatiotemporal dynamics of cropland ecosystem service interactions. However, existing research on Yanchuan’s cropland ecosystem has largely focused on land-use/cover change or the optimization of individual cropping structures. The multifunctional evolution of cropland ESs, the spatiotemporal patterns of trade-offs and synergies among these services, and the underlying coupled natural–human driving mechanisms have yet to be systematically quantified. This knowledge gap leaves territorial spatial optimization and multifunctional cropland management without a solid scientific foundation.

In light of this, the present study takes Yanchuan County as a case study. We first characterize the spatiotemporal distributions of three key cropland ESs—water yield, soil conservation, and food production. We then integrate correlation analysis, geographically weighted regression (GWR), and Shapley-value-based interpretable machine learning (GeoShapley) to systematically examine the spatiotemporal patterns of ES interactions and to quantitatively unravel their nonlinear driving mechanisms. Multi-dimensional factors, including physical geography, land use, and human activities, are considered. This study aims to deepen the scientific understanding of multifunctional interactions within cropland ecosystems in the Loess Hilly-Gully Region and to provide a reference for optimizing cropland allocation and reconciling rural livelihoods with ecological conservation following gully land consolidation.

Study area overview

Yanchuan County is located in the northern part of Shaanxi Province, on the west bank of the Yellow River (109°36′20″–110°26′44″E, 36°37′15″–37°20′55″N; Figure 1). It covers a total area of approximately 1,985 km 2, with a registered population of 180,200 and a permanent population of 136,200 in 2024. Major crops include wheat, millet, maize, legumes, jujube, and apples. In 2024, the gross regional product reached 14.618 billion CNY, and per capita disposable income of rural residents reached 15,695 CNY. The county is characterized by a large agricultural population and relatively low rural income, with leading industries comprising apples, jujubes, greenhouse agriculture, and animal husbandry. Rural livelihoods remain heavily dependent on cropland resources.

Yanchuan County possesses notable typicality and representativeness for studying cropland ecosystems in the Loess Hilly-Gully Region. First, situated in the heartland of the Loess Plateau, it exhibits highly developed loess hilly-gully landforms, with a gully density reaching 4.68 km/km 2 and a high degree of surface dissection. The high proportion of slope cropland and severe soil erosion made it a key area for the Grain for Green program. Large-scale cropland retirement substantially reduced the cropland area, with per capita basic cropland dropping below 0.08 ha, intensifying the conflict between ecological protection and food security. Second, Yanchuan was among the earliest pilot counties for gully land consolidation, beginning systematic implementation in 2011 and expanding under the Yan’an gully land consolidation major project in 2013. By 2024, the county had completed a construction scale of 3,561.59 ha, generated 801.67 ha of new cropland, and built 2,795.03 ha of high-standard basic cropland, with grain yields on consolidated land increasing substantially. The coupled implementation of the Grain for Green program on ridges and gully land consolidation on valley floors has led to a drastic reorganization of cropland quantity and spatial configuration, forming a landscape in which retired forest land on ridges coexists with high-yield check-dam land and terraces in valley bottoms. This composite cropland ecosystem is widely representative across the Loess Plateau, and Yanchuan has accordingly been selected as a typical case for investigating land-use change and the evolution of newly created cropland (Wang et al., 2013; Xiang et al., 2023; Chen et al., 2024). However, existing studies have not systematically addressed how the multiple ESs of different cropland types change over time under such drastic spatial reorganization, how trade-off and synergy relationships among these services are reconfigured, and which natural and social factors dominate their evolving interactions. Exploring these questions is essential for understanding the coupling mechanisms between cropland multifunctionality and sustainable rural livelihoods following gully land consolidation on the Loess Plateau.

Cropland ecosystem service assessment methods

The Outline of the Yellow River Basin’s Ecological Protection and High-Quality Development Plan explicitly identifies water scarcity as the principal contradiction within the basin, stipulating that water resources must be treated as the most binding constraint and that water conservation and intensive utilization must be advanced across the entire basin. Concurrently, it calls for the integrated promotion of basin-wide ecological protection, management, and restoration, highlighting soil and water conservation in the middle reaches, strengthening forest and grassland protection and restoration, enhancing soil and water conservation capacity, and promoting the development of dryland farming tailored to local conditions. Yanchuan County, situated in the heartland of the Loess Plateau and the hilly-gully region in the middle reaches of the Yellow River, epitomizes the region where these contradictions and tasks are most concentrated. On the one hand, precipitation is scarce and highly concentrated; the water resource endowment is extremely poor, far below the internationally recognized threshold of extreme water scarcity. Agricultural production is thus highly dependent on the spatiotemporal matching of limited precipitation. On the other hand, the terrain is highly dissected, with dense gullies, and the sensitivity to soil erosion is extremely high; the soil conservation function of slope cropland directly determines the preservation of the plow layer and the control of sediment delivery into the Yellow River. Soil conservation and food production therefore constitute the two most central service functions of cropland ecosystems in the study area, and their trade-off relationship is particularly sensitive under the dual disturbance of the Grain for Green program and gully land consolidation.

At the same time, the Loess Plateau, as a vital component of China’s “Two Ecological Barriers and Three Shelterbelts” ecological security strategic framework, performs a significant ecological regulation function. The Yellow River Basin also hosts a concentration of national energy and heavy-chemical industry bases; the Shaanxi section encompasses nationally designated key ecological function zones and energy-rich areas. Under the extreme scarcity of water resources in the basin means that the partitioning of precipitation among evapotranspiration, runoff, and vegetation utilization profoundly influences the trade-off and synergy relationships among ecosystem services. Annual water yield thus becomes a core indicator connecting climatic drivers and ecological responses. In light of the aforementioned regional functional positioning, resource background characteristics, and data availability and reliability, this study selects three key ecosystem services—annual water yield, soil conservation, and food production—for quantitative analysis, encompassing the representativeness of provisioning, regulating, and supporting services. All three services are quantified spatially using the InVEST model and statistical methods. The calculation principles are as follows.

Water yield

Annual water yield (WY) quantifies the amount of water generated within a landscape that becomes available for downstream use, and serves as a core indicator of water supply capacity in the water-scarce Loess Hilly-Gully Region. We applied the InVEST 3.14 Annual Water Yield module (Li et al., 2021) to simulate gridded water yield using a water balance approach based on the Budyko framework (Zhang et al., 2004). For each pixel, WY is computed as precipitation minus actual evapotranspiration (Equation 1). The InVEST output is originally in millimeters; we converted it to volume by multiplying by the pixel area, then to mass assuming a water density of 1000 kg/m 3 (1 mm of WY equals 0.9 t for a 30 m pixel). County-level totals were obtained by summing the values of all cropland pixels within the year-specific cropland mask. where WYxj is the average annual water yield for raster cell x under land cover type j (mm); AETxj denotes the annual actual evapotranspiration for the same cell and cover type (mm); and Px is the annual precipitation at cell x (mm). The evapotranspiration ratio is obtained from Equation 2:

With Equations 3–5:

Here Rxj represents the Budyko dryness index for cell x and cover type j (potential evapotranspiration/precipitation); ωx is a dimensionless parameter reflecting the plant-available water fraction; Z is an empirical seasonality constant; AWCx is the volumetric plant available water capacity (mm); Ds and Dr denote the depth of the root-restricting layer and vegetation rooting depth (mm), respectively; and PAWCx is the plant available water capacity of the soil, derived from soil texture and organic matter content as Equation 6:

Where: sand%, silt%, and clay% represent the respective proportions of sand, silt, and clay within the soil composition. Additionally, OM% indicates the percentage of soil organic matter content.

Soil conservation

Soil conservation service refers to the capacity of an ecosystem to prevent soil loss induced by rainfall and to trap sediment originating from upslope contributing areas. This function is of fundamental importance in Yanchuan County, where the combination of highly dissected terrain (gully density 4.68 km/km 2), intense summer rainstorms, and extensive historical slope cultivation has produced severe erosion. The InVEST Sediment Delivery Ratio module, which is built upon the Universal Soil Loss Equation (USLE), was employed to map soil conservation across the study area (Hamel et al., 2015). The module distinguishes two additive components: erosion mitigation, defined as the difference between potential soil loss under bare conditions and actual soil loss under current vegetation and management, and upstream sediment retention, which accounts for the mass of eroded material that is intercepted before reaching waterways (Liu Q. et al., 2024). The calculation framework is presented in Equations 7–10: where SEDRETx is the total soil conservation at cell x (t·pixel -1 ·yr -1); RKLSx is the potential erosion on bare ground; USLEx and USLEy are the actual erosion at cell x and an upslope cell y, respectively (t·ha -1 ·yr -1 or t·pixel -1 ·yr -1); SEDRx is the retained sediment quantity; SE is the sediment retention efficiency; Rx is the rainfall erosivity factor (MJ.mm·ha -1 ·hr -1 ·yr -1); Kx is the soil erodibility factor (t.ha·hr·MJ -1 ·mm -1); LSx is the dimensionless slope length and steepness factor; Cx is the cover-management factor; and Px is the support practice factor.

Food production

Food production is the most direct provisioning service of cropland, and its spatial pattern in Yanchuan has been profoundly reorganized by the dual drivers of slope-land retirement and gully-floor cropland construction. In this study, Food production (FP) was defined as grain-crop production service intensity, including maize, rice, wheat, soybean, potato, millet, and sorghum. Livestock products and non-grain economic crops, such as apples and jujubes, were excluded. County-level grain production in each year (Py, t) was used as the control total, while the approximately 1-km harvested-area data of Dai et al. (2025) were used to characterize its spatial distribution. The harvested-area grids were intersected with the corresponding 30-m cropland layer. Because positive harvested-area cells covered only a small proportion of cropland, a smoothed source-grid weight was calculated as Equation 11: where H j is the summed harvested area of the seven grain-crop groups in source grid j. The grid weight was assigned to the intersecting 30-m cropland cells, and FP intensity was calculated as Equation 12: where FP i is FP intensity in 30-m cropland cell i (t ha −1), A i is its effective cropland area (ha), w i is the corresponding source-grid weight, and C y is the set of valid cropland cells in year y. This formulation ensures (Equation 13):

The original 30-m raster was used for mapping and county-total validation. For Spearman correlation, GWR, and GeoShapley analyses, FP was aggregated to the common 1-km analysis grid using the cropland-area-weighted mean (Equation 14): where FP g is the FP intensity of 1-km analysis cell g. The aggregated layer retained the unit t ha −1. Further processing and validation details are provided in Supplementary Methods S3.

Analysis of cropland ecosystem service interactions

To identify pairwise trade-offs and synergies among water yield (WY), soil conservation (SC), and food production (FP), Spearman’s rank correlation analysis was performed using service values aggregated to the common 1-km analysis grid for each observation year. A positive and statistically significant coefficient indicated a synergistic relationship, whereas a significant negative coefficient indicated a trade-off; a non-significant coefficient indicated no statistically detectable interaction.

The corrplot package in R was used to calculate Spearman’s ρ and the associated p-values (Yang et al., 2025). The correlation coefficient was calculated as Equation 15: where ρ ranges from −1 to 1, d i is the difference between the ranks of the i-th pair of observations, and n is the number of paired 1-km grid cells. A service pair was classified as synergistic when ρ > 0 and p < 0.05, as a trade-off when ρ < 0 and p < 0.05, and as non-significant otherwise.

To characterize the spatial heterogeneity of service interactions at the same 1-km grid scale, geographically weighted regression (GWR) was used to estimate a local coefficient for each service pair (Sannigrahi et al., 2020). SC was regressed on WY for the WY–SC pair, FP on SC for the SC–FP pair, and FP on WY for the WY–FP pair. The models were implemented using the GWmodel package in R (Yang et al., 2025) with an adaptive Gaussian kernel. Bandwidths were selected by minimizing the corrected Akaike information criterion, and local pseudo t-values were calculated to assess coefficient significance. The resulting coefficients represent spatially varying statistical associations at 1-km grid cells rather than causal effects.

Driving mechanisms of ecosystem service interactions

Selection of socio-ecological driving factors

Drawing on previous studies conducted in the Yellow River Basin and considering the specific conditions of cropland ecosystems in Yanchuan County, we grouped potential driving factors into three categories: climate, topography, and human activities. A total of twelve candidate indicators were initially compiled. Following a multicollinearity diagnosis using the variance inflation factor (VIF > 5) and pairwise correlation screening (|r| > 0.7), nine indicators with low redundancy were retained (Table 1). Climatic factors comprise mean annual temperature (TMP), mean annual precipitation (PRE), and mean annual evapotranspiration (PET), which collectively govern the hydrothermal conditions for crop growth and ecosystem processes. Topographic factors include average elevation (ELE) and average slope (SL); in the highly dissected loess landscape, these variables strongly influence the spatial redistribution of water, soil, and nutrients. Surface fragmentation, terrain relief amplitude, and road accessibility were included in the initial candidate set but were excluded during screening because of high collinearity with the retained variables. Human activity indicators are represented by population density (POP), gross domestic product (GDP), the nighttime lighting index (LACC), and the human footprint intensity index (HFP). These dimensions jointly capture the pressure from agricultural land use, the intensity of land consolidation investment, and the broader socio-economic context that shapes cropland management decisions. All factors were aggregated to the 1 km analysis grid to match the spatial units of the ecosystem service interaction metrics. The source HFP dataset provides observations for 2000, 2010, and 2018. For each grid cell, a linear temporal trend was fitted to these three observations and extrapolated to 2020. The resulting 2020 HFP layer was aligned with the common 1-km analysis grid and used in the 2020 driving-factor analysis.

GeoShapley-based interpretable machine learning model

To capture nonlinear effects of socio-ecological drivers on cropland ecosystem-service patterns and interactions, we combined a machine-learning model with GeoShapley attribution (Covert et al., 2021; Li, 2024; Rao et al., 2019). Separate models were fitted for the 2020 individual ecosystem services and for the recalculated 2020 GWR local Pearson correlation coefficients of WY-SC, WY-FP, and SC-FP on the nominal 1-km grid. The substantive predictors comprised HFP, LACC, PET, POP, PRE, TMP, ELE, SL, and GDP. Residual spatial structure was represented by standardized projected coordinates, polynomial coordinate terms, and a 10 × 7 radial-basis grid, which were treated as a single grouped geographic component (GEO) during attribution. Detailed preprocessing and model settings are provided in Supplementary Methods S5–S6.

For a given sample x and feature j, the Shapley value ϕj is defined in Equation 16: where M is the complete set of features, S is a subset excluding j, and fx​(⋅) is the model prediction given the features in the subset. Shapley values satisfy the local accuracy property, as shown in Equation 17: the sum of all feature attributions equals the deviation from the model's average prediction, with ϕ0 = E [f(X)]. Shapley-based attribution was implemented using a grouped-feature decomposition, in which the coordinate-derived variables were treated collectively as the GEO component. Finally, the GeoShapley contribution of factor j at location i is defined as Equation 18:

Mapping GS j,i across the study area produces spatially explicit surfaces showing where and how strongly each driver is associated with the 2020 ecosystem-service patterns or interaction intensities. These patterns describe spatial heterogeneity and should not be interpreted as temporal changes or causal effects.

Although the InVEST simulations were run at 30 m resolution to utilize high-resolution topography and land cover data, all statistical and machine-learning analyses—including the Spearman correlation tests, GWR, and GeoShapley modeling—were conducted on a 1 km grid aligned with the native resolution of the crop harvest area data and the climatic and socioeconomic predictors. The 30 m InVEST outputs were aggregated to the 1 km grid by computing the mean value of each service within each 1 km cell before these analyses. For each response, complete observations were divided into 80% training data and an untouched 20% holdout test set using target-decile stratification and a fixed random seed. Candidate machine-learning settings were compared by five-fold cross-validation within the training data only. Model performance was evaluated using R 2, RMSE, and MAE on the untouched holdout set, together with random five-fold and spatial-block cross-validation. Detailed model settings and validation results are provided in Supplementary Methods S6 and Supplementary Tables S10-S11. Figure 2 reports the untouched holdout results, whereas spatial-block cross-validation provides the more conservative assessment of transferability to spatially separated areas.

The final machine-learning models showed different levels of predictive performance across the three ecosystem-service pairs (Figure 2). On the untouched 20% holdout test sets, the WY-SC model achieved R2 = 0.621 (RMSE = 0.217), whereas the WY-FP and SC-FP models achieved R2 = 0.994 (RMSE = 0.006) and R2 = 0.990 (RMSE = 0.008), respectively. In five-fold spatial-block cross-validation, the corresponding R2 values were 0.130, 0.926, and 0.915. These results indicate strong predictive reliability for the two FP-related interaction models but only moderate random-holdout performance and limited spatial transferability for WY-SC. Accordingly, GeoShapley attribution for WY-FP and SC-FP was interpreted with greater confidence, whereas the WY-SC attribution was treated as exploratory. The high accuracy of the FP-related models reflects prediction of the spatially smoothed GWR response surfaces and should not be interpreted as evidence of causality.

Data sources

The datasets used to quantify the three cropland ecosystem services and to analyze their driving mechanisms were assembled from multiple publicly available sources. Table 2 summarizes the primary data, their spatial and temporal characteristics, and the purposes for which each was employed. For InVEST simulations, all input rasters were resampled to a common 30 m × 30 m grid, projected to the WGS 1984 UTM Zone 49N coordinate system, and masked to the county boundary. For the subsequent statistical and machine-learning analyses—including the Spearman correlation tests, GWR, and GeoShapley modeling—all variables were aggregated to the same 1 km grid aligned with the native resolution of the crop harvest area data and the climatic and socioeconomic predictors. The aggregation was performed by computing the mean value within each 1 km cell.

Spatiotemporal differentiation of cropland use space

The cropland landscape of Yanchuan County underwent substantial spatial reorganization between 2000 and 2020, reflecting the countervailing forces of slope cropland retirement under the Grain for Green program and cropland creation through gully land consolidation. Figure 3 summarizes the changes in total cropland area and the mean elevation of cropland parcels across the three observation years.

Three salient features emerge from these data. First, total cropland area exhibited a modest but continuous decline over the 2 decades, contracting by 341.64 ha from 2000 to 2010 and by a further 761.13 ha from 2010 to 2020, yielding a cumulative net loss of 1,102.77 ha (approximately 1.1% of the 2000 baseline). Rather than indicating a simple trend of cropland abandonment, this net figure masks a far more dynamic underlying process of simultaneous loss and gain. Slope cropland on upper ridges was progressively retired and converted to forest or grassland in accordance with the Grain for Green mandates, while new terraced fields and check-dam land were incrementally established on the valley floor through gully land consolidation initiatives, partially offsetting the area lost on the slopes. The slight net reduction suggests that the pace of cropland creation in the gully floors has not yet fully compensated for the combined effects of ecological retirement programs and urbanization-related land conversion at the county margins.

The trajectory of mean cropland elevation shows a distinct vertical displacement. Between 2000 and 2010, mean elevation increased slightly, suggesting that early cropland adjustment was not simply a removal of the highest parcels. Between 2010 and 2020, mean elevation dropped sharply by over 24 m, coinciding with the scaling-up of gully land consolidation, which created cropland at much lower elevations than the slope cropland it replaced. This vertical shift reflects a broad spatial reorganization of cropland from rain-fed slopes toward valley bottoms, with ecological forest on ridges and intensive agriculture in gullies.

The spatial redistribution of cropland has direct implications for the ecosystem service functions examined in this study. The downward movement of the cropland centroid is likely to have altered the hydrothermal conditions, soil properties, and management intensity under which crops are grown. Gully-floor cropland generally benefits from deeper soils, higher soil moisture availability, better access to irrigation and mechanization, and greater protection from slope runoff, conditions that favor higher per-unit food productivity and potentially enhanced soil conservation and carbon sequestration compared with extensively managed slope parcels. Yet, these newly consolidated fields may also be subject to heightened flood risk and may entail trade-offs with water yield services at the local scale.

Spatiotemporal patterns of cropland ecosystem services

The three cropland ecosystem services—water yield (WY), food production (FP), and soil conservation (SC)—were quantified at the pixel level using the InVEST model and statistical downscaling methods, and subsequently aggregated to the county scale for each observation year (Figure 4). All three services exhibited net increases over the two-decade period, although the magnitude and trajectory of change differed substantially among them.

Water yield experienced the most pronounced increase, rising nearly tenfold from approximately 2.33 × 10 6 t in 2000 to 2.35 × 10 7 t in 2020, with the strongest acceleration occurring between 2010 and 2020. Soil conservation also showed substantial growth, more than doubling from 9.91 × 10 5 t to 2.16 × 10 6 t over the study period. Food production increased from 18,165 t in 2000 to 45,917 t in 2010, before declining to 36,164 t in 2020, resulting in a net gain of approximately 99% relative to the 2000 baseline but a reduction of about 21% from the 2010 peak. The average FP per unit cropland area followed a similar trajectory, rising from 0.181 t/ha in 2000 to 0.459 t/ha in 2010 and then declining to 0.364 t/ha in 2020.

Spatially, all three services exhibited a broadly increasing trend across the majority of the county, with the most pronounced gains concentrated in the western sector of Yanchuan County. This spatial pattern is closely associated with the heterogeneous distribution of land-use interventions across the landscape. The western portion of the county encompasses extensive tracts of the Qingjian River valley system, where the gully land consolidation project was implemented at larger scales and with higher investment intensity during the 2011–2020 period. The conversion of previously unproductive gully bottoms and floodplains into high-yield check-dam cropland and terraced fields may have contributed to enhanced soil-moisture retention, sediment trapping, and crop productivity in these areas. In contrast, the eastern and northeastern sectors, which border the Yellow River and are characterized by steeper, more dissected terrain with thinner loess cover, experienced comparatively modest service enhancements, likely constrained by both the physical limitations of slope cropland and the lower density of engineering interventions.

The decadal dynamics further reveal the phased nature of ecosystem service responses. The moderate increases observed between 2000 and 2010 primarily reflect the early effects of the Grain for Green program: vegetation recovery on retired slope cropland began to enhance soil retention and water regulation, while cropland contraction initially dampened food production gains. The accelerated growth in WY and SC after 2010 occurred within a period characterized by multiple concurrent processes, including the scaling-up of gully land consolidation, continued ecological recovery of retired cropland, and interannual climatic variability. These processes collectively restructured the spatial configuration of cropland and may have contributed to the amplification of multiple services. The post-2010 decline in total food production, despite the expansion of high-yield valley cropland, may reflect the combined effect of continued retirement of remaining slope cropland and the time lag required for newly consolidated cropland to reach full productive capacity.

Overall, the spatiotemporal patterns documented here indicate that the combined implementation of ecological restoration on ridges and agricultural intensification in gullies has generated a net enhancement of cropland ecosystem services at the county scale, with the western valley systems emerging as the primary locus of multifunctional service gains.

Cropland ecosystem service interactions

Spearman’s rank correlation analysis was applied to identify pairwise trade-off and synergy relationships among the three cropland ecosystem services—water yield (WY), food production (FP), and soil conservation (SC)—across the three observation years (2000, 2010, and 2020). The resulting correlation coefficients are summarized in Figure 5.

The GWR-based local correlation analysis revealed markedly different patterns across the three service pairs. WY–SC exhibited a stable and moderate positive association throughout the study period, with local correlation coefficients of 0.346, 0.373, and 0.377 in 2000, 2010, and 2020, respectively. The bandwidth for this service pair remained constant at approximately 20 neighbors (about 2.13 km) across all three periods, suggesting that the spatial scale over which water yield and soil conservation covary did not change substantially over time, and that their positive association is robust at a relatively fine spatial grain. This persistent positive covariation likely reflects the shared control of precipitation on both runoff generation and erosion intensity, as well as the physically consistent tendency for areas with higher water retention to also exhibit greater sediment trapping.

In contrast, the relationship between SC and FP was weakly negative in all three periods, with local correlation coefficients of −0.098, −0.038, and −0.100 in 2000, 2010, and 2020, respectively. The magnitude of these coefficients remained small (all below |0.11|), and the values fluctuated only modestly over time. This pattern indicates a persistent but weak spatial trade-off between soil conservation and food production at the 1 km scale, rather than the progressive strengthening of synergy reported in the earlier analysis. The WY–FP relationship was even weaker, with local correlation coefficients of −0.013, 0.019, and 0.034 across the 3 years. These near-zero values, combined with the direction that shifted from weakly negative to weakly positive, suggest that no meaningful or stable spatial association existed between water yield and food production over the study period. The bandwidth for FP-related models decreased substantially from 2,200 neighbors in 2000 to 900 in 2010 and 300 in 2020, indicating that the spatial scale of FP-related relationships contracted over time, which may reflect the increasingly localized influence of cropland management and land consolidation on food production patterns.

Taken together, these results indicate that the three service pairs followed distinct trajectories rather than a unified trend toward synergy. Only WY–SC showed a consistently stable and statistically meaningful positive association across all three observation years. The SC–FP pair remained in a weak negative relationship, and the WY–FP pair exhibited effectively no spatial association. The marked differences among service pairs, together with the heterogeneous spatial scales revealed by the changing bandwidths, underscore the necessity of interpreting ES interactions in a pair-specific and scale-sensitive manner. This also implies that the spatial organization of cropland under the coupled policy interventions produced differentiated effects on different service combinations rather than a uniform enhancement of multifunctional synergy.

Analysis of factors influencing the evolution of cropland ecosystem services

The GeoShapley analysis revealed distinct driving factor portfolios for the three cropland ESs (Figure 6), although model performance varied considerably among them. The WY model showed strong predictive ability (holdout R 2 = 0.956; spatial CV R 2 = 0.918), supporting reliable attribution, whereas the SC model performed moderately (holdout R 2 = 0.298; spatial CV R 2 = 0.234), suggesting that the selected predictors explain only part of its spatial variation. The FP model exhibited poor spatial generalization (holdout R 2 = 0.480; spatial CV R 2 = −0.084), indicating that its spatial patterns were dominated by spatial autocorrelation rather than by the explanatory variables; the FP results are therefore interpreted only as weak directional signals.

For water yield, temperature was the dominant natural driver with a predominantly negative effect, consistent with amplified evapotranspiration under the semi-arid climate, while GDP exerted a consistently positive influence, reflecting the role of economic investment and irrigation infrastructure. Slope showed spatially variable effects, promoting runoff generation on steep terrain but favoring retention in flat depositional areas. For soil conservation, slope was the leading positive driver—steeper areas yield greater potential erosion and thus larger absolute reductions—whereas temperature had a negative net effect, likely through accelerated organic matter decomposition and reduced vegetation cover; GDP contributed positively, indicating the benefits of investment-driven restoration and land consolidation. For food production, no factor could be robustly quantified due to the model’s weak spatial performance, though temperature appeared negatively associated and human footprint intensity weakly positive. Across all services, pronounced spatial non-stationarity was evident, confirming that factor effects vary substantially across locations and that single global interpretations are insufficient.

Building upon the individual-service analysis, we examined the drivers of the 2020 geographically weighted local Pearson r values using GeoShapley (Figures 7–9). The WY-FP and SC-FP models showed strong holdout and spatial-block performance (R 2 = 0.994/0.926 and 0.990/0.915, respectively), whereas WY-SC showed moderate holdout performance (R 2 = 0.621) but weak spatial transferability (R 2 = 0.130). The grouped GEO component accounted for 95.5%, 96.2%, and 70.5% of the summed mean absolute GeoShapley magnitude for WY-FP, SC-FP, and WY-SC, respectively. Among the nine substantive predictors, precipitation ranked first for WY-FP, nighttime light intensity ranked first for SC-FP, and GDP ranked first for WY-SC. Their mean signed contributions were negative, positive, and negative, respectively. Because GEO dominated all three models, these measured-driver effects should be interpreted as comparatively small conditional associations rather than primary causal controls; the WY-SC attribution is additionally exploratory because of its weak spatial-block performance.

The decomposition into non-geographic and geographic-interaction components indicated pronounced spatial non-stationarity, although the grouped GEO component was substantially larger than any single measured driver. Among the substantive predictors, precipitation had the largest mean absolute contribution for WY-FP and showed a negative mean signed contribution, whereas nighttime light intensity ranked first for SC-FP and showed a positive mean signed contribution. GDP ranked first for WY-SC, but this result should be treated as exploratory because the model showed weak spatial transferability. The dependence plots further indicate nonlinear and location-dependent responses; therefore, the signs of the mean contributions should not be interpreted as uniform effects across all grid cells.

Evolution of cropland ecosystem services in the Loess Hilly-Gully Region

The spatiotemporal dynamics documented in this study reveal a fundamental reconfiguration of cropland ecosystem services in Yanchuan County under the coupled implementation of the Grain for Green program and gully land consolidation. Contrary to the widely reported trade-offs between provisioning and regulating services in agricultural landscapes (Gong et al., 2024; Yang et al., 2025), our results demonstrate that the three services increased in county-level totals over 2000–2020, but their spatial relationships did not move uniformly toward synergy. Instead, the interaction results showed pair-specific patterns: a stable positive WY-SC association, a weak negative SC-FP association, and a near-zero WY-FP association.

The GWR-based analysis revealed that the three ES pairs followed distinct trajectories rather than a unified trend toward synergy. WY–SC exhibited a stable moderate positive association across the period 2000–2020, whereas SC–FP remained weakly negative and WY–FP showed no stable spatial association (Figure 5). These pair-specific patterns depart from the simplified dichotomy of trade-offs versus synergies often assumed in agricultural landscape studies (Gong et al., 2024; Yang et al., 2025). This departure from the conventional trade-off narrative can be attributed to the distinctive spatial reorganization of cropland from erosion-prone slopes to sediment-trapping valley floors (Li et al., 2013). As slope cropland was retired and replaced by forest and grassland on the ridges, the overall area under cultivation contracted only marginally, while the mean elevation of cropland dropped sharply between 2010 and 2020. This vertical displacement concentrated agricultural activity in the valley bottoms, where deeper soils, higher moisture availability, and greater protection from runoff simultaneously bolster crop yields and soil retention, while the enhanced infiltration and reduced evaporative losses in these low-lying areas amplify water yield (Zhou and Cao, 2020; Huang et al., 2022).

The underlying mechanism linking this spatial reorganization to changes in service interactions can be understood as follows: the spatial redistribution of cropland from slopes to valley floors altered local hydrological and soil processes, which in turn reshaped the relationships among services. On sloping cropland, steep gradients and thin soils promote rapid runoff generation, which entrains sediment and nutrients while leaving insufficient water for crop uptake, thereby sustaining the classic trade-off between food production and soil conservation. When these slopes are retired and cropland is re-established on valley floors through gully land consolidation, the newly created check-dam land and terraces intercept runoff and sediment from upstream slopes. This results in deeper, moisture-retentive soils that facilitate infiltration, reduce erosive surface flow, and enhance plant-available water. Consequently, the same plots that achieve higher grain yields also exhibit stronger soil conservation and greater water retention, shifting the local relationship from trade-off toward synergy. Observations from representative consolidation projects in Yanchuan County are consistent with this interpretation: after gully land consolidation, grain yields on newly created cropland increased substantially, and the improved land conditions concurrently enhanced soil and water retention functions (Qu et al., 2022). While such field evidence does not independently establish causality, the temporal alignment and spatial correspondence between consolidation intensity and the strengthening of service synergies provide reasonable support for the proposed mechanism.

The temporal phasing of these changes is particularly instructive. The moderate gains in water yield and soil conservation between 2000 and 2010 largely reflect the early eco-hydrological benefits of vegetation recovery on retired cropland (Lu et al., 2018; Ji et al., 2024), while the acceleration after 2010 coincides with the scaling-up of gully land consolidation (Qu et al., 2022). The temporal stability of the WY–SC association suggests that the coupling between water yield and soil conservation is governed by relatively persistent biophysical processes—principally the shared control of precipitation on runoff generation and erosion intensity—that are less sensitive to land-use reorganization. By contrast, the persistent weak negative SC–FP association indicates that the spatial trade-off between soil conservation and food production, although modest, has not been eliminated by the vertical reorganization of cropland. The near-zero WY–FP association further suggests that water availability and food production are spatially decoupled at the 1 km scale, likely because crop yields are controlled by a combination of soil quality, management intensity, and microclimatic conditions rather than by water yield alone.

The spatial heterogeneity revealed by the geographically weighted analysis and GeoShapley indicates that the observed associations, particularly the positive WY-SC association, were not spatially uniform The western valley systems, where gully land consolidation was most intensively implemented, showed the strongest WY–SC associations, whereas the eastern uplands exhibited weaker or more variable patterns. This spatial contrast suggests that the effects of land reorganization on ES interactions are context-dependent and cannot be generalized uniformly across the landscape. Recognizing this internal heterogeneity is essential for designing spatially targeted interventions rather than uniform management prescriptions across the entire landscape.

Driving mechanisms of cropland ecosystem service dynamics

The GeoShapley attribution results illuminate the contrasting roles of natural and anthropogenic factors in shaping individual ES patterns, although the strength of evidence varies considerably across services because of unequal model reliability. For water yield, the model performed well in both random holdout and spatial cross-validation (holdout R 2 = 0.956; spatial CV R 2 = 0.918), allowing reasonably robust attribution. Temperature was the dominant natural driver with a predominantly negative effect, consistent with the well-established role of heat stress in amplifying evapotranspiration and reducing water availability in semi-arid environments (Hu et al., 2020; Ma et al., 2022). GDP exerted a consistently positive influence, suggesting that economic investment in irrigation and land improvement has partially mitigated climatic constraints on water yield. For soil conservation, the model showed only moderate performance (holdout R 2 = 0.298; spatial CV R 2 = 0.234), so its attribution results should be interpreted with corresponding caution. Nevertheless, slope emerged as the leading positive driver, which is physically intuitive given that steeper terrain generates greater potential erosion and thus greater absolute erosion reduction (Guo et al., 2021), whereas temperature exhibited a negative net effect, likely through accelerated organic matter decomposition and reduced vegetation cover. For food production, the model failed to generalize spatially (holdout R 2 = 0.480), indicating that its apparent spatial patterns were dominated by spatial autocorrelation rather than by the explanatory variables. Therefore, robust quantitative attribution is not possible for FP, and only weak directional indications—such as a negative association with temperature and a weakly positive association with human footprint intensity—can be reported.

The interaction models yielded strong holdout and spatial-block performance for WY-FP (R 2 = 0.994 and 0.926) and SC-FP (R 2 = 0.990 and 0.915), respectively. Precipitation and nighttime light intensity were the highest-ranked substantive predictors for WY-FP and SC-FP, although their contributions were small relative to GEO. WY-SC achieved a holdout R 2 of 0.621 but a spatial-block R 2 of only 0.130. Its leading substantive predictor, GDP, and its large geographic component should therefore be interpreted as exploratory rather than as evidence of a spatially transferable mechanism. More generally, the dominance of GEO indicates that the models reproduced the smooth spatial structure of the GWR-derived response surfaces more strongly than they isolated causal effects of the measured drivers.

Extending the analysis to pairwise ES interaction intensities produced further differentiated evidence. The interaction models yielded strong holdout and spatial-block performance for WY-FP (R 2 = 0.994 and 0.926) and SC-FP (R 2 = 0.990 and 0.915), respectively. Precipitation and nighttime light intensity were the highest-ranked substantive predictors for WY-FP and SC-FP, although their contributions were small relative to GEO. WY-SC achieved a holdout R 2 of 0.621 but a spatial-block R 2 of only 0.130. Its leading substantive predictor, GDP, and its large geographic component should therefore be interpreted as exploratory rather than as evidence of a spatially transferable mechanism. More generally, the dominance of GEO indicates that the models reproduced the smooth spatial structure of the GWR-derived response surfaces more strongly than they isolated causal effects of the measured drivers (Jing et al., 2024).

It should be noted that ES interactions are likely scale-dependent (Xia et al., 2023; Yang et al., 2025); because most source datasets have a native resolution of approximately 1 km, all analyses presented here were conducted at the 1 km scale, and the inferred interaction patterns may not hold at finer or coarser spatial grains. The pronounced spatial non-stationarity observed across both individual-service and interaction-intensity models remains an important finding, but it must be interpreted in light of the model validation results. For well-performing models (WY, SC–FP, WY–FP), geographic variation in factor effects is more likely to reflect genuine spatial heterogeneity in underlying ecological and socioeconomic processes. For poorly generalizing models (FP, WY–SC), part of the apparent spatial structure may instead reflect residual spatial autocorrelation that the predictors did not adequately capture. This distinction cautions against blanket management prescriptions and underscores the need for spatially targeted interventions informed by more robust models and, ideally, field-based validation. In sum, the interaction intensity between cropland ESs is not a simple linear extension of individual service drivers but emerges from spatially heterogeneous interactions between natural gradients and anthropogenic interventions. However, the strength of the evidence varies across service pairs, and further work with improved predictors and validation data is needed before these patterns can support strong policy conclusions.

Policy implications and research prospects

The findings of this study carry several implications for the sustainable management of cropland ecosystems in the Loess Hilly-Gully Region. First, the stable positive WY–SC association, together with the differentiated patterns observed for the other service pairs, provides partial empirical support for the “Grain for Green on the ridges, gully land consolidation on the valley floor” model as a spatial strategy. However, the findings also indicate that this strategy does not uniformly generate synergies across all service combinations. Rather than framing ecological restoration and food security as competing objectives, policy should continue to promote the spatial separation and functional integration of these goals: ecological forests on steep slopes deliver regulating services, while intensified farming on consolidated valley land secures food supply, and the two zones together generate landscape-level complementarities. However, the moderate magnitude of the observed spatial associations and their localized variability caution against interpreting the county-wide pattern as evidence that all consolidated sites have achieved functional synergy. The spatial heterogeneity revealed by GWR and GeoShapley analyses underscores the need for spatially targeted measures rather than uniform prescriptions.

Second, the positive association between GDP and water yield and, to a lesser extent, soil conservation suggests that sustained investment in land engineering and irrigation infrastructure may help maintain these services. However, the effect of economic investment on food production remains uncertain because of the limited spatial generalizability of the FP model. However, the extension of GeoShapley analysis to ES interaction intensities reveals a more nuanced picture: anthropogenic factors such as GDP and population density exerted mixed and sometimes negative effects on pairwise interaction strengths, with the direction and magnitude varying considerably across locations and service pairs. This indicates that economic investment alone does not guarantee enhanced synergies; its effects on ES interactions depend on local topographic and hydrological contexts. Priority should therefore be given to areas where natural constraints are severe but engineering interventions have proven capable of shifting the balance toward stronger positive associations, as observed in the western valley systems of Yanchuan County, while in eastern uplands with steep terrain and thinner soils, complementary conservation measures may be needed to avoid exacerbating local trade-offs.

Several limitations of this study point to avenues for future research. The analysis was conducted at a single county, and while Yanchuan is broadly representative of the Loess Hilly-Gully Region, the transferability of its findings to areas with different geological substrates, precipitation regimes, or land-tenure systems remains to be tested. The exclusion of carbon sequestration—a key ecosystem service in the Loess Plateau—from the service portfolio is a notable gap; incorporating carbon storage into the interaction and attribution framework would provide a more complete picture of cropland multifunctionality and its climate-mitigation implications. The GeoShapley analysis, while powerful in revealing nonlinear and spatially heterogeneous relationships, is correlational in nature and does not establish causality; coupling machine-learning attribution with process-based models or quasi-experimental designs could strengthen causal inference. Importantly, the interpretation of spatial non-stationarity in this study currently relies on model-based analysis and indirect inference from published engineering cases in the region. Direct validation using plot-level, long-term field monitoring of soil moisture, sediment retention, and crop yield within Yanchuan County is still lacking. Future research integrating such field measurements with the spatial analytical framework adopted here would be essential to verify the hypothesized “spatial displacement–process alteration–service reorganization” mechanism. Finally, the temporal resolution of three decadal snapshots cannot fully resolve year-to-year variability driven by climatic fluctuations or capture lagged responses of ecosystem services to land-use change. Dense time-series analyses coupled with long-term field observations are needed to disentangle the effects of policy interventions from those of climatic variability and to identify early-warning signals of potential shifts in service relationships.

In summary, this study demonstrates that the spatial reorganization of cropland under coordinated ecological restoration and gully land consolidation is associated with a shift from historically conflict-ridden agriculture–environment relationships toward strengthened spatial complementarities among multiple ESs on the Loess Plateau. The results underscore the importance of considering spatial configuration—not merely area—in land-use policy, and highlight the value of interpretable machine learning in unpacking the spatially heterogeneous mechanisms that underpin ES interactions. These insights provide a scientific basis for designing more resilient and multifunctional agricultural landscapes in semi-arid regions undergoing rapid land-use transitions.

This study investigated the spatiotemporal patterns and driving mechanisms of interactions among cropland ecosystem services in Yanchuan County, a typical Loess Hilly-Gully Region, under the coupled influences of the Grain for Green Program and gully land consolidation. The main conclusions are as follows. All three cropland ecosystem services achieved net gains at the county scale over 2000–2020, although the magnitude and trajectory of change differed among them. Water yield increased nearly tenfold, soil conservation more than doubled, and food production rose from 18,165 t in 2000 to 45,917 t in 2010 before declining to 36,164 t in 2020, yielding a net gain of approximately 99% from the 2000 baseline but a reduction from the 2010 peak. The most pronounced gains were concentrated in the western valley systems, where deeper soils, better moisture conditions, and relatively higher investment in cropland improvement measures coincided with stronger gains in water yield and soil conservation. The post-2010 decline in food production likely reflects the combined effects of continued slope cropland retirement, the time required for newly consolidated cropland to reach full productive capacity, and fluctuations in agricultural management intensity.

The three cropland ES pairs exhibited markedly different interaction patterns rather than a uniform trend toward synergy. Based on the revised GWR-based local analysis, WY–SC showed a stable moderate positive association across all three periods, indicating a persistent spatial coupling between water yield and soil conservation. SC–FP exhibited a weak negative association throughout the study period, suggesting a limited but persistent trade-off tendency between soil conservation and food production at the pixel level. WY–FP displayed effectively no stable spatial association, with coefficients fluctuating around zero. These results indicate that the interactions are pair-specific and spatially heterogeneous, with stronger positive WY–SC associations in the western valley systems and weaker or more variable patterns in the eastern uplands.

The driving mechanisms of individual ESs and their pairwise interactions were both spatially non-stationary, but the strength of evidence varied across models. For water yield, temperature was the dominant negative natural driver, while GDP played a consistently positive role. For soil conservation, slope was the leading positive driver and temperature had a negative net effect, although the model’s moderate spatial generalization calls for caution. For food production, the evidence is more limited: the FP model exhibited moderate random holdout performance but weak spatial generalization, and its attribution results should therefore be interpreted as preliminary directional indications rather than definitive quantitative estimates. For pairwise interactions, the WY-FP and SC-FP models showed strong holdout and spatial-block performance, with precipitation and nighttime light intensity ranking as their leading substantive predictors, respectively. The grouped geographic component nevertheless dominated both models. The WY-SC model showed limited spatial transferability; consequently, its large geographic component and the attributed effect of GDP should be interpreted as exploratory. These findings underscore that the driving mechanisms of ES interactions are context-dependent, and that both model reliability and spatial heterogeneity must be considered in mechanistic interpretation.

Overall, the spatial reorganization of cropland under the coupled implementation of the Grain for Green program and gully land consolidation is associated with differentiated changes in ES interactions: stable positive associations for some service pairs, weak negative or negligible associations for others, and pronounced spatial heterogeneity throughout. These results highlight the importance of considering spatial configuration and pair-specific relationships in land-use policy, and demonstrate the value of interpretable machine learning for uncovering spatially heterogeneous mechanisms, while also underscoring the need for rigorous model validation and field-based evidence to support policy-relevant conclusions.