Land-use dynamics and multi-scenario projections in a mountainous ecological barrier: ecosystem service value responses and spatial sensitivity

1 Introduction

Land-use change alters ecosystem composition, spatial configuration, and the provision of ecosystem services (Foley et al., 2005; Turner et al., 2007). This influence is particularly pronounced in mountainous ecological barrier regions, where transitions among land-use categories may modify not only the spatial configuration of land use but also the continuity of forest landscapes, the composition of ecological land, and the intensity of human disturbance. These changes can subsequently affect key ecosystem services, including water conservation, climate regulation, soil retention, biodiversity conservation, and landscape maintenance (Lambin et al., 2003; Foley et al., 2005; De Groot et al., 2010; Costanza et al., 2014; Xie et al., 2017). In mountainous regions, ecological land, including forestland, grassland, wetlands, and water bodies, often constitutes a major source of ecosystem service value. Changes in forest extent, increasing landscape fragmentation, and disturbance along forest edges may alter ecological connectivity and further influence regional ecological security (Bonan, 2008; Pan et al., 2011; Gamfeldt et al., 2013; Watson et al., 2018). Scenario-based land-use projections help evaluate how alternative policy priorities alter land allocation and ecosystem service value (ESV) (Liang et al., 2021; Lou et al., 2022; Qi et al., 2024). This issue is particularly relevant in mountainous ecological barriers, where steep topographic gradients and concentrated human activity can produce marked differences in ESV responses between continuous forested mountain areas and intensively used valleys or piedmont zones (Li et al., 2021; Zhang et al., 2024).

Land-use scenario modelling generally involves two related tasks: estimating the future demand for each land-use category and allocating that demand spatially. Markov models estimate category-level demand from historical transition probabilities, whereas the PLUS model allocates projected demand by reproducing patch-based land expansion (Pontius and Malanson, 2005; Liang et al., 2021). Coupling these models with ESV assessment enables comparison of alternative land-use pathways in terms of both changes in regional ESV totals and spatial ecological outcomes (Costanza et al., 1997; Li et al., 2021; Lou et al., 2022; Jiang et al., 2023). Most scenario-based ESV studies report changes in regional totals or map the distribution of high- and low-value areas (Li et al., 2021; Lou et al., 2022; Jiang et al., 2023; Na et al., 2023). These outputs show whether a scenario produces an overall ecological gain or loss and where ESV is concentrated, but they do not indicate whether a relatively limited amount of land conversion produces a disproportionately large local response. This distinction is important in mountainous regions, where broad areas of forestland may remain stable while land conversion is concentrated in narrow river valleys, piedmont zones, and peri-urban areas (Zhang et al., 2024). To address this limitation, local ESV change needs to be evaluated relative to the intensity of land-use conversion within the same spatial unit.

The Shaanxi section of the Qinling Mountains provides a particularly informative setting for examining the contrast between regional stability and localized ecological responses. Its middle- and high-elevation areas are dominated by extensive and relatively continuous forestland, whereas cropland, settlements, transportation corridors, and recent construction-land expansion are concentrated mainly in river valleys, piedmont zones, and other low-elevation areas (Wang et al., 2024; Zhang et al., 2024). This spatial pattern creates a distinctive configuration in which a largely stable, forest-dominated regional matrix coexists with localized but intensive land-use conversion. Consequently, changes in total regional ESV may remain limited even when particular transition zones experience substantial ecological gains or losses. The Qinling Mountains therefore provide a useful case for testing whether aggregate ESV trends adequately represent the ecological consequences of spatially concentrated land conversion. Against this background, this study combines multi-scenario land-use projection with a grid-based sensitivity index to distinguish areas of high ESV from areas where ESV responds strongly to land-use conversion. Unlike conventional ESV mapping, which primarily identifies the spatial distribution of ecosystem-service value, the SI relates local relative ESV change to the proportion of land converted within the same grid cell. It therefore provides a complementary measure of the direction and relative intensity of ESV response to land-use change rather than an indicator of ecological vulnerability. Accordingly, we examined historical land-use dynamics during 2000–2020 and projected land-use patterns to 2040 under four contrasting scenarios (Figure 1). This study addresses three questions:

• What were the main land-use dynamics from 2000 to 2020, and how do projected land-use quantities and spatial patterns differ among the four scenarios in 2040?
• How do the four projected land-use pathways affect total ESV and its spatial distribution?
• Where are ESV responses most sensitive to land-use change, and how do the direction and intensity of these responses differ among scenarios?

2 Study area and data sources

2.1 Study area

The Shaanxi section of the Qinling Mountains, located in central and southern Shaanxi Province, China, was selected as a representative mountainous ecological barrier region (Wang et al., 2024; Chen et al., 2025). Figure 2 shows the geographical location of the study area. The region forms an important geographical boundary between northern and southern China and a major watershed between the Yellow River and Yangtze River basins. It is bordered by the Guanzhong Plain to the north and the Hanjiang River Basin to the south and is characterized by pronounced topographic, climatic, and ecological gradients (Chen et al., 2025). Mountains, hills, and river valleys dominate the landscape, and substantial elevation differences generate pronounced spatial variation in hydrothermal conditions and vegetation. Extensive and relatively continuous forest ecosystems support water conservation, soil retention, biodiversity maintenance, and ecological connectivity, making the region an important component of China’s ecological security system (Wang et al., 2024; Zhang et al., 2024; Chen et al., 2025). Owing to its complex natural environment, the region contains a mosaic of land-use types, including cropland, forestland, grassland, wetlands, water bodies, and construction land. Land use exhibits a clear topographic gradient: forestland is extensively distributed across the middle- and high-elevation mountains, whereas cropland and construction land are concentrated mainly in river valleys, piedmont zones, and other relatively flat areas with greater accessibility and more intensive human activity (Zhang et al., 2024; Chen et al., 2025). The juxtaposition of extensive ecological land with localized and intensive land conversion creates strong spatial contrasts in both baseline ESV and its potential response to land-use change. Land-use pressures also vary considerably within the region. The northern foothills and major river-valley corridors experience greater pressure from urban expansion, infrastructure development, and agricultural production, whereas the middle- and high-elevation mountain areas are dominated by ecological conservation functions and are subject to stricter development constraints (Wang et al., 2024; Zhang et al., 2024; Chen et al., 2025). The region is simultaneously subject to pressures from urban development and policy priorities related to cropland protection, ecological restoration, and nature conservation. These spatial contrasts make the area well suited to examining whether regional ESV trends are consistent with local responses to land conversion (Wang et al., 2024; Zhang et al., 2024).

2.2 Data sources

The analysis used three groups of spatial data. First, GlobeLand30 maps for 2000, 2010, and 2020 provided the land-use time series. Second, elevation, slope, population density, GDP, and proximity to roads, railways, settlements, and water bodies were used in the spatial allocation model. Third, land-use-specific value coefficients were applied in the ESV assessment. Table 1 lists the source, date, original resolution, and role of each dataset. GlobeLand30 maps for 2000, 2010, and 2020 were clipped to the study-area boundary and reclassified into cropland, forestland, grassland, wetland, water bodies, and construction land. The resulting 30-m maps were used to quantify historical change and to prepare the 2010 and 2020 inputs for the Markov–PLUS analysis (Chen et al., 2017; Chen and Chen, 2018). The three temporal land-use datasets were used to characterize historical land-use dynamics during 2000–2020. The 2010 and 2020 land-use datasets were further used to estimate land-use demand with the Markov model, allocate projected land use spatially with the PLUS model, and generate land-use projections for 2040 under four scenarios.

TABLE 1. Source datasets and spatial predictor variables used for land-use projection and ecosystem service value assessment.

All spatial datasets were projected to WGS 1984 UTM Zone 49 N (EPSG:32649), clipped to the study-area boundary, and aligned to the 30-m grid of the reclassified GlobeLand30 land-use dataset. Resampling of coarse-resolution datasets was conducted solely to ensure grid compatibility and did not increase their original effective spatial resolution. Distance-based predictor variables were generated using Euclidean-distance analysis in ArcGIS.

Elevation and slope were derived from a DEM accessed via the Geospatial Data Cloud. Socioeconomic predictors (population density and GDP, each at 1-km resolution) were obtained from the Data Center of the Chinese Academy of Sciences, whereas spatial proximity predictors—Euclidean distances to roads, railways, settlements, and water bodies—were generated from data obtained from the National Geomatics Database and open-access geospatial platforms. ESV coefficients adjusted for regional conditions following Costanza et al. (1997), Xie et al. (2015), Xie et al. (2017) were then applied to estimate service values for the six land-use categories. Table 1 summarizes the complete set of variables and their sources.

All spatial datasets were preprocessed in ArcGIS to ensure consistency in coordinate system, spatial extent, cell size, and raster alignment. All layers were projected to WGS 1984 UTM Zone 49 N (EPSG:32649) and clipped to the study-area boundary. The land-use datasets retained their original 30-m resolution, whereas the continuous predictor variables were resampled and aligned to the same 30-m grid to meet the input requirements of the PLUS model. Resampling the 1-km socioeconomic datasets to 30 m was performed solely to ensure grid compatibility and did not increase their effective spatial resolution. The spatial extent, cell size, number of rows and columns, and raster origin were standardized across all layers to ensure pixel-level alignment (Liang et al., 2021). Elevation and slope were derived from the DEM, and proximity variables representing distance to roads, railways, settlements, and water bodies were calculated using Euclidean-distance analysis. These variables, together with population density and GDP, were used as spatial predictors for land-use allocation in the PLUS model. The same six-category land-use datasets were used throughout the historical analysis, Markov demand projection, PLUS spatial allocation, and ESV assessment.

3 Research methods

3.1 Multi-scenario projection of land-use demand using the Markov model

The Markov model projected the demand for each land-use category in 2040 based on transition probabilities derived from historical land-use changes. The model assumes that the land-use state at a subsequent time is conditional on the current state and the transition probability matrix estimated from the preceding period (Pontius and Malanson, 2005). In this study, the 2010 and 2020 land-use datasets were used to derive the transition probability matrix for six categories: cropland, forestland, grassland, wetland, water bodies, and construction land. Scenario-specific land-use demand was then estimated for 2040 under the natural development, urban development, cropland protection, and ecological protection scenarios. These demand estimates represent conditional scenario projections rather than deterministic forecasts and served as quantitative constraints for spatial allocation in the PLUS model (Liang et al., 2021; Lou et al., 2022).

3.1.1 Scenario design and land-use demand constraints

To examine how land-use composition and spatial patterns may differ under contrasting management priorities, we developed four alternative scenarios: natural development, urban development, cropland protection, and ecological protection. Scenario differences were represented through scenario-specific land-use demand constraints and conversion cost matrices, which defined alternative pathways for construction-land expansion, cropland conservation, and ecological-land protection (Li et al., 2021; Lou et al., 2022; Jiang et al., 2023). Following the policy-oriented scenario logic adopted in recent land-use planning studies (Wu et al., 2026), the scenario assumptions were further linked to major land-management priorities relevant to the Qinling region, including development control, cropland conservation, and ecological protection. These scenarios were designed as simplified representations of contrasting planning orientations rather than direct reproductions of statutory planning zones. Accordingly, ecological conservation redlines, permanent basic farmland, and urban development boundaries served as policy contexts for scenario interpretation rather than as explicit spatial constraint layers in the present model.

3.1.1.1 Natural development scenario

The natural development scenario retained the 2010–2020 transition probabilities without additional constraints on ecological land, cropland, or construction land. It served as the baseline projection (Lou et al., 2022; Qi et al., 2024). In policy terms, this scenario represents a trend-continuation pathway against which the effects of stronger development, cropland-protection, and ecological-protection priorities can be compared.

3.1.1.2 Urban development scenario

Under urban development, construction-land demand was increased, and conversion from cropland, forestland, and grassland to construction land was permitted. Conversion of wetlands and water bodies to construction land remained restricted (Deng et al., 2021; Jiang et al., 2023; Li et al., 2021). This scenario represents a development-oriented planning context that gives greater priority to urban and infrastructure expansion, while the conversion of key ecological land remains constrained by broader territorial spatial planning and ecological-protection objectives.

3.1.1.3 Cropland protection scenario

The cropland protection scenario imposed a minimum cropland target and restricted the conversion of cropland to construction land and other categories. Construction-land growth was also reduced relative to the urban and natural development settings (Lou et al., 2022; Wang et al., 2020). This scenario broadly reflects the farmland-conservation priorities of territorial spatial planning, particularly the policy objective of maintaining cultivated land and protecting permanent basic farmland. The statutory boundaries of permanent basic farmland were not directly imposed as spatial constraints; instead, this policy priority was represented through land-demand and conversion restrictions.

3.1.1.4 Ecological protection scenario

The ecological protection scenario restricted the loss of forestland, grassland, wetlands, and water bodies and permitted selected conversions from cropland or construction land to ecological land (Wang et al., 2024; Chen et al., 2025; Zhang et al., 2024). This scenario broadly represents a conservation-oriented planning context consistent with the objectives of ecological conservation redlines, key ecological function zones, and ecological restoration in mountainous ecological barrier regions. Following the policy-oriented approach of Wu et al. (2026), these policy objectives were translated into stronger constraints on ecological-land loss and greater opportunities for conversion to ecological land, rather than being represented through the direct reproduction of statutory ecological-redline boundaries.

3.1.2 Scenario-specific land-use demand projection

After defining the four scenarios, the Markov model projected the demand for each land-use category in 2040. The 2010 and 2020 land-use datasets were overlaid to derive the transition probability matrix, and the observed land-use composition in 2020 served as the initial state. Scenario-specific demand constraints were then applied to generate the projected quantities of cropland, forestland, grassland, wetland, water bodies, and construction land under the four scenarios. These projected quantities subsequently served as quantitative targets for spatial allocation in the PLUS model. The Markov model can be expressed as follows (Pontius and Malanson, 2005):

S t + 1 = S t P

Where S t denotes the land-use area vector at time t, S t + 1 denotes the projected land-use area vector at the subsequent time step, and P is the transition probability matrix. Because the 2010 and 2020 land-use datasets were separated by 10 years, a 10-year time step was adopted.

The transition probability between land-use categories was calculated as follows:

p i j = A i j ∑ j = 1 n A i j

Where p i j denotes the probability of conversion from land-use category i to category j, A i j is the area converted from category i to category j during the study period, and n is the total number of land-use categories.

The Markov projections were further adjusted using scenario-specific quantity constraints to represent different development priorities, including urban expansion, cropland conservation, and ecological-land protection. The resulting 2040 land-use demand vectors then served as quantitative inputs for PLUS-based spatial allocation. They should be interpreted as conditional projections under predefined policy assumptions rather than deterministic forecasts of future land-use change.

3.1.3 Parameterization of scenario projections

To represent differences in land-use transition rules among the four scenario projections, we developed scenario-specific conversion cost matrices for the 2020–2040 period. In the PLUS model, the conversion cost matrix regulates whether one land-use category can be converted into another during spatial allocation. These matrices incorporated the management priorities associated with each scenario into the spatial allocation of projected land use (Jiang et al., 2023; Li et al., 2021; Liang et al., 2021). Given the dominance of forestland in the study area and the importance of ecological-land continuity for regional ecosystem services, the conversion rules were designed to emphasize differentiated constraints on construction-land expansion, cropland conservation, and ecological-land protection (Jiang et al., 2023; Li et al., 2021; Zhang et al., 2024).

In the conversion cost matrix, each row denotes the source land-use category, and each column denotes the target land-use category. The matrix used binary values of 1 and 0: a value of 1 indicates that conversion between two land-use categories is permitted, whereas a value of 0 indicates that the conversion is prohibited or constrained. The six land-use categories were coded as a–f: a represents cropland, b forestland, c grassland, d wetland, e water bodies, and f construction land. The corresponding conversion cost matrices are presented in Table 2.

TABLE 2. Conversion cost matrices for land-use change under different scenarios during 2020–2040.

a, b, c, d, e, and f represent cropland, forestland, grassland, wetland, water bodies, and construction land, respectively; 1 indicates that conversion is permitted, whereas 0 indicates that conversion is prohibited or constrained.

Table 2 shows clear differences in land-use conversion rules among the four scenarios. Under the natural development scenario, historical land-use transition patterns were largely retained. Conversions among cropland, forestland, and grassland remained possible, while construction land continued to expand in line with the previous growth trend. The urban development scenario prioritized construction-land expansion. In this scenario, cropland, forestland, and grassland were allowed to convert to construction land. However, conversions from wetlands and water bodies to construction land were restricted to limit urban encroachment on ecologically important areas. The cropland protection scenario focused on maintaining cropland area. Conversions from cropland to forestland, grassland, wetland, water bodies, and construction land were strictly limited, whereas conversions from forestland and grassland to cropland were permitted. This setting reflects the priority given to cropland conservation and the maintenance of agricultural land. The ecological protection scenario emphasized the conservation of ecological land, including forestland, grassland, wetlands, and water bodies. Conversions from these ecological land categories to construction land were constrained, while some cropland and construction land were allowed to shift toward ecological land such as forestland and grassland. This setting represents ecological restoration and the expansion of ecological land under the ecological protection scenario (Zhang et al., 2024; Chen et al., 2025). These conversion matrices therefore served as simplified modelling representations of the contrasting policy and planning orientations described in Section 3.1.1. Together with scenario-specific land-use demand, they translated differences in development priorities, cropland conservation, and ecological protection into comparable land-use transition constraints, rather than reproducing the legally defined boundaries of individual policy instruments.

3.2 Spatial allocation of projected land use using the PLUS model

PLUS spatially allocated the Markov-derived land demand. The model related land-use expansion observed between 2010 and 2020 to eight predictor variables to estimate the development probability of each land-use category. The cellular-automata module then allocated the scenario-specific demand subject to the corresponding conversion matrix (Liang et al., 2021). The land expansion analysis module estimates the spatial associations between observed land-use expansion and selected predictor variables. The resulting development probabilities, together with scenario-specific land-use demand and conversion rules, then guide the generation of spatially explicit land-use projections. Compared with conventional cellular automata models, PLUS better reproduces patch formation and the spatial expansion of multiple land-use categories, making it suitable for alternative land-use scenarios (Liu et al., 2017; Luan and Liu, 2022; Zhao and Shao, 2023).

Using the 2010 and 2020 land-use datasets, the PLUS model identified the expansion areas of different land-use categories during 2010–2020. The land expansion analysis module linked the observed expansion of each land-use category during 2010–2020 with selected spatial predictor variables to estimate land-use development probabilities. Considering the mountainous topography of the Shaanxi section of the Qinling Mountains, its ecological conservation function, and the influence of human activities, eight driving variables were selected for land-use spatial simulation: elevation, slope, distance to roads, distance to railways, distance to settlements, distance to water bodies, population density, and GDP (Liu et al., 2017; Liang et al., 2021). These variables represented terrain constraints, accessibility, hydrological proximity, population pressure, and economic intensity. Together, these predictor variables represent three broad dimensions relevant to land-use allocation: terrain conditions, transportation and settlement accessibility, and socioeconomic intensity. These variables were selected because they are spatially explicit, consistently available across the study area, and commonly used in regional land-use simulation to represent major biophysical constraints, accessibility gradients, and socioeconomic pressures. The predictor set captures the principal spatial factors associated with observed land-use expansion during 2010–2020 while maintaining consistency in spatial coverage and resolution. Policy-related differences among the four projections were represented primarily through scenario-specific land-use demand and conversion rules rather than spatially explicit policy layers.

For each scenario, the projected demand vector and conversion matrix were entered into PLUS to generate the 2040 land-use map. Model performance was assessed by hindcasting the 2020 pattern from the 2010 map; comparison with the observed 2020 map yielded an overall accuracy of 92.03% and a Kappa coefficient of 0.780 (Congalton, 1991; Foody, 2002; Pontius and Millones, 2011; Olofsson et al., 2014). These metrics were used only to evaluate the model’s ability to reproduce an observed historical land-use pattern under the calibrated configuration. They were not interpreted as evidence of predictive validity for the 2040 projections, which remain conditional on the historical transition structure, selected predictors, scenario-specific demand, and conversion rules.

3.3 Ecosystem service value assessment

Ecosystem service value (ESV) assessment provides a useful framework for quantifying ecosystem services and evaluating the ecological consequences of land-use change. Because land-use categories differ in their capacity to supply ecosystem services, changes in land area and spatial conversion among categories can alter the magnitude, composition, and spatial distribution of regional ESV. The valuation of ecosystem services was initially developed by Costanza et al. Later, Xie et al. adapted and refined the ecosystem service value equivalent factors according to the characteristics of China’s terrestrial ecosystems, making this approach more suitable for regional-scale assessments of land-use-change-induced ecological effects in China (Costanza et al., 1997; Xie et al., 2015; 2017).

Based on the land-use classification scheme, the ESV assessment included six land-use categories: cropland, forestland, grassland, wetland, water bodies, and construction land. Ecosystem services were grouped into four primary classes: provisioning, regulating, supporting, and cultural services. These classes included 11 secondary service functions: food production, raw material production, water supply, gas regulation, climate regulation, environmental purification, hydrological regulation, soil conservation, nutrient cycling maintenance, biodiversity conservation, and aesthetic landscape. The coefficients show substantial differences in unit-area ecosystem service values among land-use categories. Water bodies and wetlands generally have the highest values per unit area, followed by forestland, whereas grassland and cropland show comparatively lower values. The ESV coefficient for construction land was set to 0 because its capacity to provide ecosystem services is limited (Xie et al., 2015; 2017). This setting does not imply that construction-land change is ecologically neutral. Assigning a zero coefficient to construction land also excludes ecosystem services provided by urban green spaces embedded within built-up areas, such as parks and greenbelts, and may therefore slightly overestimate the ESV loss associated with conversion to construction land. Expansion of construction land can indirectly reduce regional ESV by replacing land-use categories that provide ecosystem services, including cropland, forestland, grassland, wetland, and water bodies. In addition, the coefficient for water supply from cropland is negative, reflecting water consumption associated with agricultural production (Xie et al., 2015; 2017). Table 3 presents the ecosystem service value coefficients used in this study, and ESV was calculated as follows:where ESV represents the total ecosystem service value of the study area, with a unit of 1 0 8 yuan; A i represents the area of land-use type i, with a unit of km²; V C i represents the ecosystem service value coefficient per unit area of land-use type i, with a unit of yuan·hm^-2·a⁻¹; and n represents the number of land-use types. Since the unit of A i is km², whereas the value coefficient is expressed in yuan·hm^-2·a⁻¹, the formula is multiplied by 100 to convert the area unit. The same ESV coefficients were applied to the 2020 baseline and all 2040 scenario projections, so differences among scenarios reflected changes in land-use composition and spatial allocation rather than changes in valuation coefficients over time. Accordingly, the valuation framework represents differences among land-use categories but does not differentiate ecological conditions or vegetation types within the same category.

E S V = 1 1 0 8 ∑ i = 1 n A i × 100 × V C i

TABLE 3. Ecosystem service value coefficients for the Shaanxi section of the Qinling Mountains (yuan·hm^-2·a⁻¹).

To examine the contribution of individual ecosystem services, the value of each service was calculated using the following formula:where E S V f denotes the value of the f-th ecosystem service function, with a unit of 1 0 8 yuan; V C i f denotes the value coefficient per unit area for the f-th ecosystem service function of land-use type i, expressed in yuan·hm^-2·a⁻¹.

E S V f = 1 1 0 8 ∑ i = 1 n A i × 100 × V C i f

To map the spatial distribution of ESV, the corresponding value coefficients were assigned to each land-use category in the 30 m × 30 m land-use raster (Nelson et al., 2009; Burkhard et al., 2012). For spatial analysis, the 30-m ESV raster was aggregated into 240 m × 240 m grid cells, each comprising up to 64 valid 30-m pixels. The ESV of each grid cell was calculated as the sum of the ESV contributions of its constituent pixels. The same grid system was subsequently used to calculate land-use conversion intensity and the spatial sensitivity index. The 240-m grid corresponds to an 8 × 8 aggregation of the original 30-m land-use pixels. This scale was selected to balance the retention of localized land-use transitions in narrow valleys, piedmont zones, and peri-urban areas with the need to reduce excessive sensitivity to isolated pixel-level changes. It also provides a common spatial unit for comparing ESV change and land-use conversion intensity across scenarios.

3.4 Spatial sensitivity of ESV responses to land-use change

We defined a grid-based sensitivity index (SI) as the ratio of relative ESV change to the proportion of pixels that changed land-use category within each 240 m × 240 m grid cell. The sign of SI indicates an ESV gain or loss, whereas its absolute value indicates response intensity. The index enables comparison among scenarios but does not measure ecological vulnerability or ecosystem resilience.

For grid cell g under scenario s, the relative change in ESV was calculated as:where R s , g E S V is the relative ESV change rate; E S V s , g 2040 is the projected ESV of grid cell g under scenario s in 2040; and E S V g 2020 is the baseline ESV of the same grid cell in 2020.

R s , g E S V = E S V s , g 2040 − E S V g 2020 E S V g 2020

The intensity of land-use conversion was calculated as:where R s , g L U is the proportion of valid 30-m pixels within grid cell g whose land-use category changed between 2020 and 2040; N s , g c h a n g e is the number of changed pixels under scenario s; and N g v a l i d is the total number of valid pixels in the grid cell. A complete interior grid cell contains 64 pixels, whereas cells along the study-area boundary may contain fewer valid pixels. A pixel was classified as changed if its land-use category differed between the 2020 baseline and the corresponding 2040 projection.

R s , g L U = N s , g c h a n g e N g v a l i d

The sensitivity index was then calculated as:

S I s , g = R s , g E S V R s , g L U

The sign and absolute magnitude of S I s , g have distinct interpretations. A positive SI indicates that land-use conversion is associated with an increase in grid-level ESV, whereas a negative SI indicates an ESV loss. An SI value of zero indicates that land-use conversion produces no net change in grid-level ESV. The absolute magnitude, ∣ S I s , g ∣, represents response intensity, with larger absolute values indicating stronger ESV responses to a given intensity of land-use conversion. Therefore, a negative SI with a large absolute magnitude represents a strong negative response rather than low sensitivity.

SI was calculated only for grid cells in which land-use conversion occurred and the baseline ESV was nonzero. Grid cells with R s , g L U = 0 were classified as unchanged and excluded from the SI calculation. When land-use conversion occurred but E S V g 2020 = 0, the relative ESV change rate was undefined; these cells were therefore treated as undefined and excluded from the SI calculation and scenario-level statistical summaries. Unchanged and undefined cells were treated separately from valid SI values in the spatial analysis.

For cartographic comparison, valid SI values from all four scenarios were displayed using the same continuous diverging color scale centered at zero. Positive responses were represented in red and negative responses in blue, with greater color intensity indicating larger ∣ S I ∣values. The same display range was applied to all four scenario maps to ensure direct spatial comparison. Raw SI values were retained for numerical analysis and were not converted into categorical sensitivity classes.

To examine whether grid cells with very limited land-use conversion disproportionately influenced SI estimates, we conducted an additional threshold-based robustness test. The main analysis retained valid grid cells containing at least one changed 30-m pixel, whereas the robustness test retained only cells containing at least two changed pixels ((N^{change}_{s,g} \geq 2)) within each 240 m × 240 m grid cell. The SI formulation itself was unchanged; only the minimum land-use conversion threshold used for inclusion was increased. This threshold-based test evaluates the influence of very limited within-grid land-use conversion on SI estimates; it does not constitute a multi-grain sensitivity analysis of the 240-m aggregation scale. We compared scenario-level SI statistics under the two thresholds in terms of mean, standard deviation, minimum, and maximum values and recorded the number of grid cells satisfying each threshold. Detailed results are reported in Supplementary Table S1.

4 Results

4.1 Historical land-use dynamics during 2000–2020

Historical land-use dynamics during 2000–2020 were examined in terms of spatial distribution, area change, and transitions among land-use categories using the 2000, 2010, and 2020 datasets (Figures 3, 4; Table 4).

TABLE 4. Changes in land-use areas in the Shaanxi section of the Qinling Mountains during 2000–2020.

Forestland remained the dominant land-use category throughout the study period, forming extensive and spatially continuous patches across the middle- and high-elevation mountainous areas (Figure 3). Cropland was concentrated primarily in river valleys, piedmont zones, peri-urban areas, and other relatively flat locations. Grassland occurred mainly as scattered patches between forestland and cropland, whereas wetlands and water bodies occupied comparatively limited areas along rivers, reservoirs, and low-lying zones. Construction land was concentrated around towns and transportation corridors and expanded noticeably between 2000 and 2020.

Table 4 shows that forestland decreased slightly from 66,717.06 km² in 2000 to 66,574.20 km² in 2020, corresponding to a net loss of 142.86 km². Its share of the study area declined from 78.05% to 77.88%. Cropland experienced the largest net reduction, decreasing from 13,837.66 km² to 13,233.87 km², a loss of 603.79 km². Grassland decreased by 73.34 km² over the same period. In contrast, construction land increased from 346.13 km² to 1,050.10 km², representing a net gain of 703.97 km². Water bodies and wetlands increased by 102.41 and 13.84 km². The magnitude of land-use change differed between the two subperiods. During 2010–2020, construction land increased by 542.85 km², substantially more than during 2000–2010, while cropland decreased by 521.84 km². Thus, most of the construction-land increase and cropland loss observed during the full study period occurred after 2010. The land-use transition diagram shows that persistence of existing land-use categories dominated the transition pattern (Figure 4). Forestland accounted for the largest retained area. New construction land originated mainly from cropland, with smaller areas converted from forestland and grassland. Conversions among cropland, forestland, and grassland were also observed, while wetlands and water bodies experienced relatively limited gains and losses.

Historical land-use dynamics were characterized by a relatively stable forest-dominated regional land-use structure, continued cropland loss, and accelerated construction-land expansion after 2010. The accelerated construction-land expansion during 2010–2020 occurred within a broader period of intensified urbanization and increasingly explicit ecological development controls. At the national level, the National New-type Urbanization Plan (2014–2020) provided a major policy framework for continued urbanization and spatial development (CPC Central Committee and State Council, 2014), while the regulatory framework for ecological protection and development control in the Qinling Mountains was revised in 2017 and 2019 (Standing Committee of the Shaanxi Provincial People’s Congress, 2019). These overlapping development and conservation pressures provided a relevant policy context for the observed land-use changes, although the present analysis does not attribute these transitions to any single policy.

4.2 Multi-scenario land-use projections for 2040

4.2.1 Historical hindcasting performance

Model performance was evaluated using the observed 2010 land-use map as the initial condition and by comparing the simulated 2020 land-use map with the observed 2020 map. The validation yielded an overall accuracy of 92.03% and a Kappa coefficient of 0.780. These results indicate that the calibrated model reproduced the observed 2020 land-use pattern with reasonable agreement. The hindcasting exercise served as a historical consistency check rather than validation of the predictive accuracy of the 2040 scenario projections.

4.2.2 Projected land-use structure

Table 5 summarizes the projected areas of each land-use category under the four scenarios in 2040. Forestland remained the dominant category under all four scenarios, accounting for more than 77% of the study area. However, substantial differences emerged in construction-land expansion, cropland retention, and changes in ecological land.

TABLE 5. Land-use structure in 2020 and under four scenarios in 2040.

Under natural development, cropland and grassland decreased by 353.06 and 45.32 km² relative to the 2020 baseline, respectively, whereas forestland increased by 237.94 km². Construction land increased by 160.05 km², while wetlands and water bodies increased slightly by 0.17 and 0.22 km², respectively. The projected changes partly followed the historical trend, with continued construction-land expansion and cropland and grassland loss, but with a modest increase in forestland. Urban development produced the largest construction-land area, reaching 2,502.94 km² in 2040, an increase of 1,452.84 km² from the 2020 baseline. Cropland decreased by 993.87 km² to 12,240.00 km², representing the largest cropland loss among the four scenarios, while forestland declined by 432.24 km². Under cropland protection, cropland increased to 13,450.00 km², 216.13 km² above the 2020 baseline, representing the largest projected cropland area among the four scenarios. Construction land remained essentially unchanged at 1,050.10 km², while forestland declined by 524.20 km² to 66,050.00 km² and grassland increased by 308.07 km². Under ecological protection, forestland and grassland increased by 325.80 and 243.27 km² relative to 2020, respectively. Wetlands increased by 9.00 km², while water bodies increased slightly by 0.44 km². Construction land decreased by 2.10 km², representing the only scenario in which construction-land area declined relative to 2020. The ecological protection scenario produced the largest net increase in the combined area of forestland, grassland, wetlands, and water bodies.

The main differences among scenarios concerned construction-land change, cropland retention, and ecological-land conservation: urban development produced the largest construction-land expansion, cropland protection retained the largest cropland area while maintaining construction land near the 2020 level, and ecological protection produced the largest net gain in ecological land while slightly reducing construction land.

4.2.3 Projected spatial patterns

The calibrated Markov–PLUS configuration spatially allocated the scenario-specific land-use demand to generate the projected land-use maps for 2040 (Figure 5). Across all scenarios, forestland remained widely distributed throughout the mountainous areas, cropland remained concentrated in river valleys and areas of relatively gentle terrain, and construction land expanded mainly around existing urban patches.

Under the urban development scenario, projected construction land expanded most extensively outward from existing urban areas and was concentrated mainly in river valleys, piedmont zones, and transportation-accessible locations. Construction-land expansion frequently occurred along the margins of cropland and forestland. Under the cropland protection scenario, more cropland patches were retained in river valleys and gently sloping areas, while construction-land expansion was less extensive than under the natural and urban development scenarios. Under the ecological protection scenario, forestland and other ecological land formed the most spatially continuous pattern, while construction-land expansion was most restricted. Increases in forestland, grassland, wetlands, and water bodies occurred mainly around existing ecological-land patches and localized transition zones.

The natural development scenario showed an intermediate spatial pattern. Construction land continued to expand outward from existing urban patches, but the extent of expansion was smaller than under urban development. The overall distributions of forestland, cropland, and grassland remained broadly similar to the 2020 baseline.

The strongest spatial contrast occurred between urban development, which produced the most extensive construction-land expansion, and ecological protection, which retained the most spatially continuous ecological land.

4.3 Ecosystem service value responses to scenario projections

4.3.1 Total and land-use-specific ESV responses

Table 6 summarizes total and land-use-specific ecosystem service values for the 2020 baseline and the four 2040 scenario projections.

TABLE 6. Ecosystem service value assessment results (10⁸ yuan).

Total ESV in 2020 was 2060.87 × 10⁸ CNY. Relative to the 2020 baseline, total ESV increased by 0.42%–2069.58 × 10⁸ CNY under ecological protection and by 0.24%–2065.72 × 10⁸ CNY under natural development. In contrast, total ESV decreased by 0.58%–2048.97 × 10⁸ CNY under cropland protection and by 0.84%–2043.60 × 10⁸ CNY under urban development. Ecological protection therefore yielded the highest projected ESV, whereas urban development yielded the lowest. Forestland contributed the largest share of total ESV in the 2020 baseline and under all four scenario projections. Under natural development, the ESV of forestland increased from 1912.66 × 10⁸ to 1919.49 × 10⁸ CNY, offsetting decreases in the ESV of cropland and grassland. Under urban development, the ESV of forestland and cropland decreased to 1900.24 × 10⁸ and 60.52 × 10⁸ CNY, respectively, and total ESV reached its lowest level among the four scenarios. Under cropland protection, the ESV of cropland increased to 66.51 × 10⁸ CNY and the ESV of grassland increased to 28.97 × 10⁸ CNY, while the ESV of forestland decreased to 1897.60 × 10⁸ CNY; total ESV remained below the 2020 baseline. Under ecological protection, the ESV of forestland increased to 1922.02 × 10⁸ CNY, while the ESV of grassland and wetlands increased to 28.56 × 10⁸ and 1.82 × 10⁸ CNY, respectively.

Differences in total ESV among scenarios therefore primarily reflected projected changes in forestland and other ecological land. Construction-land expansion affected ESV indirectly by replacing land-use categories with nonzero ESV coefficients.

4.3.2 Spatial patterns of ESV and scenario-induced changes

Figure 6 presents the spatial distribution of ESV under the four scenarios, whereas Figure 7 shows the absolute ESV gains and losses between the 2020 baseline and each 2040 scenario. These absolute changes should be distinguished from the sensitivity index in Section 4.4, which standardizes ESV change by the intensity of land-use conversion.

Across all scenarios, high-ESV areas were concentrated in extensively forested mountainous areas, whereas relatively low-ESV areas occurred mainly in river valleys, piedmont zones, and peri-urban areas with greater proportions of cropland and construction land. The broad regional pattern remained similar across scenarios, but localized gains and losses differed substantially.

Under urban development, negative ESV changes were most evident in river valleys, piedmont zones, peri-urban areas, and forest margins where construction land replaced cropland, forestland, or grassland. Localized positive changes occurred where water bodies, wetlands, or other ecological land increased. Under cropland protection, negative changes were more widely dispersed and occurred mainly where forestland or other ecological land with higher ESV coefficients was converted to cropland or other land-use categories with lower ESV coefficients. Positive changes were comparatively limited. Under ecological protection, positive ESV changes were most spatially extensive and coincided primarily with increases in forestland, grassland, and wetlands, together with localized increases in other ecological land. Negative changes were relatively localized. Under natural development, positive and negative ESV changes coexisted. Positive changes occurred mainly where forestland and other ecological land increased, whereas negative changes were associated primarily with construction-land expansion and localized losses of cropland and grassland.

Thus, regional ESV totals masked substantial localized gains and losses. The ecological protection scenario showed the most extensive positive ESV changes, whereas the urban development scenario showed pronounced negative changes in areas experiencing land-use conversion.

4.4 Spatial sensitivity of ESV responses to land-use change

Figure 8 presents the spatial sensitivity index for valid grid cells under the four scenarios. Positive SI values indicate increases in ESV, negative values indicate losses, and the absolute magnitude of SI represents response intensity. Only grid cells with valid SI values are interpreted below. Grid cells without land-use conversion and those with undefined SI were displayed separately and excluded from interpretation.

Valid SI values showed pronounced spatial variation. Strong positive and negative responses occurred mainly in river valleys, piedmont transition zones, peri-urban areas, forest margins, and localized areas of ecological-land conversion.

Under the urban development scenario, positive and negative responses were spatially interspersed. Negative responses were prominent in river valleys, piedmont zones, peri-urban areas, and localized forest margins, particularly where cropland, forestland, or grassland was converted to construction land. Positive responses occurred locally where wetlands, water bodies, or other ecological land increased. Under the cropland protection scenario, negative SI values were spatially dominant and were concentrated mainly in areas where forestland was converted to land-use categories with lower ESV coefficients. Positive-response cells were less widespread. Under the ecological protection scenario, positive SI values were the most spatially extensive. Strong positive responses occurred mainly in piedmont transition zones, areas adjacent to river valleys, and locations where cropland or construction land was converted to forestland, grassland, wetlands, or water bodies. Negative-response cells were comparatively limited and localized. Under the natural development scenario, positive and negative responses coexisted. Positive responses occurred mainly where forestland and other ecological land increased, whereas negative responses were associated with construction-land expansion and localized losses of cropland and grassland.

Urban development showed the broadest negative-response pattern, ecological protection showed the broadest positive-response pattern, cropland protection was dominated by negative responses, and natural development showed a mixture of positive and negative responses.

The threshold-based robustness test showed that, when the minimum land-use conversion threshold was increased from one to two changed 30-m pixels, 47.49%–68.83% of the grid cells with detected land-use conversion were retained across the four scenarios (Supplementary Table S1). The sign of the mean SI remained unchanged under all scenarios. Mean SI changed from 0.3448 to 0.3025 under natural development, from 0.5859 to 0.4892 under urban development, from −1.2767 to −1.2582 under cropland protection, and from 1.1925 to 1.1058 under ecological protection. Standard deviations decreased under natural development, urban development, and cropland protection and increased only slightly under ecological protection. The range of extreme values also narrowed substantially under the first three scenarios, particularly under urban development, whereas the maximum value under ecological protection remained unchanged. These results indicate that the direction and overall interpretation of scenario-level SI remained stable under the stricter conversion threshold, although some extreme values were sensitive to grid cells with very limited land-use conversion.

A transition-based attribution further identified the land-use conversions underlying negative SI responses (Supplementary Table S2). Under natural development, gross negative ESV change within negative-SI grid cells was mainly associated with forestland-to-cropland conversion (55.69%), cropland-to-construction land conversion (21.95%), and forestland-to-grassland conversion (20.98%). Under urban development, forestland-to-construction land conversion was the dominant transition, accounting for 95.42% of gross negative ESV change, followed by cropland-to-construction land conversion (3.53%). Under cropland protection, forestland-to-grassland and forestland-to-cropland conversions accounted for 57.28% and 42.71%, respectively. Negative responses under ecological protection were highly localized, and the gross negative ESV change within these cells was associated with water-body-to-wetland conversion.

5 Discussion

5.1 Ecological trade-offs among alternative land-use pathways

The four scenarios revealed clear trade-offs among construction-land expansion, cropland retention, and ecological-land conservation. Relative to the 2020 baseline, urban development produced the largest construction-land increase of 1,452.84 km² and reduced total ecosystem service value (ESV) by 0.84%. Cropland protection increased cropland by 216.13 km² while maintaining construction land near the 2020 level, but reduced total ESV by 0.58%. Ecological protection increased total ESV by 0.42% and slightly reduced construction land by 2.10 km², whereas natural development increased total ESV by 0.24%, with continued cropland loss but a modest increase in forestland. No scenario simultaneously maximized construction-land area, cropland area, and aggregate ESV. These results describe trade-offs among land-use priorities and their aggregate ESV outcomes rather than formally quantified trade-offs among individual ecosystem services (Bennett et al., 2009; Deng et al., 2021; Jiang et al., 2023). Because construction land was assigned an ESV coefficient of zero and the analysis did not independently evaluate food security, economic output, housing, infrastructure, or employment, the scenario comparison reflects ecological value rather than overall social welfare.

The absolute magnitude of ESV should, however, be compared cautiously across studies because estimates depend strongly on study extent, land-use classification, valuation coefficients, and the ecosystem services included. The 2020 baseline ESV estimated here (2060.87 × 10⁸ CNY) is lower than the 968.33 billion CNY yr^-1 reported in an earlier valuation of vegetation ecosystem services across the broader Qinling–Daba Mountains, but the latter covered a larger geographical domain and adopted a different valuation scope. Therefore, the numerical difference should not be interpreted directly as a difference in ecological performance. A more meaningful point of consistency across studies is the dominant contribution of forest ecosystems and the concentration of high ecosystem-service supply in relatively intact mountainous areas.

Urban development resulted in the lowest total ESV and the broadest negative-response pattern. Negative responses were concentrated in river valleys, piedmont zones, peri-urban areas, and forest margins where construction land replaced cropland, forestland, or grassland. The transition-based attribution further showed that forestland-to-construction-land conversion accounted for 95.42% of the gross negative ESV change within negative-SI grid cells, while cropland-to-construction-land conversion contributed a further 3.53%. This indicates that the negative response under urban development was driven primarily by the replacement of relatively high-ESV forestland with construction land rather than by conversion extent alone. Because SI standardizes ESV change by conversion intensity, these locations are not simply those with the largest converted area; they are places where a given level of conversion was associated with a comparatively large ESV loss. Similar concentrations of ecological loss have been reported in rapidly urbanizing transition zones (Deng et al., 2021; Jiang et al., 2023; Zhao and Shao, 2023). The projected regional ESV decline of 0.84% under urban development is relatively modest compared with historical ESV losses reported for more intensively urbanized regions. For example, Dai et al. (2021) reported a 30.92% decline in ESV in Chengdu between 2003 and 2018, while Yang et al. (2021) estimated a 24.72% decline in the Pearl River Delta between 1990 and 2017. These percentages are not directly comparable with the present scenario projection because the studies differ in time span, valuation coefficients, land-use classification, and treatment of construction land. Nevertheless, the contrast suggests that the relatively small regional loss projected for the Qinling study area reflects the persistence of an extensive forest-dominated landscape and the spatial concentration of development pressure, rather than an absence of substantial local ecological effects. This distinction is consistent with the strong localized negative responses revealed by the SI analysis (Dai et al., 2021; Yang et al., 2021).

Cropland protection retained the largest cropland area, but total ESV remained 0.58% below the 2020 baseline. The attribution results indicate that this negative response was associated almost entirely with losses of forestland: forestland-to-grassland conversion accounted for 57.28% of gross negative ESV change within negative-SI grid cells, while forestland-to-cropland conversion accounted for 42.71%. Thus, the ecological cost of this pathway arose not simply from cropland expansion itself, but from the conversion of relatively high-ESV forestland to land-use categories with lower ESV coefficients. This result does not diminish the importance of food production, which was not quantified independently in this study. Rather, it indicates that the ecological outcome of cropland protection depends on the broader land reallocation induced by the protection constraints, particularly whether relatively high-ESV forestland is replaced by land-use categories with lower ESV coefficients. Retaining existing contiguous cropland is therefore likely to impose a lower ecological cost than compensating for cropland loss through the conversion of forestland, wetlands, or water bodies (Bennett et al., 2009; He et al., 2020).

Ecological protection produced the highest total ESV, a slight decrease in construction-land area, and the most extensive positive-response pattern. Negative SI responses were comparatively rare and localized. Within these negative-SI cells, the gross negative ESV change was associated with water-body-to-wetland conversion. This should not be interpreted as evidence of ecological degradation caused by wetland expansion; rather, it reflects the coefficient-based valuation framework used here, in which water bodies have a higher unit ESV coefficient than wetlands. The positive ESV gains were associated with the retention or expansion of forestland, grassland, wetlands, and water bodies, consistent with the regulating and supporting services provided by intact forest and aquatic ecosystems (Ouyang et al., 2016; Watson et al., 2018; Wang et al., 2024). This outcome partly reflects the demand constraints and conversion rules imposed in the ecological protection scenario; it should therefore be interpreted as the ecological consequence of a conservation-oriented pathway rather than as a complete socioeconomic optimum.

Natural development increased regional ESV by 0.24%, but localized losses remained evident around construction-land expansion and areas of cropland, grassland, or localized forestland loss. Unlike those under the urban development scenario, these negative responses were distributed across a wider range of transition types. Forestland-to-cropland conversion contributed 55.69% of the gross negative ESV change within negative-SI grid cells, followed by cropland-to-construction land (21.95%) and forestland-to-grassland conversion (20.98%). This broader distribution of contributions indicates that negative responses under natural development arose from several concurrent land-use transition pathways rather than from a single dominant conversion process. The threshold-based robustness test further supports these scenario-level interpretations: excluding grid cells with only one changed 30-m pixel did not alter the direction of mean SI under any scenario, although the reduced magnitude of several extreme values indicates that the tails of the SI distribution remain sensitive to very low land-use conversion intensity.

The preferred pathway depends on the management objective: urban development provides the largest construction-land area, cropland protection retains the most cropland, and ecological protection yields the highest aggregate ESV. A future optimization scenario could combine compact construction-land growth, retention of existing cropland, and protection of ecological land in locations with strong ESV responses. Because such a combined scenario was not simulated here, it remains a proposition for future testing rather than an identified optimum.

5.2 Spatial priorities for land-use management

The spatial results suggest three broad management priorities: maintaining the continuity of ecological land, limiting land-use conversion in development-transition zones, and retaining existing cropland in suitable terrain (Wang et al., 2024; Zhang et al., 2024). These priorities should not be interpreted as statutory zoning because ecological conservation redlines, permanent basic farmland, and urban development boundaries were not included in the analysis. Accordingly, the scenario results provide comparative evidence on alternative land-management and planning orientations rather than spatial prescriptions for existing statutory planning zones.

In continuous forest areas, wetlands, water bodies, and headwater zones, conversion to construction land should be limited where possible. Areas that maintained high ESV across all scenarios may warrant priority protection, whereas locations with strong positive responses under ecological protection may be suitable for targeted restoration. Limiting dispersed construction and maintaining vegetated buffers along forest edges and rivers could reduce further fragmentation (Ouyang et al., 2016; Watson et al., 2018; Wang et al., 2024). Because landscape connectivity was not quantified, these recommendations concern the spatial continuity of ecological land rather than measured ecological connectivity.

In river valleys, piedmont zones, peri-urban areas, transportation corridors, and forest margins, a key concern is ESV loss associated with construction-land expansion. Infill development and redevelopment within existing built-up areas could help reduce outward expansion. Locations with strong negative SI values may warrant closer evaluation before further land-use conversion, particularly near forest edges and riparian zones (Deng et al., 2021; Jiang et al., 2023). These locations represent strong ESV responses rather than independently verified ecological-vulnerability zones. This spatial pattern is broadly consistent with previous studies of the Qinling Mountains. Yu et al. (2023) found that ecosystem-service-sensitive areas were concentrated mainly at lower elevations, whereas areas with high ecosystem-service supply were located primarily in the mountainous hinterland. More recently, Ye et al. (2026) identified major river corridors and urban–rural transition zones in the northern Qinling foothills as areas where high ecosystem-service value coincides with elevated ecological risk, while high-value ecological areas were associated more strongly with forest- and grassland-dominated mountainous landscapes. These findings support the distinction observed here between relatively stable high-ESV mountain interiors and more responsive low-elevation transition zones. The present SI analysis further refines this pattern by identifying locations where a given intensity of land-use conversion is associated with a disproportionately strong ESV response, rather than identifying high- or low-ESV areas alone (Yu et al., 2023; Ye et al., 2026). Compared with previous studies that primarily mapped ecosystem-service supply, sensitivity, or ecological risk, the marginal contribution of this study lies in linking local ESV responses explicitly to the intensity of land-use conversion, thereby distinguishing areas of high ecosystem-service value from areas where relatively limited land-use change may produce disproportionately strong ecological responses.

In valleys and gently sloping areas with concentrated cropland, existing contiguous cropland should be preferentially retained, and its conversion to construction land should be reduced. Cropland compensation should, where possible, avoid the conversion of forestland, wetlands, and water bodies. Steep or fragmented cropland may be considered for restoration where this is consistent with statutory plans. Because cropland quality and productivity were not evaluated directly, permanent protection priorities should be informed by the present results together with independent data on soil quality, agricultural productivity, and farmland regulation (He et al., 2020; Wang et al., 2020).

The mapped responses support compact development, retention of existing cropland, and targeted protection or restoration of ecological land. These priorities should be assessed together with agricultural production, socioeconomic demand, and statutory planning constraints because ESV represents only one dimension of land-use performance.

5.3 Limitations and future research directions

The 2040 land-use maps should be interpreted as conditional scenario projections rather than deterministic forecasts. Model performance was evaluated through historical hindcasting using the 2010 land-use pattern to reproduce the observed 2020 pattern. Although the resulting overall accuracy (92.03%) and Kappa coefficient (0.780) indicate reasonable agreement with the historical reference map, this validation demonstrates only the model’s ability to reproduce an observed historical pattern under the calibrated configuration. It does not establish predictive accuracy for 2040. Moreover, Kappa is reported here as a conventional summary measure to facilitate comparison with previous land-use modelling studies, but it should not be interpreted as a comprehensive indicator of simulation validity because it compresses different sources of disagreement into a single statistic. The 2040 projections remain conditional on the historical transition matrix, selected predictor variables, scenario-specific land-use demand, and conversion rules. Future policy adjustments, infrastructure development, climate change, socioeconomic dynamics, and other drivers not represented explicitly in the model may therefore produce land-use trajectories that differ from those simulated here. Accordingly, the four scenarios are intended to support comparison of alternative land-use pathways rather than provide deterministic forecasts of future land-use patterns.

The ESV assessment is also subject to uncertainty arising from the use of land-use-specific static value coefficients. The equivalent-factor approach assigns a uniform coefficient to all pixels within the same land-use category and therefore does not represent intra-class differences in ecosystem condition, vegetation type, stand structure, management intensity, or vertical zonation. This limitation is particularly relevant in the Qinling Mountains, where forestland spans pronounced elevational and climatic gradients and may include substantially different forest communities. As a result, the use of a single forestland coefficient may smooth spatial heterogeneity: ecosystem services may be overestimated in degraded or relatively low-function forest areas and underestimated in intact or relatively high-function forest areas. It may also underestimate changes in ecosystem-service provision that occur within the same land-use category, because changes in forest condition or composition are not represented unless the land-use class itself changes. The direction of bias in the regional total cannot be determined a priori because it depends on the relative extent and condition of different forest types. Similarly, assigning a zero ESV coefficient to construction land excludes ecosystem services supplied by urban green spaces and may therefore overestimate ESV losses associated with urban expansion. Accordingly, the ESV estimates should be interpreted primarily as standardized comparisons of ecosystem-service values among scenarios rather than as precise measurements of ecosystem-service provision at individual locations. Future studies could incorporate vegetation type, forest condition, elevation-dependent coefficients, remote-sensing indicators, or process-based ecosystem-service models to better represent intra-class heterogeneity.

The PLUS simulation is also constrained by predictor selection. The eight variables represent major terrain, accessibility, hydrological, and socioeconomic conditions, but do not explicitly include ecological restoration policies, climate-related changes, or tourism development. Policy priorities were partly represented through scenario-specific demand and conversion rules, but their spatial implementation was not directly modeled. This may lead to underestimation of land-use changes in areas where policy, tourism, or climate effects are strong. Future studies should incorporate spatially explicit policy, climate, and tourism variables.

The SI is a relative response measure, not an index of ecological vulnerability or resilience. Its value depends on the 240-m grid, baseline ESV, and the proportion of land-use conversion within each grid cell. Aggregation may smooth narrow or fragmented features, while cells with low baseline ESV or very limited land-use conversion may yield large ratios from modest absolute changes. To evaluate the latter effect, a threshold-based robustness test was conducted by increasing the minimum land-use conversion criterion from one to two changed 30-m pixels per 240-m grid cell. The direction of mean SI remained unchanged across all four scenarios, although several extreme values decreased in magnitude, indicating that the scenario-level interpretation was stable while the tails of the SI distribution remained sensitive to very low land-use conversion intensity (Supplementary Table S1). Cells with zero baseline ESV were excluded because relative change was undefined; this may omit gains where construction land was projected to convert to ecological land. Although the 240-m grid was selected to balance local spatial detail and the stability of grid-level comparison, the present robustness test addresses only sensitivity to very limited within-grid conversion and does not evaluate sensitivity to the spatial aggregation scale itself. Future analyses should therefore examine alternative grid sizes and grid origins. Comparing SI with an absolute metric, such as ESV change per unit area converted, would also help assess the influence of low baseline values.

Future work should therefore test the robustness of the projected patterns using expanded predictor sets, dynamic socioeconomic and climate scenarios, spatially explicit policy and tourism data, field observations, and more process-based ecosystem-service models.

6 Conclusion

Taking the Shaanxi section of the Qinling Mountains as a representative mountainous ecological barrier region, this study examined historical land-use dynamics during 2000–2020, projected land-use patterns for 2040 under four alternative scenarios, and assessed the associated ecosystem service value (ESV) responses and spatial sensitivity. The main conclusions are as follows:

• From 2000 to 2020, the regional land-use structure remained dominated by forestland, which decreased slightly by 142.86 km². Cropland experienced the largest net loss, decreasing by 603.79 km², whereas construction land increased by 703.97 km². Most of the construction-land expansion and cropland loss occurred during 2010–2020, indicating more intensive land-use conversion during the second decade. Persistence of existing land-use categories dominated the transition pattern, and new construction land originated mainly from cropland, with smaller areas converted from forestland and grassland. Historical land-use dynamics were therefore characterized by a relatively stable forest-dominated regional land-use structure accompanied by localized cropland loss and intensified construction-land expansion.
• The four 2040 scenarios retained a common forestland-dominated regional structure but differed substantially in construction-land expansion, cropland retention, and changes in ecological land. The urban development scenario resulted in the largest construction-land area and the greatest cropland loss. The cropland protection scenario retained the largest cropland area, although this was accompanied by reductions in some ecological-land categories. The ecological protection scenario generated the largest combined increase in forestland, grassland, wetlands, and water bodies, while slightly reducing construction-land area. The natural development scenario partly continued the historical pattern, with construction-land expansion and losses of cropland and grassland, but with a modest increase in forestland. Spatial differences among scenarios were concentrated mainly in river valleys, piedmont zones, and peri-urban areas: urban development produced the most extensive outward expansion of construction land, whereas ecological protection maintained the greatest spatial continuity of ecological land. These scenarios reveal clear trade-offs among accommodating development demand, retaining cropland, and conserving ecological land.
• Total ESV in 2020 was 2060.87 × 10⁸ CNY. Relative to this baseline, projected ESV increased by 0.42%–2069.58 × 10⁸ CNY under ecological protection and by 0.24%–2065.72 × 10⁸ CNY under natural development, but decreased by 0.58%–2048.97 × 10⁸ CNY under cropland protection and by 0.84%–2043.60 × 10⁸ CNY under urban development. Forestland contributed the largest share of ESV in the 2020 baseline and under all four scenarios. High-ESV areas were concentrated mainly in continuous forest-dominated mountain areas, whereas scenario-induced gains, losses, and strong ESV responses occurred primarily in river valleys, piedmont transition zones, peri-urban areas, forest margins, and localized areas of ecological-land conversion. The urban development scenario showed the most pronounced negative-response pattern, whereas ecological protection showed the most spatially extensive positive responses. Cropland protection showed localized negative responses where ecological land was converted to cropland, and natural development exhibited a mixed pattern of positive and negative responses.

Scenario rankings differed by management objective: urban development provided the largest construction-land area, cropland protection retained the most cropland, and ecological protection yielded the highest ESV. Regional totals nevertheless masked strong local responses in river-valley, piedmont, and peri-urban transition zones. These locations may warrant explicit consideration in future land-use optimization, while recognizing that the results depend on the scenario assumptions, fixed ESV coefficients, and the 240-m analytical grid.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fenvs.2026.1942894/full#supplementary-material

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Footnotes

Edited by: Merja H. Tölle, University of Kassel, Germany

Reviewed by: Yangyang Zhang, Chinese Academy of Sciences (CAS), China

Lei Wu, Shanxi Normal University, China