Water yield is estimated using Budyko’s water–heat coupling balance principle (Lin et al., 2021) through the water production module of the InVEST model. Higher water production indicates a greater water supply. The specific calculation formula is as follows:

In the formula, Y_(x) represents the annual water yield (mm) of a land use type; AET_(x) is the annual actual evapotranspiration (mm) of the grid unit; and P_(x) is the annual precipitation (mm) of the grid unit.

Habitat quality plays an important role in maintaining biodiversity. Based on the InVEST model, the main parameters were set with reference to previous research (Zhou et al., 2024), expert interview results, and the model’s user manuals, while also considering the characteristics of the study area. The calculation principle is as follows:

In the formula, Q_xj represents the habitat quality index of grid x in land use type j; H_j is the habitat suitability of land use type j; D_xj is the habitat degradation degree of grid x in land use type j; k is a semi saturation constant; and Z is the default model parameter.

Soil conservation services are evaluated using the revised universal soil loss equation (RUSLE) (Xiang et al., 2025). Greater soil erosion corresponds to higher sediment output and lower soil conservation. The specific calculation is performed using the following formula:

In the formula, A_c is the amount of soil conservation (t/hm²), determined by the difference between potential erosion (A_p) and actual erosion (A_r). R is the rainfall erosivity factor (MJ mm hm⁻² h⁻¹); K is the soil erodibility factor (t hm² h hm⁻² MJ⁻¹ mm⁻¹); LS is the terrain factor; C is the vegetation coverage factor; P is the factor of soil and water yield measures. In the model, the rainfall erosivity factor R is calculated from annual rainfall, the soil erodibility factor K is calculated based on the content of sand, silt, clay, and organic carbon in the soil, and C is calculated using the normalized vegetation index.

Carbon sequestration services are evaluated based on the four major carbon pools in the InVEST model. Carbon stocks are calculated by multiplying the carbon density of each land use type (Wang et al., 2018; Meng and Wu, 2023). The basic principle is as follows:

In the formula, C_above, C_below, C_dead, C_soil, and C_z represent aboveground carbon storage, underground carbon storage, dead organic matter carbon storage, soil carbon storage, and total carbon storage, respectively.

Marxan is a system protection planning model based on the simulated annealing method. It delineates the scope of protected areas under certain cost conditions (Ban et al., 2009) and helps in constructing, designing, and evaluating the spatial planning of protected areas (including land, marine, and freshwater systems). In this study, a 5-km² unit was used as the research unit, and the protection targets were set as 70%, 50%, and 30% of the study area, respectively (Woodley et al., 2019). The number of software iterations is set to 1,000,000, and the boundary length was corrected while keeping the remaining parameters unchanged. The model objective function is:

In the formula, T denotes the value of the objective function; Cost represents the total cost of the planning unit (PU); and Boundary is the length of the boundary of the protected area system. The boundary length correction coefficient (BLM) is a parameter that determines the aggregation of the protected area system. Penalty refers to the penalty incurred for failing to achieve the protection goal, calculated based on the protection cost of the planning unit. The protection feature penalty factor (CFPF) is used to emphasize the relative importance of different protection features.

The average habitat quality index of the Huaihe River Basin over the past 33 years was 0.399, 0.370, 0.363, 0.396, 0.383, 0.351, 0.349, and 0.349, respectively, and the habitat quality showed a fluctuating downward trend (Table 2). Habitat quality showed a downward trend in 1990–2000 and 2005–2023, with the most significant decline in 2010–2015, which may be related to human activities. The upward trend in 2000–2005 may be related to the Interim Regulations on Water Pollution Prevention and Control in the Huaihe River Basin promulgated in 1995, which reduced wastewater discharge and improved watershed habitat. From 1990 to 2023, the habitat quality was mainly at a relatively low level. The land types in 1990 (70.2%), 1995 (69.6%), 2000 (69.5%), 2005 (69.0%), 2010 (68.3%), 2015 (67.6%), 2020 (67.1%), and 2023 (66.9%) were all cropland.

On the whole, the habitat quality of the Huaihe River Basin is low, and there are patterns in its spatial distribution. From 1990 to 2023, areas with low habitat quality were scattered and expanded over time. The relatively low areas are mainly distributed in the plains of the Huaihe River Basin; the proportion of medium and relatively high areas is low, and the coverage change is not obvious, mainly concentrated in the central part of the Shandong Plain, the southern part of the Jiaodong Hills, and the southern part of the Dabie Mountains. The high areas are distributed in the water bodies and the mountainous regions at the edge of the Huaihe River Basin (Figure 3).

This change shows that the habitat quality is highly consistent with the geological characteristics of the basin. The land types corresponding to areas with high habitat quality in the Huaihe River Basin are grassland, forest, and water. These areas are at high altitude, steep in terrain, and less disturbed by human activities. Areas with relatively low habitat quality are mainly concentrated in cropland, mostly plains and gentle slopes, which are easily affected by human activities. Areas with low habitat quality are mainly concentrated in construction land, primarily due to highly concentrated urbanization and industrialization activities. Therefore, the disturbance of human activities is one of the main reasons for the decline in habitat quality.

From 1990 to 2023, the total carbon storage in the Huaihe River Basin showed a slight upward trend, increasing by 0.91 × 10⁸ t over 33 years (Table 3). Cropland carbon storage decreased from 32.7 6 × 10⁸ to 29.17 × 10⁸ t, but its contribution remained dominant. Construction land carbon storage increased from 2.76 × 10⁸ to 6.73 × 10⁸ t. This may be related to an increase in carbon stocks caused by policy interventions, such as ecological restoration projects including returning farmland to forests, returning grasslands, and natural forest protection. However, due to the expansion of construction land from 1990 to 2015, land types with higher carbon density (cropland, forest) are reduced, offsetting the substantial increase in total carbon.

From the spatial distribution of carbon storage in the Huaihe River Basin, the high carbon density values are mainly concentrated in the eastern part of Yantai, the eastern part of Wuhu, Rizhao, the northern part of Qingdao, and the southern part of the Dabie Mountains, while low values are concentrated near water bodies and construction land, showing a distribution pattern of high in the north and south and low in the middle (Figure 4). Based on the analysis of the distribution characteristics of land use types, the ecological land, such as forest and cropland, is consistent with the spatial layout of high-value areas, and the construction land and water are consistent with the low-value layout. Therefore, the total amount of carbon storage is closely related to land use type, and changes in land use have a significant impact on carbon storage.

According to the time change trends of water yield and rainfall in the Huaihe River Basin from 1990 to 2023, the annual precipitation was 810.70, 615.00, 727.03, 800.07, 689.29, 727.19, 859.94, and 659.27 mm, respectively. The annual average water depth was 265.77, 234.86, 208.32, 265.15, 219.20, 216.89, 266.32, and 204.03 mm, respectively. Both water yield and precipitation in the basin showed a “W”-shaped fluctuating downward trend (Figure 5). The water yield decreased by 57.45, 48.26, and 62.29 mm in 1990–2000, 2005–2015, and 2020–2023, respectively. In 2000–2005 and 2015–2020, it showed an upward trend, with increases of 56.83 and 49.43 mm. Therefore, the precipitation in 2020 was the highest, and the water yield reached its peak. In 2023, precipitation is low, and the average water yield falls into a trough.

The spatial distribution of water yield services in the Huaihe River Basin closely corresponds to precipitation, showing a generally positive correlation, and there is no obvious spatial correlation with actual evapotranspiration, indicating that the model simulates water production with high accuracy (Figure 6). High water yield values are mainly concentrated in areas of cropland, grassland, and construction land in the southeastern part of the basin, where the precipitation is higher, and vegetation coverage is dense. Strong transpiration and water vapor accumulation in these areas promote precipitation formation. Low water yield values are concentrated in the central and northeastern waters, where precipitation is lower, and transpiration is vigorous.

Soil conservation, from the perspective of temporal distribution, showed a pattern of decrease–increase–decrease over the past 33 years, with an overall downward trend. During this period, it decreased from 37.23 to 30.15 t, a reduction of 7.09 t or 19.0%, reflecting a continuous decline in soil conservation function. The range of soil conservation was 0–12,368.9 t, with an average value of 31.34 t. The area of high-value soil conservation increased, while the other value areas changed slowly. From a spatial distribution perspective, soil conservation exhibited relatively large spatial differences, generally showing the pattern of “high in the south and low in the north”. Most high-value areas are concentrated in grassland and higher-altitude regions (Figure 7), which are less affected by human activities. These areas have higher vegetation coverage and lower actual soil erosion, resulting in better soil conservation.

ArcGIS was used to assign the results of ecosystem services in 1990 and 2023 to the corresponding vector layers in county-level units, and local Moran index analysis was performed. High–low agglomeration and low–high agglomeration showed a trade-off relationship, and high–high agglomeration and low–low agglomeration showed a synergistic relationship. This was done to understand the spatial trade-offs/synergies of the service functions in the basin (Figure 8).

The high–high clustering area is concentrated in the mountainous hills in the south and west of the Huaihe River Basin. The terrain in this area is undulating, the precipitation is prone to cause surface runoff, and the coverage of forest and grassland is high, which is conducive to vegetation interception and soil water storage. At the same time, forest and grassland provide a superior environment for biology, vegetation photosynthesis is enhanced, and soil organic matter is continuously accumulated, so that the ecosystem service level in the south and west is higher than that in other areas, forming a “high–high” cluster; the low–low clustering is concentrated in the plain areas of cropland and construction land distribution in the eastern and western regions, mainly cropland. The cropland has been plowed for a long time, and the vegetation coverage changes with the seasons. The soil is exposed for a long time, and its soil retention capacity becomes weak. The habitat quality is affected by pesticides and fertilizers. The cropland reduces the soil water storage capacity, and many types of ecological service functions are synergistically degraded in the region, so that the ecosystem service function of itself and the surrounding area is at a low level, forming “low–low clustering”. High–low clustering sporadically appears in the transition zone from mountain to plain or plain forest and grassland coverage areas. The area has a high vegetation coverage, and the ecosystem service level, such as carbon storage and water yield, is high. However, due to being surrounded by cropland and construction land, it is easy to be disturbed by surrounding human beings, forming a high–low clustering with high itself and low around. Low–high clustering appears in the transition zone from plain to mountainous area, and the utilization type is mostly cropland, which makes the ecosystem service functions such as carbon storage, water yield, and soil conservation low, which may be affected by the radiation effect of the surrounding high–high clustering area. However, due to the limitation of land use type, it is difficult to effectively improve, forming a low–high clustering with low itself and high surrounding. Between 1990 and 2023, due to the expansion of construction land and the polarization phenomenon, the clustering of soil conservation, habitat quality, and water yield showed a decreasing trend, while the protection of returning farmland to forest and natural forest increased the clustering of high carbon storage and decreased other clustering.

Spearman correlation analysis identified a total of six correlations between the four ecosystem services and explored the correlation between the ecological functions in five different periods (Figure 9). Each ecosystem promotes the other, and the relationship between them is dominated by synergy. Among them, the relationship between habitat quality and soil conservation was strongly synergistic, and the average degree of synergy was above 0.45. The relationships between carbon storage and habitat quality, carbon storage and soil conservation, and carbon storage and water yield were moderately synergistic, and the average degree of synergy was between 0.21 and 0.4. The relationship between habitat quality-water yield and soil conservation-water yield showed a low synergistic relationship, and the average degree of synergy was between 0 and 0.2. Habitat quality-soil conservation showed a strong synergy. The lush vegetation in the mountainous and hilly areas of the basin provided habitats for animals and plants, reduced soil erosion, and created a positive feedback loop, enhancing both habitat and soil retention. Carbon storage and other ecosystem service functions maintain a moderate synergy, indicating that the other three ecosystem services are strongly affected by other related factors. The low synergy degree of habitat quality-water yield and even the trade-off in 2000 was mainly due to the fact that the higher habitat areas were mainly in mountainous and hilly areas, while the water yield capacity was not as good as that in plain areas, which made it difficult to achieve synergy between the two. In order to reduce soil erosion, the construction of terraces, engineering soil reinforcement, and other measures may reduce surface runoff to a certain extent, making the soil conservation-water yield synergy low.

From 1990 to 2023, among the ecosystems in the Huaihe River Basin, carbon storage-soil conservation, carbon storage-water yield, and habitat quality-soil conservation showed a downward trend, which was related to urbanization development and transitional development. Habitat quality-water yield, and soil conservation-water yield are on the rise. Policy support has improved water quality in the basin. At the same time, the development of projects such as returning polders to lakes and soil and water conservation is conducive to the simultaneous improvement of habitat and soil conservation.

Irreplaceability provides a protection priority sequence for all planning units, which facilitates assigning protection levels to each unit (Figure 10). Usually, areas with irreplaceability values above 80 are designated as first-level protection priority areas, which have the highest protection value. Areas with values in the range (60,80] are considered second-level protection priority areas, which high protection value. Areas with values irreplaceable in the range (40,60] are designated as third-level protection priority areas (Guo et al., 2018; Shi et al., 2018; Crouzeilles et al., 2015), which have moderate protection value and are overlaid with wetland protection areas, nature reserves, and soil erosion prevention and control zones. Among the 30% of the targets, only the third-level priority protection areas are located in water bodies, mountainous, and hilly areas, accounting for a small proportion of the total area. The 50% target secondary protected areas cover 2.09% of the study area and are situated near designated natural protected areas. This may be because nature reserves provide a stable ecological environment, and establishing priority protected areas nearby helps construct ecological corridors. Proximity to nature reserves facilitates resource sharing and reduces management costs. Overall, third-level protection areas account for 90% of the total area. Its essence is to provide light protection for the entire area while strengthening protection in key regions. The first-level priority area at 70% of the target is largely consistent with the second-level priority area at 50%. The protection priority of the core ecological area is highly stable and continuous, remaining the key protection under varying levels of protection effort. The second-level priority area has expanded significantly, covering more than 80% of the total study area, indicating that the areas with certain ecological value, though noncore, are included in the secondary protection level. The third-level protection area overlaps with soil erosion remediation zones, creating a spatial synergy between ecological protection and restoration. While maintaining ecological functions, each protection area should also undertake ecological restoration tasks, such as soil erosion control, thereby linking protection with restoration.

According to the division of different protection targets, it can be seen that 50% is the optimal result when the protection target is 50%. The reasons are as follows: covering more than 80% of the core ecological protection areas, the ecosystem services are better; since the entire study area is protected, it is easy to form the ecological radiation effect; the overlap with the existing nature reserve domain is high, so the maximum protection area is generated when the cost is the lowest.

Based on Zonal Statistics as Table software, the population density, per capita GDP, construction land area, and human footprint index were counted as cost factors. The evaluation results of ecosystem service function were standardized, and the areas with functional values greater than the average were analyzed as protection objects, so that the minimum cost could be used to achieve maximum biodiversity protection (Figure 11). The optimal protection area with a protection target of 30% was mainly concentrated in mountainous hills and high-biodiversity areas. When the protection target was low, the core area of ecological protection was preferred, and the maximum ecological protection benefit could be obtained with the minimum economic cost. The optimal protection area of 50%–70% protection target continued to expand to plains, waters, and wetlands, shifting from scattered distribution to contiguous development and reflecting the expansion from the core area to sub ecological value areas.