Determining the driving mechanism of ecosystem services provided by cropland in China

Abstract

Cropland provides the base for the rising food demand. Recent studies highlighted the more systematic cropland management from ecosystem services (ESs) perspective, rather than solely depending on the provisioning services. However, the driving mechanism of cropland related ESs has not been well understood, hindering the co-benefits of sustainable crop-management practices. Taking the world’s largest grain producer, China, as a study case, this work contributed to determine the dynamic changes and driving mechanism of cropland related ESs by integrating InVEST model, geographically weighted regression and partial least squares-structural equation modeling (PLS-SEM). It was found that, from 2000 to 2020, the increased water purification (WP) was along with the decreased grain production (GP), soil retention (SR), carbon storage (CS) and water retention (WR) which mainly occurred in southern China. Aggregation of ESs were identified with strong spatial heterogeneity, while CS-WR and WR-WP almost showed persistent trade-offs. The direct effects of climate on WR and WP, vegetation on GP, terrain on SR and human activities on CS were found to be most significant. While the opposite effects of climate on GP during the study period were also observed. The positive effect of human activities on GP and WR was lessen by the indirect effect of terrain, highlighting the synergy management of interaction mechanisms among drivers. Finally, implications of these findings were discussed.

Graphical abstract

Keywords

Cropland ecosystem services / driving mechanism / dynamic changes / synergy management / trade-offs

Highlight

	● Southern China is the vital area for cropland ESs but with substantial decline.
	● CS-WR and WR-WP exhibit persistent trade-offs relationship.
	● WR and WP of cropland were direct impacted by climate most.
	● The positive effect of human activities on GP and WR was lessen by the indirect effect of terrain.

1 Introduction

Food, feed and fiber, classified into provisioning ecosystem services (ESs) from cropland, are essential for human survival. Generally, the cropping practices are determined by these provisioning ESs^[1]. Recent studies have recognized the more important role of cropland ESs, such as water retention (WR), carbon storage (CS), water purification (WP), and soil retention (SR), in supporting the UN Sustainable Development Goals (SDGs) related to food security, climate mitigation, and water conservation^[2]. China faces substantial challenges to realize these ESs given that 20% of the global population is fed by only 8% of the global arable land and 6% of water resources, indicating the significant anthropogenic pressures on agroecosystem functions^[3]. To respond the land-system sustainability emergency, China enacted large-scale programs related to cropland since 2000, such as the Grain for Green Program, the Natural Forest Conservation Program, and the Comprehensive Agricultural Development Program^[4]. However, the dynamic changes and the related driving mechanism of cropland ESs still remain unknown, hindering optimization strategies to navigate the sustainability trade-offs.

Intensive studies have attempted to bring ESs into decision making^[5–^11]. Specified to cropland ESs, various evaluation frameworks have been well proposed which can be generally grouped into physical flows based approaches^[12], ecological economics based accounting^[13] and monetary based method^[14]. Understanding the relationship among various ESs can prepare systematic implications to avoid unintended consequences. For example, Hu et al.^[15] demonstrated that total nitrogen export and water yield can be greatly improved by paddy land-to-dry land (PLDL) program in the Erhai River Basin. In the Dongting Lake Basin, Zheng et al.^[16] indicated that agricultural land use intensity had greater impacts on the trade-offs between grain production, water purification than land use types. Tian et al.^[17] demonstrated that the increased cultivated land can increase water yield greatly, but decrease the soil retention in the Zhangcheng District. Ma et al.^[18] found, in the Jianghuai Ecological Economic Zone, that crop production was a trade-off with carbon storage, habitat quality and water purification under the returning farmland to lakes policy.

Factors impacting cropland ESs, such as planting mode^[19], policy intervention^[20], landscape composition and configuration^[21], land use change^[22], crop diversity^[23] and climate changes^[24], have also been examined. These studies generally used correlation analysis and regression method, however, may ignore the combined effects of these factors due largely to the impacts of these factors are intertwined. The driving factors can have mutual influences on each other. For example, the terrain enhances cropland ES value indirectly through its impact on vegetation, human activities, and soil characteristics^[25]. To this content, systems approach integrating both socioeconomic and environmental drivers should be adopted to address the coupled human and natural systems^[26]. Recently, partial least squares-structural equation modeling (PLS-SEM) has been used to examine the driving mechanism of ESs, attributed to its ability to detect both the direct and indirect effects of drivers, as well as their complex interactions^[27,^28]. Their results show that beyond the direct drivers, such as climate, urbanization and geomorphic features, the interaction of drivers (such as economic development and meteorological) also shows the indirectly considerable impacts on ESs. Determining these mechanisms therefore provides more systematic picture toward sustainable management of ESs.

The above studies provide insightful implications for cropland management based on ESs. However, they are generally restricted to relatively small scale and short-term data, lacking profound insights into long-term trends at the macro level. Especially in China where the terrain and climatic conditions are of great difference, a comprehensive understanding of the changes in cropland ESs is particularly urgent. In addition, most studies focus on a few driving factors, by using established methods including correlation analysis and geographic detector modeling, while failing to examine the interactions among multiple factors. The insightful implications cannot be raised without the in-depth understanding of the interactions among multiple factors, especially given the significant heterogeneities in the cropland system of China. Under such a circumstance, this paper aims to answer two questions: what is the dynamic changes of cropland ESs in China, and what is the driving mechanism of cropland ESs from both ecological and socioeconomic perspective. The study framework (Fig.1) was used to answer these questions. We first evaluate the ESs provided by cropland in China including grain production, water retention, carbon storage, soil retention and water purification in 5-year epochs from 2000 to 2020. The aggregation of cropland ESs and their related relationships are examined by adopting G_i^* statistics method and geographically weighted regression (GWR), respectively. Finally, PLS-SEM is used to reveal the driving mechanism of the cropland ESs.

2 Method

2.1 Introduction of the cropland in China

The cropland in China was preliminarily classified as paddy or dry land^[29]. This was divided into nine agricultural zones that have been officially identified in China as listed in Tab.1. Due to limited data sources in Hong Kong, Macao and Taiwan, the study excluded these regions. There are significant heterogeneities in the cropland system of China in terms of ecological condition. For example, the highest annual precipitation and average temperature are found in SC, reaching to 1487–1909 mm and 21.8–25.1 °C, respectively. While the annual precipitation in the NASR is only between 146 and 391 mm. Also, the lowest temperature is observed in NCP ranging from 5.5 to 9.2 °C.

2.2 Data source and processing

Seven categories of data were used in this study are as follows. (1) Land use/land cover (LULC) data. The LULC data covering 2000, 2005, 2010, 2015, and 2020 were from Resources and Environment Science and Data Center^[29]. The LULC descriptions were segmented using a two-level classification system in this study, based on their natural properties of land resources. (2) Meteorological data. Annual average precipitation was also derived from Resources and Environment Science and Data Center^[30], Evapotranspiration was from the Institute of Tibetan Plateau Research Chinese Academy of Sciences^[31]. Temperature was from the National Tibetan Plateau Data Center^[32]. (3) Topographic data. Digital elevation model (DEM) was obtained from Geospatial Data Cloud site. Elevation and slope were based on DEM data by ArcGIS. (4) Soil data. Data on the content of sand, silt, clay, organic, pH and bulk were taken from the Harmonized World Soil Database version 1.2^[33]. (5) Vegetation data. The normalized difference vegetation index (NDVI) came from the MOD13A3 remote sensing data within the MODIS data set^[34]; enhanced vegetation index (EVI) and leaf area index (LAI) came from Earth Data^[35,^36]. (6) Socioeconomic data. Gross domestic product (GDP) and population density (POP) came from Resources and Environment Science and Data Center^[37,^38]. Human footprint (HFP) and night light index (NTL) were taken from Mu et al.^[39] and Zhong et al.^[40], respectively. (7) The division of nine agricultural zones was from Resources and Environment Science and Data Center. The detailed data source and relevant information are shown in Tables S1–S4. All data were processed by ArcGIS software, with a uniform resolution of 1 km × 1 km.

2.3 Evaluating the cropland ecosystem services

This study used the ArcGIS platform for spatialization of grain yield statistics. The equivalence factor method, emergy analysis and InVEST model are commonly adopted to evaluate the physical flows of ESs. Of these, the InVEST model is a powerful tool with the advantages of spatial visualization effect and relatively easy access to data^[41]. Therefore, the InVEST model was used to assess cropland ESs.

2.3.1 Grain production

GP is essential for both food security and regional sustainability^[42]. A significant linear correlation is revealed between NDVI and GP^[43]. Based on GP statistics, GP was spatially assigned to cropland using wholly positive NDVI values in agricultural district units, with production data obtained from National Bureau of Statistics, as:

(1)G_i=G_{sum}\times \frac{{NDVI}_i}{{NDVI}_{sum}}

where, G_i (t) is the grain production of the ith cropland grid, G_sum (t) is the total grain production of the study area, NDVI_i is the annual NDVI value of the ith cropland grid, and NDVI_sum is the sum of NDVI value for cropland.

2.3.2 Water retention

WR indicates the water reserved in ecosystems over a given period of time^[44]. WR was expressed as evapotranspiration and runoff subtracted from precipitation as:

(2){WR}_i=Y_i-{Runoff}_i

(3)Y_i=\left(1-\frac{{AET}_i}{P_i}\right)\times P_i

(4){Runoff}_i=P_i\times C_i

where, for each pixel i, WR_i (mm·yr⁻¹) is water retention, Y_i (mm·yr⁻¹) is the annual output of water, P_i (mm·yr⁻¹) is the annual precipitation, AET_i (mm·yr⁻¹) is the actual annual evapotranspiration, Runoff_i (mm·yr⁻¹) is the annual average runoff and C_i is the runoff coefficient.

2.3.3 Soil retention

SR was estimated using the revised universal soil loss equation as:

(5)D_i=R_i\times K_i\times {LS}_i\times C_i\times P_i

(6){SR}_i=R_i\times K_i\times {LS}_i\times (1-C_i\times P_i)\times {SDR}_i

(7)SD{R}_i=\frac{SD{R}_{max}}{1+\exp \left({\displaystyle \frac{I{C}_0-I{C}_i}{k}}\right)}

where, D_i is the annual soil loss, SR_i is the soil retention, R_i (MJ·mm·(ha·h)⁻¹) is the rainfall erosivity, K_i (t·ha·h·(MJ·ha·mm)⁻¹) is the soil erodibility, LS_i is the slope length-gradient factor, C_i is the crop-management factor, and P_i is the support practice factor. SDR_i is determined by the connectivity index IC, SDR_max is the highest theoretical SDR, while the calibration parameters IC₀ and k determine the shape of the relationship between SDR and IC.

2.3.4 Carbon storage

The carbon storage and sequestration model in the InVEST model was used to evaluate CS as:

(8)C_{total}=C_{above}+C_{below}+C_{soil}+C_{dead}

where, C_total is the total carbon storage, C_above, C_below, C_soil, and C_dead are the aboveground carbon stock, underground carbon stock, soil carbon storage, and dead organic matter carbon stock, respectively.

2.3.5 Water purification

The Nutrient Delivery Ratio model can evaluate the nutrient retention ability of an ecosystem, providing a proxy for WP function. The nitrogen export value is selected as the indicator for evaluating WP. Higher nitrogen export represents lower WP capacity^[45]. The calculations were:

(9)X_{{export}_{tot}}={\sum }_i{X}_{{export}_i}

(10)X_{{export}_i}={load}_{surf,i}\times {NDR}_{surf,i}+{load}_{subs,i}\times {NDR}_{subs,i}

where, X_{{export}_i} and X_{{export}_{tot}} are the nutrient export at the pixel and watershed scales, respectively, {load}_{surf,i} and {load}_{subs,i} are nutrients transported by surface and shallow subsurface runoff, respectively, NDR is nutrient delivery ratio for a pixel i, {NDR}_{surf,i} is nutrients transported by surface flow, and {NDR}_{subs,i} is the subsurface flow of nutrients.

The results of ESs were standardized from 0 to 1. Considering computational costs and validity of results^[46], we built a 20 km × 20 km grid, and the mean values of ESs at the scale were computed using the Zonal Statistics toolbox in ArcGIS.

2.4 Determining the spatiotemporal aggregation of cropland ES

Spatial analysis was used to detect regions with both high and low value for ESs. We used G_i^{\ast } statistics method presented by Getis and Ord^[47] to perform spatial analysis as:

(11)G_i^{\ast }=\frac{\sum _j^n{W}_{ij}{x}_j}{\sum _j^n{x}_j}

(12)Z(G_i^{\ast })=\frac{G_i^{\ast }-E(G_i^{\ast }\text{)}}{\sqrt{Var(G_i^{\ast })}}

where, G_i^{\ast } is the spatial aggregation index for grid i, W_{ij} is the weight matrix of grids i and j, n is the total amount of grids, x_j is the attribute value for grid i,Z(G_i^{\ast }) is defined as the degree of significance of the aggregation index, E(G_i^{\ast }) is the mathematical expectations, and Var(G_i^{\ast }) is the variances coefficient of G_i^{\text{*}}.

2.5 Identifying the relationships among cropland ecosystem services

Common methods used in the study of ES trade-offs and synergies include multivariate regression tree^[48] and correlation analysis^[49]. While GWR evolved from the standard regression framework and can test spatial non-stationarity^[50], due to its ability to enable the spatial evaluation. Therefore, GWR was used to identify spatial interactive correlations between pairs of ESs, including CS-WP, CS-SR, CS-WR, CS-GP, SR-WP, SR-GP, WR-WP, WR-SR, WR-GP, and WP-GP. Positive and negative regression coefficients between ES pairs indicate synergy and trade-off, respectively. The package “GWmodel” was run in R 4.2.1 to perform GWR^[51] as:

(13)y_i={\beta }_0({\mu }_i,v_i)+{\sum }_{k=1}^p{\beta }_k{x}_{jk}({\mu }_i,v_i)x_{jk}+{\epsilon }_i

where, ({\mu }_i,v_i) is the spatial location of point i, p is the number of argument, y_i is the number of dependent variable, x_{jk} is independent variables, {\epsilon }_i is random errors; {\beta }_0({\mu }_i,v_i)is the intercept at i, and {\beta }_k({\mu }_i,v_i) is the regression coefficient.

2.6 Determination of the driving mechanism of cropland ecosystem services

Structural equation models (SEM) can be typically divided into two types, namely covariance-based SEM (CB-SEM) and PLS-SEM^[52]. Compared to CB-SEM, PLS-SEM is more suitable for exploratory research, due to that its low requirement for sample distribution, and it being less susceptible to multicollinearity bias and can be possessed efficiently based on the related software^[53,^54]. This study examined the relationships between variables by PLS-SEM based on SmartPLS 3.

To build the PLS-SEM of cropland ESs, we selected multiple driving factors according previous studies^[25,⁵⁵–^57] in this study. As shown in Tab.2, PLS-SEM includes both latent and observed variables. The drivers identified and incorporated into the model were changes in climate (ET, PRE and Tempt), soil (BULK, PH and SAND), vegetation (EVI, LAI and FVC), terrain (DEM and Slope) and human activity (GDP, HFP, NTL, POP and Paddy). All variables were first normalized. Meanwhile, the multicollinearity of variables was detected by variance inflation factor with this factor always below 10 (Table S5), which met the criteria^[58].

3 Results

3.1 Dynamic changes of the ecosystem services provided by cropland in China

Fig.2 shows the change in GP. This had an increasing trend from 462 Mt in 2000 to 669 Mt in 2020, with the annual growth rate at 1.87%. In terms of GP per unit area, an increase of 48.0% was found in 2020 from 258 t·km⁻² in 2000. All districts had upward trends in GP except SC (−8.99 Mt), while NCP and NASR had the highest annual growth rate of 4.83% and 4.21%, respectively. MLYP, HHHP and NCP together contributed over 67% to the national GP during the study period.

Fig.2 shows that the total WR increased from 291 × 10⁹ m³ in 2000 to 357 × 10⁹ m³ in 2020, in which the most notable increases were observed in MLYP (28.8 × 10⁹ m³), SBSR (18.0 × 10⁹ m³), and NCP (17.8 × 10⁹ m³). However, WR decreased in SC (−14.7 × 10⁹ m³). During the study period, MLYP provided the largest contribution to total WR, with the shares and amounts ranging from 47.2% and 137 × 10⁹ m³ to 46.6% and 166 × 10⁹ m³, respectively. This was followed by YGP and SC. QTP and LP both had increased in WR, but together this only represented less than 3% of total WR. In terms of unit area, the WR decreased from 163 × 10³ m³·km⁻² in 2000 to 153 × 10³ m³·km⁻² in 2005 and increased to 204 × 10³ m³·km⁻² in 2020.

During the study period, the capability of WP declined considerably because nitrogen export increased by 7.22% (Fig.2). The nitrogen export initially decreased from 332 Gg in 2000 to 314 Gg in 2010, and then increased to 362 Gg in 2015, and finally decreased to 356 Gg in 2020. The largest contribution to total nitrogen export came from MLYP, with the share ranging from 31.4% to 30.3%. From 2000 to 2020, the largest annual growth rate was observed in NCP (1.77%), followed by QTP (1.60%) and NASR (1.06%). While SC decreased most (−1.01%), followed by YGP (−0.26%). From unit area perspective, the upwards trend was found for nitrogen export, with values ranging from 186 kg·km⁻² in 2000 to 203 kg·km⁻² in 2020.

Fig.2 shows that SR had a fluctuating trend. It first decreased from 33.2 Gt in 2000 to 30.7 Gt in 2005, then rose to 33.0 Gt in 2010, and finally dropped to 29.5 Gt in 2020. Between 2000 and 2020, the largest and lowest annual growth rate in SR were observed in QTP (4.25%) and SC (−1.55%), respectively. During the study period, YGP and SBSR together contributed over 58.5% of total SR. The SR unit area decreased from 18.5 kt·km⁻² in 2000 to 16.8 kt·km⁻² in 2020, with the annual growth rate at −0.48%. The main growth occurred in QTP (31.2 kt·km⁻²), while SR decreased the most in YGP (−13.3 kt·km⁻²) and SC (−10.4 kt·km⁻²).

Finally, the changes in CS are shown in Fig.2. Total CS had a downward trend from 2.46 Gt in 2000 to 2.41 Gt in 2020, representing an annual decline of −0.10%. The zones with the highest annual growth rate were located in QTP (1.15%) and NASR (0.51%). In NASR, CS increased notably with 38.62 Mt, while that in MLYP (−35.09 Mt) and HHHP (−24.91 Mt) had a decreasing trend during the study period. The unit area of CS remained relatively steady, but increased slightly from 1.37 kt·km⁻² in 2000 to 1.374 kt·km⁻² in 2020.

3.2 Aggregation of ecosystem services provided by cropland in China

Fig.3 shows dynamic change of high- and low-intensity areas of cropland ESs at different confidence levels from 2000 to 2020, respectively. The spatial distribution of high- and low-intensity areas of GP changed greatly between 2000 and 2020. The main high-intensity areas of GP function were distributed in HHHP and MLYP, with a gradual expansion toward NCP. Meanwhile, NASR, LP and QTP remained low-intensity areas and expanding to SBSR, SC and YGP. Between 2000 and 2020, the high-intensity areas in NCP increased from 0.00% to 98.0% and low-intensity areas in YGP increased from 0.10% to 93.9%. The distribution pattern of WR can be characterized by high-intensity south and low-intensity north pattern. The high-intensity areas in WR were mainly located in SC, MLYP and eastern YGP and SBSR. The high- and low-intensity areas in YGP and SBSR were relatively variable from 2000 to 2020. The spatial distribution of WP and CS were similar, demonstrating a high-intensity north and low-intensity south pattern. For WP, the proportions of high- and low-intensity areas in HHHP and QTP changed significantly from 2000 to 2020. For CS, NCP had the greatest change in the proportion of high-intensity areas, which increased from 71.7% in 2000 to 100% in 2005, and finally decreased to 80.1% in 2020. The high-intensity areas in NASR and low-intensity areas in SC remained flat, both closed to 100%. The spatial pattern of SR fluctuated over 20 years in SBSR, MLYP, YGP, SC, and QTP, with the remaining regions being stable. The high-intensity areas of MLYP, SC, and YGP changed in almost same trend. The proportion of high-intensity areas in NCP and HHHP was unchanged from 2000 to 2020, remaining at 0%.

3.3 Trade-offs and synergies among ecosystem services

The results from GWR analysis indicated that the trade-offs and synergies of cropland ESs in China show heterogeneity (Fig.4). During the study period, we found that CS-WR and WR-WP mainly had a trade-off relationship, while CS-SR had a synergy relationship. In addition, SR-GP and WR-GP generally had a synergy relationship, but had a trade-off relationship in HHHP, and the latter also in MLYP and HHHP. Between 2000 and 2020, CS-GP mainly had spatial trade-offs, while regions in HHHP, MLYP, and SBSR consistently showed synergies in 2010, 2015 and 2020. The trade-offs and synergies for WR-SR were more variable. The ratio of spatial trade-offs was higher than that of spatial synergies for WP-GP and SR-WP in some years, indicating that they mainly had spatial trade-offs. While their spatial trade-off distributions all included the regions NASR, QTP, and YGP. It should be noted that CS-WP had more spatial synergies than spatial trade-offs, mainly in HHHP, MLYP, LP, NASR, NCP, and SC.

3.4 Driving mechanism of cropland ecosystem services

Factor loadings were between 0.474 and 1.00, and significant (p < 0.001), indicating better reflection of latent variables (Tables S6–S7). Most of values for composite reliability were above 0.7 with an average variance extracted above 0.5, which showed that the model reliability and validity were well established^[59,^60] (Tables S8–S12). Predictive relevance of latent variable was greater than 0, indicating predictive validity of model^[61]. The goodness-of-fit in this study was greater than 0.36, revealing the overall performance of the entire model^[62]. Overall, the above indicators demonstrated that the PLS-SEM model built in this study was acceptable.

3.4.1 The direct driving factors of cropland ecosystem services

Fig.5 shows the results of PLS-SEM analysis in 2020, and the results in other studied years are provided in Table S13–S16. It can be found that from 2000 to 2020, vegetation, human activities and soil all positively influenced GP, in which vegetation had the greatest effect. Also, the impact intensity of vegetation started at 0.535 in 2000, increased to 0.671 in 2010, and then decreased to 0.565 in 2020, having a trend of initial increase followed by decrease. The climate effect on GP was positive in 2000 and 2005, but negative since 2010. The effect of terrain on GP was not significant in 2000 and negative from 2005 to 2020.

Climate and human activities had consistent positive effects on WR, with climate having the greatest effect evidenced by path coefficients over 0.6 during the study period. This is not surprising due to WR is determined by evapotranspiration and precipitation which are largely impacted by climate. Soil and vegetation had a slight negative influence on WR. The impact of vegetation changed from −0.015 in 2000 to −0.083 in 2015, and finally returned to −0.015 in 2020. However, the effect of soil showed greater fluctuation over time. It increased from −0.174 to −0.145 in 2000–2010, declined to −0.198 in 2015, then rose slightly to −0.191 in 2020. Terrain had only a slight positive influence on WR in 2000 but it was negative in 2005–2020. Climate, terrain, human activities and vegetation influenced WP negatively, with climate having the greatest effect. Soil was the only driver positively impacting WP from 2000 to 2015, and became the insignificant driver in 2020.

Terrain, vegetation and climate had a consistent positive influence on SR, with terrain having the greatest effect evidenced by path coefficients over 0.5. The impact of terrain was relatively strong, but its intensity had an overall downward trend over time, gradually weakening from 0.670 in 2000 to 0.521 in 2020. Human activity had a slight positive influence on SR from 2000 to 2015, but it was negative in 2020. Soil had a slightly negative impact on SR in 2015, with a path coefficient of −0.036. Human activities, climate and vegetation had a consistent negative influence on CS, which was the most affected by human activities. The effect of terrain on CS was slightly positive in 2010 with a path coefficient of 0.011, but not significant in the other years. The effect of soil on CS was positive from 2000 to 2020. The impact of soil increased from 0.162 in 2000 to 0.197 in 2005, declined thereafter to 0.140 in 2015, and finally rose again to 0.180 in 2020.

3.4.2 Interaction mechanism among the driving factors

Fig.5 shows that the positive effect of human activities on GP and WR was lessened by the indirect effect of terrain. The positive effects of human activities on climate enhanced the positive effect of climate on WR, and lessened the negative effect of climate on GP in 2010, 2015 and 2020. Given the indirect effect of climate, the positive effect of soil on GP was reduced, while the negative effect of soil on WR was enhanced. For GP and WR, climate had a positive effect on vegetation, while terrain and soil negatively affected vegetation.

As shown in Fig.5, the negative effect of climate on soil weakened the positive effect of soil on CS and WP. The negative effect of vegetation on CS and WP was reduced by the positive effect of climate on vegetation and increased by the negative effect of soil on vegetation. Human activities lessen the negative effects of climate on CS and WP, and terrain increased the negative effects of human activities on CS and WP.

Human activities enhanced the positive effects of climate on SR (Fig.5). Climate weakened the positive effect of soil on SR except in 2015. The negative effect of terrain on human activities weakened the positive effect of human activities on SR except in 2020. Climate enhanced the positive effect of vegetation on SR, whereas soil weakened it.

4 Discussion

4.1 Comparison with the similar studies related to ecosystem service evaluation

In Table S3, we compare our study to the similar studies to demonstrate the robustness of our results related to ES evaluation. In our study, the quality of GP can be guaranteed due to it was evaluated mainly based on officially published statistical data and NDVI. The amounts of water yield in our study are close to the natural river runoff data reported in Water Resources Bulletin of China^[63]. We found that the spatial distributions of SR in our study are similar to that of Rao et al.^[64] based on the USLE model. Given that our study presents the first evaluation specified to cropland ESs at national scale, we compared our results with the previous at regional scale. The data in Table S3 shows that our results for CS are close to that in Qinghai, China^[65]. Finally, our results show that regions with higher GP have lower WP, which is accordance with the concept that higher GP indicates higher consumption of nitrogen fertilizers, weakening water purification capacity^[66].

4.2 Realizing the importance of cropland related ecosystem services

Our results are consistent with the understanding that there are significant amounts of regulating ESs from cropland which have been generally ignored^[67]. This inadequacy can hamper a more holistic understanding of sustainable crop production. It has been emphasized that cropland related ESs should be considered in decision-making^[68]. Embracing ES-based on knowledge can provide more systematic insights rather than solely relying on crop production inputs, grain outputs and environmental impacts. The results of our study, which reveal the dynamic changes of cropland ESs in China, can provide important information for enhancing sustainability. For example, the most substantial decreases of ESs occurred in SC. Meanwhile, there were zones of aggregation for some important ESs. Therefore, targeted efforts should be undertaken to solve this dilemma.

Management practices in China generally place emphasis more on the productive functions of cropland, while ignoring its comprehensive role in ESs, and the related relationships as well as driving mechanism, which is critical to inform sustainable development strategies. For example, our study showed that regions with higher GP often have lower WP in HHHP and NCP. Nitrogen management strategies should be optimized, and measures such as precision fertilizer application and enhanced soil nutrient monitoring can effectively reduce nutrient losses while improving crop yields and quality^[69]. The National High-Standard Farmland Construction Plan mainly focuses on improving the production capacity of agricultural products, but neglects comprehensive consideration from the perspective of ESs. Therefore, future versions of this plan should attach importance to the specific process of water purification, analyze the complex relationships between ESs, and promote the farmland construction to balance production and ecology. To address the decline of cropland ESs in some regions, governments should introduce relevant policies to encourage farmers to adopt sustainable agricultural practices, forming a virtuous cycle of food security and environmental sustainability^[70]. Therefore, it is necessary to protect and appropriately use cropland, reduce the negative impact of agricultural practices on the ecosystem, and leverage positive ES functions^[71,^72].

Our study demonstrates that the high-intensity regions of GP were concentrated in HHHP and MLYP, due to their relatively high proportion of cropland area. The high-intensity regions expanded toward NCP, where the largest expansion of cropland area in China was found in this area^[73]. In addition, the NCP and HHHP were the low-intensity regions of SR, which is due to expansion and inappropriate use of cropland^[74,^75]. In regions of GP aggregation, emphasis should be placed on enhancing measures, such as cropland protection and agricultural technology promotion, to ensure stable production while maintaining other ecological functions. The distribution of WR was characterized by the high-intensity south and low-intensity north pattern, paralleling precipitation. This is because WR is largely affected by the climate^[76,^77]. Our results showed that the spatial distribution of WP, with a high-intensity north and low-intensity south pattern, could stem from the loss of nitrogen and phosphorus from cropland causing water pollution^[78]. Fertilizer application can enhance crop productivity^[79], but inappropriate use may exacerbate pollution.

Cropland made significant contributions to the achievement of SDG 2, which aims to achieve zero hunger. Cropland, as the carrier of grain production, ensures food supply of humanity^[72]. Additionally, cropland is linked to poverty eradication (SDG 1), and effective use of cropland and increased crop yields can elevate farmer incomes^[80]. China has given considerable attention to improving the crop production and reducing the related environmental impacts. For example, the farmland redlines policy and high-quality cropland construction, which were also adopted to ensure food security, focused on quantity of rural land. The Cultivated Land Quality Program principally aims to address soil quantity decrease, nutrient imbalances and pollution. However, these policies seldom acknowledge the important role of ESs from cropland. Fortunately, the Chinese government recently requires to enhance the ability of cropland to sequestrate carbon, which not only aids in addressing climate change but also indirectly promotes crop growth, providing support for achieving sustainable development goals related to climate action (SDG 13). We argue that well understanding of the ecosystem functions behind the ESs that affect crop production can also improve social benefits of rural people, especially under China’s rapid urbanization. For example, integrating CS from cropland into China’s carbon trading market is one of the potential ways to increase the farmer incomes.

4.3 Understanding the direct effects of factors on cropland ecosystem services

The results of our PLS-SEM analysis show that terrain and vegetation had the largest negative and positive impacts on GP, respectively. However, these results are consistent with some previous studies. For example, Jiang et al.^[66] found cropland in high altitude had relatively low GP due to the low soil quality and difficulties of cropland management. However, Hu et al.^[25] reported that vegetation had positive impacts on GP in Jiangxi Province. It is not surprising that water yield is positively and largely impacted by climate^[81]. Zhao et al.^[82] found that, among the climate observed variables, precipitation and evapotranspiration had large impacts on WR which is consistent with our findings. Our results also reveal the negative impacts of climate and terrain on WP, which can be attributed to the fact that climate change can modify the conditions and processes of rivers, and thereby influence water quality^[83], whereas flows of pollutants are determined by the elevation within terrain^[66]. Terrain was also found to have the largest impact on SR in our study, which is consistent with Rao et al.^[64] and Zhu et al.^[84]. These researchers demonstrated that terrain affects SR by transforming sediment runoff characteristics and transport processes. It is widely agreed that more intensive agriculture production activities can decrease CS in cropland^[85]. Human activities, which include HFP, share of paddy land and POP in this study, reflect the intensity of agriculture production. Therefore, it can be inferred that human activities can negatively impact cropland-related CS.

It should be noted that the effects of drivers on cropland ESs were not always maintained positive or negative throughout the studied period. There in our results demonstrated that the effects of climate on GP changed from positive in 2000 and 2005 to negative in 2010, 2015, and 2020. This may be due to the fact that the structure of the agrometeorological phenomena in China changed substantively after 2010^[86]. There are considerable studies that have aimed to reveal the complex impacts of climate on crop production^[87,^88]. However, their results lead to divergent conclusions. While our study provides evidence of dynamic changes of the climate effects, which can also inform the construction of climate smart agriculture.

4.4 Toward synergistic management of cropland ecosystem services

The indirect effects of factors and the interaction among the studied factors were examined based our proposed PLS-SEM model. Comparing with the studies solely focusing on the direct effects of drivers which may lack the effective solutions to sustainability challenges^[89,^90], our results demonstrated that there are complex interactions among drivers, which can inform more systematic strategies toward sustainable cropland management. For example, we found that vegetation had both positive effects on GP and SR, and negative effects of soil on vegetation were also found. It is widely acknowledged that soil types with varying levels of fertility and organic matter can have significant impacts on vegetation growth^[91,^92], however, the interactions among soil-vegetation-GP/SR should also be considered. Also, we found that the negative effect of climate on soil reduced the positive effect of soil on CS and WP. Consequently, we highlight the importance of systematic management of cropland due to the fact that the related drivers and environmental flows are generally complexly interconnected in coupled human and natural systems^[26]. There was one study that used a PLS-SEM model to study cropland ESs in Jiangxi, which aggregated four ESs into the cropland ES value based on deemed monetary value. This aggregation approach may indicate important implications for stakeholders related to economic values, however, still hinder our understanding on the interactions among ESs and further the co-benefit strategies. For example, our study demonstrated that trade-offs with CS-WR and WR-WP, whereas the human activities had the different effects on these ESs, which provide important information on synergistic governance.

4.5 Raising targeted measures for improving ecosystem services from cropland

There are many opportunities to use our findings to guide ES-based management strategies toward sustainable development. Actually, governments in China have raised policy interventions to enhance ecosystem function, such as gross ecosystem product accounting^[9], ecological redline policy^[5], and large-scale programs^[4]. However, there still lacks of targeted national measurement related to cropland ESs.

In China, the livelihoods of farmers rely heavily on the ecosystem products but ecosystems in rural areas are generally fragile^[93]. More recently, great efforts have been made to end poverty in China. The related policies, however, generally focus on implementing rural infrastructure construction, creating employment opportunities and special transfer payment, which may aggravate the imbalance between the supply of and demand for ESs. This dual challenge can be solved by adopting ES-based approaches. Some successful regional measurements enacted in China, such as payment for ES schemes, in which the beneficiaries directly pay producers for the delivery of ESs. For example, the PLDL program in Beijing, which improved the local water quantity and quality with economic gain for both ES beneficiaries and providers^[6].

There are prerequisites for using ES-based solutions. ES monitoring systems should be first equipped to gather accurate information. Also, a coordinated vision for cropland ES design and implementation should be raised to avoid separate policies and government agencies for separate ESs. Market-based instruments, such as a trading market, can be adopted to strengthen cross regional cooperation between areas with oversupply and shortage of ESs due to the fact that most ESs are typically public goods. Our results also highlight that there would be trade-offs at times when improving cropland related ESs. Therefore, adopting ES-based approaches should simultaneously consider interactions among multiple sectors to reduce the risks of reducing one problem while exacerbating another.

4.6 Limitations and prospects

Despite its merits there are still some limitations in this study. First, this study focused on five important ESs, while other ESs from cropland, such as crop pollination and cultural services, are ignored, mainly due to data deficiencies. This may limit our understanding of a more holist picture of cropland ESs. However, our study still provides a scientific framework to understand the driving mechanism of cropland ESs. Also, some limitations are still inherent in the model applied. For example, the InVEST model is sensitive to data changes and simplifies the process, which is detailed in its user guide^[41]. The GWR model runs with a uniform bandwidth setting for all variables, which may lead to bias and noise in the model^[94]. The indicators available for assessing the overall effectiveness of PLS-SEM are limited, which may hinder the understanding of the fit of the model. Finally, the results in this study may not reflect the regional variation due to that some localized data were not available.

There are many directions that future studies could usefully take. For example, ecosystem disservices from cropland (e.g., carbon emissions, competition for water and pollution emissions) can also be combined into our proposed framework, especially given that considerable environmental impacts have been induced from cropland. Also, the interaction between cropland ESs and specific agricultural practices can be examined further to provide more practical implications. By optimizing farming systems and implementing appropriate fertilizer and irrigation application, cropland ESs can be improved. In response to intensifying climate change, efforts can be made to examine and establish management models for cropland ecosystems that adapt to such changes. For example, linking future development scenarios (such as SSP and RCP) can provide insightful results for sustainable crop management. By exploring interdisciplinary research pathways, theoretical support can be provided for developing cropland protection measures, which in turn promotes the positive functioning of cropland ESs. Finally, exploring the related beneficiaries and providers of cropland ESs at the multiple scales can build the foundation to enact payments for ES schemes^[95].

5 Conclusions

Bringing consideration of ESs into cropland management has been recognized to move toward sustainable development, understanding the driving mechanism of cropland ESs can provide systematic strategies of realizing co-benefits. This work uses the world’s largest grain producer, China, as a case study. The dynamic changes and driving mechanism of the ESs provided by cropland in China were revealed by integrating InVEST model, a spatial analysis tool and PLS-SEM model. The results indicate that SR and WP decreased at national scale, and CS-WR and WR-WP had persistent trade-offs relationships. Special attention should be given to SC given that there were high-intensity areas for some important ESs and substantial ESs decline. WR and WP were most directly impacted by climate, whereas vegetation, terrain and human activities had the most significant direct effects on GP, SR and CS, respectively. The complex interaction mechanisms among driving factors found highlight the need for synergy management of cropland ESs. Finally, it is concluded that there are opportunities for development of sustainable cropland ES management.

References

[1] Power A G. Ecosystem services and agriculture: Tradeoffs and synergies.Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences, 2010, 365(1554): 2959–2971

[2] Cumming G S, Buerkert A, Hoffmann E M, Schlecht E, von Cramon-Taubadel S, Tscharntke T. Implications of agricultural transitions and urbanization for ecosystem services.Nature, 2014, 515(7525): 50–57

[3] Wang S T, Bai X M, Zhang X L, Reis S, Chen D L, Xu J M, Gu B J. Urbanization can benefit agricultural production with large-scale farming in China.Nature Food, 2021, 2(3): 183–191

[4] Bryan B A, Gao L, Ye Y Q, Sun X F, Connor J D, Crossman N D, Stafford-Smith M, Wu J G, He C Y, Yu D Y, Liu Z F, Li A, Huang Q X, Ren H, Deng X Z, Zheng H, Niu J M, Han G D, Hou X Y. China’s response to a national land-system sustainability emergency.Nature, 2018, 559(7713): 193–204

[5] Bai Y, Wong C P, Jiang B, Hughes A C, Wang M, Wang Q. Developing China’s Ecological Redline Policy using ecosystem services assessments for land use planning.Nature Communications, 2018, 9(1): 3034

[6] Zheng H, Robinson B E, Liang Y C, Polasky S, Ma D C, Wang F C, Ruckelshaus M, Ouyang Z Y, Daily G C. Benefits, costs, and livelihood implications of a regional payment for ecosystem service program.Proceedings of the National Academy of Sciences of the United States of America, 2013, 110(41): 16681–16686

[7] Bateman I J, Harwood A R, Mace G M, Watson R T, Abson D J, Andrews B, Binner A, Crowe A, Day B H, Dugdale S, Fezzi C, Foden J, Hadley D, Haines-Young R, Hulme M, Kontoleon A, Lovett A A, Munday P, Pascual U, Paterson J, Perino G, Sen A, Siriwardena G, van Soest D, Termansen M. Bringing ecosystem services into economic decision-making: land use in the United Kingdom.Science, 2013, 341(6141): 45–50

[8] Goldstein J H, Caldarone G, Duarte T K, Ennaanay D, Hannahs N, Mendoza G, Polasky S, Wolny S, Daily G C. Integrating ecosystem-service tradeoffs into land-use decisions.Proceedings of the National Academy of Sciences of the United States of America, 2012, 109(19): 7565–7570

[9] Ouyang Z Y, Song C S, Zheng H, Polasky S, Xiao Y, Bateman I J, Liu J G, Ruckelshaus M, Shi F Q, Xiao Y, Xu W H, Zou Z Y, Daily G C. Using gross ecosystem product (GEP) to value nature in decision making.Proceedings of the National Academy of Sciences of the United States of America, 2020, 117(25): 14593–14601

[10] Xu W H, Xiao Y, Zhang J J, Yang W, Zhang L, Hull V, Wang Z, Zheng H, Liu J G, Polasky S, Jiang L, Xiao Y, Shi X W, Rao E M, Lu F, Wang X K, Daily G C, Ouyang Z Y. Strengthening protected areas for biodiversity and ecosystem services in China.Proceedings of the National Academy of Sciences of the United States of America, 2017, 114(7): 1601–1606

[11] Zheng H, Wang L J, Peng W J, Zhang C P, Li C, Robinson B E, Wu X C, Kong L Q, Li R N, Xiao Y, Xu W H, Ouyang Z Y, Daily G C. Realizing the values of natural capital for inclusive, sustainable development: Informing China’s new ecological development strategy.Proceedings of the National Academy of Sciences of the United States of America, 2019, 116(17): 8623–8628

[12] Zhang W, Ricketts T H, Kremen C, Carney K, Swinton S M. Ecosystem services and dis-services to agriculture.Ecological Economics, 2007, 64(2): 253–260

[13] Shah S M, Liu G Y, Yang Q, Wang X Q, Casazza M, Agostinho F, Lombardi G V, Giannetti B F. Emergy-based valuation of agriculture ecosystem services and dis-services.Journal of Cleaner Production, 2019, 239: 118019

[14] Zabala J A, Martínez-Paz J M, Alcon F. A comprehensive approach for agroecosystem services and disservices valuation.Science of the Total Environment, 2021, 768: 144859

[15] Hu Y N, Peng J, Liu Y X, Tian L. Integrating ecosystem services trade-offs with paddy land-to-dry land decisions: a scenario approach in Erhai Lake Basin, southwest China.Science of the Total Environment, 2018, 625: 849–860

[16] Zheng H N, Peng J, Qiu S J, Xu Z H, Zhou F, Xia P, Adalibieke W. Distinguishing the impacts of land use change in intensity and type on ecosystem services trade-offs.Journal of Environmental Management, 2022, 316: 115206

[17] Tian Y L, Xu D X, Song J X, Guo J H, You X J, Jiang Y. Impacts of land use changes on ecosystem services at different elevations in an ecological function area, northern China.Ecological Indicators, 2022, 140: 109003

[18] Ma S, Qiao Y P, Jiang J, Wang L J, Zhang J C. Incorporating the implementation intensity of returning farmland to lakes into policymaking and ecosystem management: a case study of the Jianghuai Ecological Economic Zone, China.Journal of Cleaner Production, 2021, 306: 127284

[19] Zhuang M F, Geng Y, Pan H Y, Wu F, Wu D. Ecological and socioeconomic impacts of payments for ecosystem services—A Chinese garlic farm case.Journal of Cleaner Production, 2021, 285: 124866

[20] Zheng W W, Li S M, Ke X L, Li X Y, Zhang B B. The impacts of cropland balance policy on habitat quality in China: a multiscale administrative perspective.Journal of Environmental Management, 2022, 323: 116182

[21] Martin E A, Dainese M, Clough Y, Báldi A, Bommarco R, Gagic V, Garratt M P D, Holzschuh A, Kleijn D, Kovács-Hostyánszki A, Marini L, Potts S G, Smith H G, Al Hassan D, Albrecht M, Andersson G K S, Asís J D, Aviron S, Balzan M V, Baños-Picón L, Bartomeus I, Batáry P, Burel F, Caballero-López B, Concepción E D, Coudrain V, Dänhardt J, Diaz M, Diekötter T, Dormann C F, Duflot R, Entling M H, Farwig N, Fischer C, Frank T, Garibaldi L A, Hermann J, Herzog F, Inclán D, Jacot K, Jauker F, Jeanneret P, Kaiser M, Krauss J, Le Féon V, Marshall J, Moonen A C, Moreno G, Riedinger V, Rundlöf M, Rusch A, Scheper J, Schneider G, Schüepp C, Stutz S, Sutter L, Tamburini G, Thies C, Tormos J, Tscharntke T, Tschumi M, Uzman D, Wagner C, Zubair-Anjum M, Steffan-Dewenter I. The interplay of landscape composition and configuration: new pathways to manage functional biodiversity and agroecosystem services across Europe.Ecology Letters, 2019, 22(7): 1083–1094

[22] Jiang L, Wang Z Z, Zuo Q T, Du H H. Simulating the impact of land use change on ecosystem services in agricultural production areas with multiple scenarios considering ecosystem service richness.Journal of Cleaner Production, 2023, 397: 136485

[23] Renard D, Tilman D. National food production stabilized by crop diversity.Nature, 2019, 571(7764): 257–260

[24] Su Y, Gabrielle B, Makowski D. The impact of climate change on the productivity of conservation agriculture.Nature Climate Change, 2021, 11(7): 628–633

[25] Hu B F, Xie M D, He R B, Shi Z, Zhou Y, Ni H J, Li H Y. Spatio-temporal evolution of cropland ecosystem services value and its spatially varying dominate over the past two decades: a case study in Jiangxi Province in southern China.Journal of Cleaner Production, 2023, 427: 139228

[26] Liu J G, Mooney H, Hull V, Davis S J, Gaskell J, Hertel T, Lubchenco J, Seto K C, Gleick P, Kremen C, Li S X. Systems integration for global sustainability.Science, 2015, 347(6225): 1258832

[27] Liu Q, Qiao J J, Li M J, Huang M J. Spatiotemporal heterogeneity of ecosystem service interactions and their drivers at different spatial scales in the Yellow River Basin.Science of the Total Environment, 2024, 908: 168486

[28] Li W S, Wang L Q, Yang X, Liang T, Zhang Q, Liao X Y, White J R, Rinklebe J. Interactive influences of meteorological and socioeconomic factors on ecosystem service values in a river basin with different geomorphic features.Science of the Total Environment, 2022, 829: 154595

[29] Xu X L, Liu J Y, Zhang S W, Li R D, Yan C Z, Wu S X. China multi-period land use and land cover change remote sensing monitoring data set (CNLUCC). Resource and Environmental Science Data Registration and Publishing System, 2018

[30] Xu X L. Monthly spatial interpolation dataset of Chinese meteorological elements. Resource and Environmental Science Data Registration and Publication System, 2022

[31] Peng S Z. 1-km monthly potential evapotranspiration dataset in China (1901–2022). National Tibetan Plateau Data Center, 2022

[32] Peng S Z. 1-km monthly mean temperature dataset for china (1901–2023). National Tibetan Plateau Data Center, 2019

[33] Food, Agriculture Organization of the United Nations (FAO), International Institute for Applied Systems Analysis (IIASA), International Soil Reference, Information Center (ISRIC), Institute of Soil Science-Chinese Academy of Sciences (ISS-CAS), Joint Research Centre (JRC). (FAO, IIASA, ISRIC, ISS-CAS, JRC). Harmonized World Soil Database Version 1.1. FAO, IIASA, ISRIC, ISS-CAS, JRC, 2009

[34] Didan K. MODIS/Terra Vegetation Indices Monthly L3 Global 1km SIN Grid V006. NASA EOSDIS Land Processes DAAC, 2015

[35] Didan K. MODIS/Terra Vegetation Indices 16-Day L3 Global 250m SIN Grid V006. NASA EOSDIS Land Processes DAAC, 2015

[36] Myneni R, Knyazikhin Y, Park T. MODIS/Terra Leaf Area Index/FPAR 8-Day L4 Global 500m SIN Grid V006. NASA EOSDIS Land Processes DAAC, 2015

[37] Xu X L. Chinese GDP Spatial Distribution Kilometer Grid Dataset. Resource and Environmental Science Data Registration and Publishing System, 2017

[38] Xu X L. China Population Spatial Distribution Kilometer Grid Dataset. Resource and Environmental Science Data Registration and Publishing System, 2017

[39] Mu H W, Li X C, Wen Y N, Huang J X, Du P J, Su W, Miao S X, Geng M Q. A global record of annual terrestrial Human Footprint dataset from 2000 to 2018.Scientific Data, 2022, 9(1): 176

[40] Zhong X Y, Yan Q W, Li G E. Development of long time series nighttime light dataset of China (2000–2020). Journal of Global Change Data & Discovery, 2022, 6(3): 416−424 (in Chinese and English)

[41] Natural Capital Project. InVEST 3.13.0 User’s Guide. Stanford University, University of Minnesota, Chinese Academy of Sciences, The Nature Conservancy, World Wildlife Fund, and Stockholm Resilience Centre, 2022

[42] Zhang D, Min Q W, Liu M C, Cheng S K. Ecosystem service tradeoff between traditional and modern agriculture: a case study in Congjiang County, Guizhou Province, China.Frontiers of Environmental Science & Engineering, 2012, 6(5): 743–752

[43] Kuri F, Murwira A, Murwira K S, Masocha M. Predicting maize yield in Zimbabwe using dry dekads derived fromremotely sensed vegetation condition index.International Journal of Applied Earth Observation and Geoinformation, 2014, 33: 39–46

[44] Ouyang Z Y, Zheng H, Xiao Y, Polasky S, Liu J G, Xu W H, Wang Q, Zhang L, Xiao Y, Rao E M, Jiang L, Lu F, Wang X K, Yang G B, Gong S H, Wu B F, Zeng Y, Yang W, Daily G C. Improvements in ecosystem services from investments in natural capital.Science, 2016, 352(6292): 1455–1459

[45] Sun X, Li F. Spatiotemporal assessment and trade-offs of multiple ecosystem services based on land use changes in Zengcheng, China.Science of the Total Environment, 2017, 609: 1569–1581

[46] Liu Z T, Wu R, Chen Y X, Fang C L, Wang S J. Factors of ecosystem service values in a fast-developing region in China: insights from the joint impacts of human activities and natural conditions.Journal of Cleaner Production, 2021, 297: 126588

[47] Getis A, Ord J K. The analysis of spatial association by use of distance statistics.Geographical Analysis, 1992, 24(3): 189–206

[48] Ndong G O, Villerd J, Cousin I, Therond O. Using a multivariate regression tree to analyze trade-offs between ecosystem services: application to the main cropping area in France.Science of the Total Environment, 2021, 764: 142815

[49] Li G Z, Cai J. Spatial and temporal differentiation of mountain ecosystem service trade-offs and synergies: a case study of Jieshi Mountain, China.Sustainability, 2022, 14(8): 4652

[50] Gao J B, Zuo L Y. Revealing ecosystem services relationships and their driving factors for five basins of Beijing.Journal of Geographical Sciences, 2021, 31(1): 111–129

[51] Gollini I, Lu B B, Charlton M, Brunsdon C, Harris P. GWmodel: an R package for exploring spatial heterogeneity using geographically weighted models.Journal of Statistical Software, 2015, 63(17): 1–50

[52] Schuberth F, Rosseel Y, Rönkkö M, Trinchera L, Kline R B, Henseler J. Structural parameters under partial least squares and covariance-based structural equation modeling: a comment on Yuan and Deng (2021).Structural Equation Modeling, 2023, 30(3): 339–345

[53] Huang Y, Wei C L, Yan H, Hao R J. Ability modeling for spacecraft control system based on partial least square structural equation model approach. Chinese Space Science and Technology, 2024, 44(2): 98−108 (in Chinese)

[54] Hair J F, Risher J J, Sarstedt M, Ringle C M. When to use and how to report the results of PLS-SEM.European Business Review, 2019, 31(1): 2–24

[55] Yuan J, Li R, Huang K. Driving factors of the variation of ecosystem service and the trade-off and synergistic relationships in typical karst basin.Ecological Indicators, 2022, 142: 109253

[56] Xia H, Yuan S F, Prishchepov A V. Spatial-temporal heterogeneity of ecosystem service interactions and their social-ecological drivers: implications for spatial planning and management.Resources, Conservation and Recycling, 2023, 189: 106767

[57] Peng H J, Zhang X S, Ren W, He J. Spatial pattern and driving factors of cropland ecosystem services in a major grain-producing region: A production-living-ecology perspective.Ecological Indicators, 2023, 155: 111024

[58] Hair J F, Sarstedt M, Ringle C M, Mena J A. An assessment of the use of partial least squares structural equation modeling in marketing research.Journal of the Academy of Marketing Science, 2012, 40(3): 414–433

[59] Nunnally J C, Bernstein I H. Psychometric Theory. New York: McGraw-Hill, 1994

[60] Fornell C, Larcker D F. Evaluating structural equation models with unobservable variables and measurement error.Journal of Marketing Research, 1981, 18(1): 39–50

[61] Stone M. Cross-validatory choice and assessment of statistical predictions.Journal of the Royal Statistical Society. Series B, Statistical Methodology, 1974, 36(2): 111–133

[62] Tenenhaus M, Amato S, Vinzi V E. A global goodness-of-fit index for PLS structural equation modelling.The XLII SIS Scientific Meeting, 2004, 1(2): 739–742

[63] Ministry of Water Resources (MWR) of the People’s Republic of China. Water Resources Bulletin of China. Available at MWR website on July 9, 2021

[64] Rao E M, Ouyang Z Y, Yu X X, Xiao Y. Spatial patterns and impacts of soil conservation service in China.Geomorphology, 2014, 207: 64–70

[65] Yan K, Wang W F, Li Y H, Wang X F, Jin J X, Jiang J, Yang H Q, Wang L J. Identifying priority conservation areas based on ecosystem services change driven by Natural Forest Protection Project in Qinghai Province, China.Journal of Cleaner Production, 2022, 362: 132453

[66] Jiang S, Cheng X Y, Yu S Y, Zhang H B, Xu Z H, Peng J. Elevation dependency of ecosystem services supply efficiency in great lake watershed.Journal of Environmental Management, 2022, 318: 115476

[67] Halperin S, Castro A J, Quintas-Soriano C, Brandt J S. Assessing high quality agricultural lands through the ecosystem services lens: Insights from a rapidly urbanizing agricultural region in the western United States.Agriculture, Ecosystems & Environment, 2023, 349: 108435

[68] Gawith D, Hodge I. Focus rural land policies on ecosystem services, not agriculture.Nature Ecology & Evolution, 2019, 3(8): 1136–1139

[69] Liu Y Z, Zhuang M H, Liang X, Lam S K, Chen D L, Malik A, Li M Y, Lenzen M, Zhang L Y, Zhang R, Zhang L X, Hao Y. Localized nitrogen management strategies can halve fertilizer use in Chinese staple crop production.Nature Food, 2024, 5(10): 825–835

[70] Cai S Y, Zhao X, Pittelkow C M, Fan M S, Zhang X, Yan X Y. Optimal nitrogen rate strategy for sustainable rice production in China.Nature, 2023, 615(7950): 73–79

[71] Wang J Y, Li G Y, Fang X Q, Su M Y, Cao Y. Exploring the connotation and research framework of cultivated land ecosystem services from the perspective of social-ecological systems. Acta Ecologica Sinica, 2024, 44(15): 6881−6891 (in Chinese)

[72] Zhou Y, Li X H, Liu Y S. Cultivated land protection and rational use in China.Land Use Policy, 2021, 106: 105454

[73] Zhong H M, Liu Z J, Wang J Y. Understanding impacts of cropland pattern dynamics on grain production in China: an integrated analysis by fusing statistical data and satellite-observed data.Journal of Environmental Management, 2022, 313: 114988

[74] Tang Q, Hua L, Cao Y F, Jiang L, Cai C F. Human activities are the key driver of water erosion changes in northeastern China.Land Degradation & Development, 2024, 35(1): 62–75

[75] Hu X Y, Dong C, Zhang Y. Dynamic evolution of the ecological footprint of arable land in the Yellow and Huaihai Main grain producing area based on structural equation modeling and analysis of driving factors.Ecological Informatics, 2024, 82: 102720

[76] Hou J X, Pan H H, Du Z Q, Wu Z T, Zhang H. Spatiotemporal analysis of water ecosystem services of the Yellow River Basin in Shanxi Province. Arid Land Geography, 2024, 47(6): 1047-1060 (in Chinese)

[77] Zhou L Y, Cui W L, Yang F. Spatiotemporal variations and driving forces analysis of ecosystem water conservation in coastal areas of China.Ecological Indicators, 2024, 162: 112019

[78] Zhao L Q, Fu D G, Li T, Yuan X Q, Wang S C, Liu C E, Duan C Q. Temporal variations of N and P losses via surface runoff from Chinese farmland after fertilisation.Soil & Tillage Research, 2025, 246: 106338

[79] Gu B J, Zhang X M, Lam S K, Yu Y L, van Grinsven H J M, Zhang S H, Wang X X, Bodirsky B L, Wang S T, Duan J K, Ren C C, Bouwman L, de Vries W, Xu J M, Sutton M A, Chen D L. Cost-effective mitigation of nitrogen pollution from global croplands.Nature, 2023, 613(7942): 77–84

[80] Tian F Y, Wu B F, Zeng H W, Watmough G R, Zhang M, Li Y R. Detecting the linkage between arable land use and poverty using machine learning methods at global perspective.Geography and Sustainability, 2022, 3(1): 7–20

[81] He L J, Zheng B F, Wan W, Xie Z Y, Song X, Zhu J Q, Liu Z. Trade-off and synergy of ecosystem services in the Yangtze River Economic Belt and its driving factors. Environmental Sciences, 2024, 45(6): 3318−3328 (in Chinese)

[82] Zhao G L, Tian S M, Liang S, Jing Y C, Chen R X, Wang W W, Han B. Dynamic evolution trend and driving mechanisms of water conservation in the Yellow River Basin, China.Scientific Reports, 2024, 14(1): 26304

[83] Li L, Knapp J L A, Lintern A, Ng G H C, Perdrial J, Sullivan P L, Zhi W. River water quality shaped by land–river connectivity in a changing climate.Nature Climate Change, 2024, 14(3): 225–237

[84] Zhu M Y, He W M, Zhang Q F, Xiong Y Z, Tan S D, He H M. Spatial and temporal characteristics of soil conservation service in the area of the upper and middle of the Yellow River, China.Heliyon, 2019, 5(12): e02985

[85] Intergovernmental Panel on Climate Change (IPCC). Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge, UK: Cambridge University Press, Cambridge, 2013

[86] Zhao Y J, Zheng R H, Zheng F L, Zhong K Y, Fu J X, Zhang J Q, Flanagan D C, Xu X M, Li Z. Spatiotemporal distribution of agrometeorological disasters in China and its impact on grain yield under climate change.International Journal of Disaster Risk Reduction, 2023, 95: 103823

[87] Mirón I J, Linares C, Díaz J. The influence of climate change on food production and food safety.Environmental Research, 2023, 216(P3): 114674

[88] Kukal M S, Irmak S. Climate-driven crop yield and yield variability and climate change impacts on the US Great Plains agricultural production.Scientific Reports, 2018, 8(1): 3450

[89] Mina M, Bugmann H, Cordonnier T, Irauschek F, Klopcic M, Pardos M, Cailleret M. Future ecosystem services from European mountain forests under climate change.Journal of Applied Ecology, 2017, 54(2): 389–401

[90] Li D L, Cao W F, Dou Y H, Wu S Y, Liu J G, Li S C. Non-linear effects of natural and anthropogenic drivers on ecosystem services: Integrating thresholds into conservation planning.Journal of Environmental Management, 2022, 321: 116047

[91] Qiu H H, Zhang J Y, Han H R, Cheng X Q, Kang F F. Study on the impact of vegetation change on ecosystem services in the Loess Plateau, China.Ecological Indicators, 2023, 154: 110812

[92] Li W W, Yan D H, Weng B S, Lai Y Q, Zhu L, Qin T L, Dong Z Y, Bi W X. Nonlinear effects of surface soil moisture changes on vegetation greenness over the Tibetan Plateau.Remote Sensing of Environment, 2024, 302: 113971

[93] Wang Y, Bilsborrow R E, Zhang Q, Li J F, Song C H. Effects of payment for ecosystem services and agricultural subsidy programs on rural household land use decisions in China: Synergy or trade-off.Land Use Policy, 2019, 81: 785–801

[94] Huang Z W, Li S Y, Peng Y H, Gao F. Spatial non-stationarity of influencing factors of China’s county economic development base on a multiscale geographically weighted regression model.ISPRS International Journal of Geo-Information, 2023, 12(3): 109

[95] Liu J G, Yang W. Integrated assessments of payments for ecosystem services programs.Proceedings of the National Academy of Sciences of the United States of America, 2013, 110(41): 16297–16298

RIGHTS & PERMISSIONS

The Author(s) 2025. Published by Higher Education Press. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0)