Exploring Trade-Offs and Synergies in Social–Ecological System Services across Ecological Engineering Impact Regions: Insights from South China Karst

Abstract

Karst ecosystems have become complex social–ecological systems (SESs) as a result of the interventions of large-scale ecological restoration programs, and the ecosystem services (ESs) that provide regional well-being can, to some extent, be described as social–ecological system services (S–ESs). Understanding the relationships among multiple S–ESs and exploring their drivers are essential for effective ecological management in karst areas, especially in regions differently affected by ecological engineering programs. Taking South China Karst (SCK) as a study area, we first identified two regions as comparative boundaries, namely significant engineering impact regions (SEERs) and non-significant ecological engineering impact regions (NEERs). Then we used ES assessment models, Spearman correlation, and optimal parameter geographic detector to identify the supply capacity, trade-offs/synergies, and their drivers of six types of S–ESs in SEERs and NEERs. The findings included: (1) SEERs were predominantly concentrated in the central and southern SCK regions, accounting for 33.98% of the total SCK area, with the most concentrated distribution observed in Guizhou and Guangxi. (2) Within the entire SCK, six S–ESs maintained a relatively stable spatial distribution pattern over time, with the most pronounced increase in soil conservation and a slight decrease in water retention, and the S–ES hotspots were more concentrated within the SEERs. (3) Most S–ES pairs within SEERs were optimized synergistically, with lower trade-off intensity and higher synergy intensity compared to NEERs. (4) S–ES pairs were affected by the interactions between the natural and socio-economic factors, with land use changes playing a crucial role, and natural factors were difficult to predict but cannot be ignored. Based on the results, we propose different SES sustainable development suggestions, with a view to providing theoretical support for the optimization of SES functions and the consolidating of integrated ecological construction.

1. Introduction

Ecosystem services (ESs) encompass products and benefits derived from nature that are crucial for human survival and development [1,2]. They serve as a link between natural environmental systems and socio-economic systems [3,4], and are categorized by the Millennium Ecosystem Assessment (MEA) into provisioning, regulating, supporting, and cultural services [5]. Ecosystem services research is rooted in ecological processes but is inherently a social endeavor [6]. The ecosystems crucial for human survival and development have evolved into complex social–ecological systems, rather than purely natural ecosystems. Social–ecological systems (SESs) are complex, dynamic natural ecosystems and social systems that are interdependent and closely linked, and constitute a system of human–earth relations [7]. Based on the discussed concept of SESs, some current studies argue that the concept of ES puts the benefits that humans obtain from ecosystems at the center, and places ecosystems in a subordinate position to serve human societies, ignoring the objective fact that humans are originally part of natural ecosystems [8], and dismisses the interaction between societies and ecosystems [9]. Therefore, some studies have proposed such terms as “social–cultural ecosystem services” [10], “social–ecological services” [11], “social–ecological system service” [12], or “society–ecosystem service” [13], as they more accurately reflect the interrelationship between social systems and ecological systems.

Relevant environmental policies strive to optimize multiple ESs through effective management [14], whereas the intertwined ecological and spatial dynamics often prevent the simultaneous maximization of all services [15]. Diverse human preferences can determine the feedback and interactions among ESs, manifesting as spatial trade-offs or synergies: trade-offs arise when one or several ESs are increased at the cost of others; synergies arise when multiple ESs are increased at the same time [16,17]. These relationships can respond to the external or internal changes (drivers) within a system [18]. Investigating and analyzing the trade-offs and synergies among different ESs has become central to current ES research, offering insights into optimal and suboptimal ecological management decisions [19,20]. As the research progresses, the concept of ES trade-offs/synergies has evolved from a simple understanding of objective laws to a decision-making tool for social–ecological systems [19]. Numerous researchers have quantified the ES trade-offs and synergies using various methods and scales including spatial mapping [21,22], statistical analysis [23,24], scenario analysis [25,26], and model simulation [27]. Most of the studies have measured the ES trade-offs and synergies for the entire study area. However, this approach cannot entirely capture the nuances of these relationships on a microscopic scale. Increasingly, there has been a growing focus on examining the spatial heterogeneity of service relationships, which helps to identify the leading functions of different regions and effectively allocate the environmental resources [28]. The transformation of ES relationships has been generally affected by certain factors such as climate, vegetation type, terrain, geomorphological features, land-use modes and management measures, urbanization, and ecological engineering [29]. Studies on the driving mechanisms behind ES relationships are moving towards integrated land use management [30]. Several researchers have leveraged the changes in ES interactions to assess the overall impact of ecological engineering programs [31,32]. In conclusion, analyzing multiple non-linear relationships between ESs and identifying driving factors, characteristics, and response rates at different scales is crucial for implementing appropriate policies [33]. Although many studies have been conducted to demonstrate that ecological restoration improves ESs, there is a lack of systematic evaluation of the effectiveness of restoration actions in enhancing ESs [34]. Moreover, most of the studies on ES trade-offs/synergies have focused on the local scale, with fewer studies at the macro scale [35].

Karst is an important component of the world’s fragile ecological zones, and extensive human activities, such as mining, agricultural and grazing activities, as well as urbanization and deforestation, have led to severe ecosystem degradation in many karst areas around the world [36]. In recent years, the phenomenon of karst ecosystem degradation has been reported in Southwest China [37], the Mediterranean [38], and the Caribbean [39], with studies related to the degradation of karst landscapes, the pollution of underground water resources, and karst disasters, etc. [40]. However, studies on large-scale ecological management of karst ecosystem degradation have mainly focused on China. China has approximately 3.44 million km² of karst area, accounting for approximately 36% of its total land area and 15.6% of the global 2.2 × 10⁷ km² karst area [41,42]. South China Karst (SCK) is one of the largest contiguous karst areas on Earth, with an exposed area of 5.5 × 10⁵ km² [43,44]. Due to the human disturbance, most of these areas are currently characterized by secondary coppices and shrubs, and some have experienced karst desertification (KD), which is difficult to restore [45]. Accompanied by severe soil erosion, extensive exposure of basement rocks, and a sharp decline in soil productivity, KD is gradually encroaching on production, living, and ecological spaces, severely threatening and limiting regional sustainable development [43,46,47,48]. In response to various social-ecological issues, the Chinese government has actively promoted karst desertification control (KDC) at the national level since 2000 [49]. Over the past two decades, ecological restoration programs in the SCK have notably increased vegetation cover, especially in karst areas [50,51]. However, this can only be a preliminary step in the restoration process. Recent studies at local and regional scales have begun to focus on ecosystem functions and services [52].

Following the implementation of several ecological engineering programs, the improvements of ecosystem functions in SCK are evident. These include enhancements in soil conservation [53], soil quality [54], and carbon sequestration and oxygen release [55,56]. However, there are also certain adverse impacts, such as reduced water yield [57] and declines in habitat quality and biodiversity [58]. These changes have increased the complexity of the relationships between ESs within the SCK. Therefore, a quantitative analysis of these relationships and their social-ecological drivers is crucial for enhancing the capacity of karst landscapes to provide benefits to society over time [36]. Numerous studies have focused on these objectives, utilizing either the entire SCK or typical karst regions as study areas. They have examined the impacts of ecological engineering on ES trade-offs and synergies [59], identified the key factors influencing these trade-offs and synergies to guide effective ecological protection and restoration policies [60], and predicted scenarios for ES trade-offs and synergies to inform the priority restoration recommendations and land use development [61,62]. Further research integrating social and ecological factors is essential to elucidate the underlying mechanisms of ES relationships in the context of ecological restoration, particularly through comprehensive studies across the entire SCK.

Given the large-scale ecological restoration measures in SCK, we adopted the term “social–ecological system service (S–ES)” to highlight the importance of human intervention in this region. Based on the above background, we first identified the significant engineering impact regions (SEERs) and non-significant ecological engineering impact regions (NEERs) in SCK. We then quantitatively analyzed the supply capacity, trade-offs/synergies, and trade-off/synergistic driving mechanism for six S–ESs: food production (FP), carbon sequestration and oxygen release (CS&OR), water retention (WR), soil conservation (SC), habitat quality (HQ), and ecological recreation (ER). This study focused on four specific objectives: (1) identifying different ecological engineering impact regions; (2) illustrating the spatiotemporal change characteristics of the S–ESs from 2000 to 2020; (3) exploring the S–ES trade-offs/synergies relationships and social–ecological drivers within two ecological engineering impact regions; and (4) providing some suggestions for the sustainable management of SES in SCK. The overall concept and data processing flow were depicted in a comprehensive analytical framework (Figure 1), organized into four main sections corresponding to the study objectives.

2. Materials and Methods

2.1. Study Area

SCK is situated between 97°21′–117°19′ E and 20°13′–34°19′ N, spanning eight provinces (cities) in Southwest China (Figure 2). Centered on the Yunnan–Guizhou Plateau and southern hilly regions, it covers approximately 1.94 × 10⁶ km², with the karst formations occupying 27.36% of this area [59]. It experiences an average annual temperature above 15 °C, with the annual precipitation exceeding 1100 mm, occurring predominantly during warm periods. The region features a special geological background, with strong karst erosion, dense distribution of mountain hazards, and a small ecological capacity, making the ecological environment very fragile. At the same time, this region is home to 100 million people living at a population density of 1.5 times the national average, and human–land conflicts are extremely acute [63]. It is one of the vulnerable areas with the most serious regional environmental problems and socio-economic conflicts in China.

2.2. Data Sources

The data in our study include natural ecological aspects such as land use, digital elevation model (DEM), soil properties, precipitation, and temperature on the one hand; on the other hand, socio-economic data from statistical yearbooks, government websites, etc. were also incorporated. Details of the data used and their sources are given in Table S1 of the Supplementary Materials.

2.3. Identification of Ecological Engineering Impact Regions

Three conditions were simultaneously established to identify the ecological engineering impact areas: (1) in sub-region 1, the Normalized Difference Vegetation Index (NDVI) value demonstrated a significant improvement (verified by the F test), with the slope trend higher than the regional average; this step was processed by Slope trend detection; (2) in sub-region 2, there was no significant correlation between NDVI and climate change (annual average temperature and annual precipitation), so this step was processed by Partial correlation detection in Matlab (v. 2016a) software; and (3) sub-region 3 fell within the boundaries of ecological programs implemented in SCK since 2000, with the main data sources including the “Dataset of Eco-efficiency Assessment of Major Ecological Projects in China, 2000–2019” [64] and the research of Shao et al. [65]. Overlaying these three sub-regions could yield significant ecological engineering impact regions (SEERs), while non-significant ecological engineering impact regions (NEERs) consisted of the regions that did not fulfill these conditions. The above identification methods and results were validated in similar research by Chen et al. [31]. The SEERs in this study primarily refer to regions with significant vegetation cover improvement through human governance over the past 20 years. The details of the slope trend and partial correlation detection formulas are provided in the Supplementary Materials.

2.4. Social–Ecological System Services Assessment

Given the naturally fragile characteristics of the region and insights from previous studies in SCK [57,58,59], this study evaluated six types of S–ESs. Here, S–ES is simply a name proposed based on the context of the study; the connotation and categorization methodology were not substantially changed, and the classification in our study still refers to the MEA [5]. The assessment methods of provisioning services (food production), regulating services (carbon sequestration and oxygen release, and water retention), supporting services (soil conservation and habitat quality), and cultural services (ecological recreation) are described below.

2.4.1. Food Production (FP)

FP serves as a fundamental provisioning service that is crucial for ensuring the resident production and livelihood, thereby contributing significantly to food security. Previous studies have established a strong linear correlation between food provision and the NDVI [66]. Accordingly, grain production for agriculture was allocated to the cultivated land based on the NDVI values [67]. The calculation formula for FP is as follows:

$$G_x={\displaystyle \frac{{NDVI}_x}{{NDVI}_{sum}}}\times G_{sum}$$

(1)

where $G_x$ (t/ha) is the food production of grid x; $G_{sum}$ (t/ha) is the total output of food production (rice, wheat, corn, soybean, and potato crops) for each city; ${NDVI}_x$ is the NDVI value in the xth grid; and ${NDVI}_{sum}$ is the total value of the NDVI for cultivated land.

2.4.2. Carbon Sequestration and Oxygen Release (CS&OR)

CS&OR plays a crucial role in regional climate regulation. According to the reaction equation for the vegetation photosynthesis and respiration, the formation of 1 g of dry matter involves the release of 1.63 kg of CO₂ and production of 1.2 kg of O₂. The volume of CS&OR can be calculated from NPP (net primary productivity) based on this ratio [68,69]. The NPP data were derived from the MODIS17A3 datasets, which were processed through splicing, formatting, reprojection, and other operations (the details are provided in Table S1 of the Supplementary Materials). The equation for calculating CS&OR is:

$${CS\&OR}_x={NPP}_x\times 1.63+{NPP}_x\times 1.2$$

(2)

where ${CS\&OR}_x$ (g C/m²) is the carbon sequestration and oxygen release of grid x; and ${NPP}_x$ (g C/m²) is the NPP value of grid x.

2.4.3. Water Retention (WR)

WR is a critical factor for the climate regulation in karst regions, as determined by water yield (WY) and a corresponding correction coefficient [70]. According to the InVEST model (v. 3.14.0), the WY adheres to the principles of the water cycle, incorporating precipitation, plant transpiration, surface evaporation, root depth, and soil depth. The adjustments included the topographic index, soil saturation water conductivity, and flow velocity coefficient. The details regarding the setup and calculations are available in the Supplementary Materials:

$${WR}_x=\mathrm{m}\mathrm{i}\mathrm{n}(1,{\displaystyle \frac{249}{Velocity}})\times \mathrm{m}\mathrm{i}\mathrm{n}(1,{\displaystyle \frac{0.9\times TI}{3}})\times \mathrm{m}\mathrm{i}\mathrm{n}(1,{\displaystyle \frac{Ksat}{300}})\times W{Y}_x$$

(3)

$$W{Y}_x=(1-{\displaystyle \frac{{AET}_x}{P_x}})\times P_x$$

(4)

where ${WR}_x$ (mm) is the water retention of grid x; $Velocity$ is the velocity coefficient without dimension; $TI$ is topographic index without dimension; $Ksat$ (cm/d) is the soil saturated water conductivity [71]; WY_x (mm) is the annual water yield of grid x; ${AET}_x$ (mm) is the average annual evapotranspiration of grid x; and $P_x$ (mm) is the average annual precipitation of grid grid x. The details regarding the setup and calculations are available in the Supplementary Materials.

2.4.4. Soil Conservation (SC)

Soil degradation is a prominent issue associated with KD. Therefore, the selection of SC as a representative supporting service is crucial for karst areas. SC was calculated by the RULSE model, known for its simplified structure, appropriate parameters, and effective simulation capabilities [72]. The equation for calculating SC is as follows:

$$SC=R\times K\times LS\times (1-C\times P)$$

(5)

where $SC$ (t/ha) represents soil conservation; $R$ represents the rainfall erosion factor [73]; $K$ represents the soil erosion factor [74,75]; $LS$ indicates the slope length slope factor [76,77]; $C$ represents the vegetation cover and crop management factor [78]; and $P$ represents the soil retention and water retention measure factor [79,80]. The details regarding the settings and calculations are available in the Supplementary Materials.

2.4.5. Habitat Quality (HQ)

Karst areas can face various ecological challenges including KD and biodiversity degradation. HQ is defined as the capacity to support biodiversity and can be assessed using the Habitat Quality module of InVEST, which can determine the HQ of land cover type j in each raster cell x within the study area. The habitat degradation index can quantify the stress level experienced by each cell grid. Higher values indicate greater stress, whereas lower values indicate better HQ on a scale from 0 to 1. Following the methodology proposed in references [81,82], the equation for calculating HQ is as follows:

$$Q_{xj}=H_j\left[1-\left({\displaystyle \frac{D_{xj}^z}{D_{xj}^z+k^z}}\right)\right]$$

(6)

where $Q_{xj}$ is the habitat quality of pixel x in landscape category j; $D_{xj}$ is the stress level of grid x in landscape category j; k is the half-saturation constant; z is the normalized constant; and $H_j$ is the habitat suitability of landscape category j. Details regarding the settings and calculations are available in the Supplementary Materials.

2.4.6. Ecological Recreation (ER)

Karst areas exhibit diverse landforms and offer a range of ecological landscapes and recreational opportunities. The evaluation of the ecological recreation involves utilizing an appropriate index framework based on the recreational potential and opportunity [83]. The principle of index universality and transferability was employed to select the most widely used indicators, which were then applied to evaluate the supply potential of the ER, considering the actual conditions of SCK. Six key indicators were adopted, including the degree of naturalness [84], landscape diversity, vegetation coverage, connectivity of natural attractions, elements of rivers and lakes, and accessibility level [85], and the data were processed as shown in Figure 3. The equation for calculating the ER is as follows:

ER_x = (R_1x + R_2x + R_3x + R_4x + R_5x + R_6x)/6

(7)

where R_1x is the degree of naturalness of grid x determined by the Hemeroby index; R_2x is the landscape diversity of grid x; R_3x is the vegetation coverage of grid x; R_4x is the connectivity of natural attractions of grid x; R_5x represents the elements of rivers and lakes of grid x; and R_6x is the accessibility level of grid x.

2.4.7. Hotspot Identification (Getis–Ord Gi*)

S–ES hotspots are areas where the provisioning capacity of one or more services is relatively outstanding within a given region, while cold spots are areas of weakness. With the help of the Getis–Ord Gi* module in ArcGIS (v. 10.8) software, we analyzed the S–ESs concentration of high value areas (hot spots) and low value areas (cold spots) at the county (district) scale, based on observing each S–ES distribution trend from 2000 to 2020. The equation is as follows:

$$G^{*}={\displaystyle \frac{\sum _{i=1}^n{Q}_{ij}-\overline{a}\sum _{i=1}^n{Q}_{ij}}{\sqrt{\frac{\sum _{i=1}^n{a}_i^2}{n}}-a^{-2}\sqrt{\frac{\sum _{i=1}^n{Q}_{ij}^2-{\left(\sum _{i=1}^n{Q}_{ij}\right)}^2}{n-1}}}}$$

(8)

where G* is the agglomeration index of grid i; a_i is the attribute of grid i; Q_ij is the weight matrix; n is the total number of cells; ā is the average of the total S–ES of all pixels.

2.5. Trade-Off or Synergy Analysis

• We employed the Spearman’s non-parametric correlation analysis, a quantitative method to measure the variation and strength of these interactions, to assess trade-offs and synergies among S–ES pairs [86]. Using ArcGIS (v. 10.8) software, we established a 5 km × 5 km grid to encompass the six types of S–ESs and collected sample points with varying S–ES values. Subsequently, using Origin (v. 2001) software, we applied Spearman’s correlation analysis to determine the correlation coefficients:

  $$r={\displaystyle \frac{\sum i(x_{ij}-\overline{x})(y_{ij}-\overline{y})}{\sqrt{\sum {i(x_{ij}-\overline{x})}^2\sum {i(y_{ij}-\overline{y})}^2)}}}$$

  (9)

  where $r$ is the correlation coefficient and the value domain [−1, 1]. If r > 0, there is a synergistic relationship between the two S–ESs, and r < 0 indicates that there is a trade-off between the two S–ESs. Moreover, $x_{ij}$ and $y_{ij}$ represent the data values of different types of S–ESs.
• The root mean square deviation (RMSD) can quantify the intensity of dispersion of an individual ES from the standard deviation of the average ES and has been widely used for the quantitative analysis of trade-offs. At this point, the meaning of trade-offs is not limited to characterizing negative trade-off relationships, but can also effectively express the degree of imbalance in the rate of isotropic change among ES synergies [62,87]. The calculation formula is as follows:

  $$RMSD=\sqrt{{\displaystyle \frac{\sum _{i=1}^n({ES}_i-\overline{ES})2}{n-1}}}$$

  (10)

  where $RMSD$ is the balance strength of the S–ES; n is the number of samples; ${ES}_i$ is the value of the ith standardized S–ES; and $\overline{ES}$ is the expected value of all types of S–ESs.

To eliminate the dimensional effects, each type of S–ES was standardized using the following equation:

$$E_X={\displaystyle \frac{E_x-E_{min}}{E_{max}-E_{xmin}}}$$

(11)

where $E_X$ is the standardized value for each S–ES, with its value ranging from 0 to 1; $E_x$ is an observed S–ES value; and $E_{xmin}$ and $E_{max}$ are the minimum and maximum observed S–ES values, respectively.

2.6. Driving Factor Detection

• Driving factors: Identifying the driving mechanisms is crucial for assessing the likelihood of trade-offs or synergies between ESs [88]. DEM, slope, lithology, and landform type were selected as some of the natural variables to represent the unique topography and geomorphology in karst regions; annual average temperature, annual precipitation, potential evaporation, and NDVI were widely used as climate and vegetation factors; light density, population density, and Gross Domestic Product were widely used to estimate socio-economic activities; and land use change reflected changes in human activities (

  Table 1. Driving factors types and discretization methods (taking factor-discretization methods of FP-CS&OR pairs in SEERs as an example).

  Factor Types	Factors	Data Types	Discretization Methods
  Natural	DEM (De)	Continuous	Equal breaks
  	Slope (Sl)	Continuous	Quantile breaks
  	Landform type (Lf)	Classified	-
  	Lithology (Li)	Classified	-
  	Annual average temperature (Te)	Continuous	Standard deviation breaks
  	Annual precipitation (Pr)	Continuous	Standard deviation breaks
  	Potential evaporation (Pe)	Continuous	Natural breaks
  	NDVI (Nd)	Continuous	Natural breaks
  Socio-economic	Light density (Lt)	Continuous	Quantile breaks
  	Population density (Po)	Continuous	Quantile breaks
  	Gross Domestic Product (GDP) (Gd)	Continuous	Quantile breaks
  	Land use change (Lu)	Classified	-

  ).
• Optimal parameter geographic detector (OPGD): The traditional geographic detector requires the manual setting for discretizing the continuous data, which is susceptible to inaccurate discretization and subjective factors. The OPGD addresses this issue by calculating the q-values across various grading methods and intermittent numbers, ensuring robust results across different spatial scales. Furthermore, it can extract the geographic features and spatial explanatory variables to reveal underlying patterns [89]. The GD package in R (v. 4.3.1) software was utilized with methods such as equal breaks, natural breaks, quantile breaks, geometric breaks, and standard deviation breaks. The classification series ranging from five to ten classes was set, and the spatial scale yielding the highest q-value was selected as the geodetector analysis parameter. The data types and discretization methods for the different factors are also listed in Table 1.
• Single-factor Detection: The core of factor detection lies in determining whether an independent variable x affects a dependent variable y by examining whether their spatial distributions converge. The calculation formula is as follows:

  $$q=1-{\displaystyle \frac{{\sum }_{h=1}^L{N}_h{\sigma }_h^2}{N_{}{\sigma }^2}}$$

  (12)

  where q is the explanatory rate of the influencing factors, with a value ranging from 0 to 1; the larger the value of q, the stronger the explanatory rate; x and y variables are superimposed to form an L-layer in the y-direction, which is denoted by h = 1, 2, …, L; N_h and N are the sample sizes of the sub-region h and the whole region, respectively; and σ²_h and σ² are the dispersion variances of sub-region h and whole region y, respectively.
• Interaction Factor Detection: The interaction detector assesses the interaction effects of two overlapping control variables by evaluating the relative importance of their interactions. It examines five types of interactions, including non-linear weakening, univariate weakening, bivariate enhancement, independent enhancement, and non-linear enhancement [90].

3. Results

3.1. Characteristics of Ecological Engineering Impact Regions

The extracting process of SEERs and NEERs is shown in Figure 4. SEERs were predominantly concentrated in the central and southern SCK regions, with the most concentrated distribution observed in Guizhou and Guangxi, which was consistent with the findings of Wang et al. [91]. These regions coincided with areas of high concentration and intensity of ecological engineering projects. NEERs were concentrated in the Sichuan Basin, upland areas of western Sichuan, eastern Hunan and Hubei, and regions along the southeastern coast of Guangdong, where no ecological engineering projects were implemented. In addition, most of these regions were non-karst. SEERs accounted for 33.98% of the total SCK area, more than a third of the total area. Notably, the SEERs in our study do not encompass all of the regions presenting enhanced ecological restoration effects. Instead, they primarily denote the regions in which substantial improvements in vegetation cover have been achieved through human management over the past two decades.

3.2. Spatio-Temporal Dynamic Characteristics of S–ESs

During the study period, the mean values of FP (+8.06%), CS&OR (+12.48%), SC (+22.64%), and ER (+1.29%) exhibited the increasing trends, while WR exhibited a slight decrease (−4.7%) and HQ decreased marginally (−0.9%). Spatial distribution patterns of each S–ES were broadly similar in 2000 and 2020 (Figure 5). The FP hotspots were concentrated primarily in the non-karst areas, such as the Sichuan Basin and the eastern part of Hubei. The CS&OR values were lower in the north and higher in the south, with prominent hotspots observed in south Yunnan, western Guangxi and most parts of Guangdong. The SC hotspots were predominantly high along the edge of the Sichuan Basin and the border area between Guizhou and Guangxi. The WR distribution corresponded closely with the rainfall patterns (Figure S2 in the Supplementary Materials), presenting higher values in the east and lower values in the west. Spatially, low HQ values were notably clustered in the Sichuan Basin and eastern Hubei, presenting an inverse distribution pattern compared with areas with high FP values. The ER hotspots were concentrated in most parts of Guizhou and Guangxi, and central Guangdong. Overall, the FP and HQ showed opposite trends spatially, with the remaining S–ESs showing large differences in spatial distribution.

The standard values of S–ESs within the entire SCK, SEERs, and NEERs are provided in Figure 6, which shows that CS&OR, WR, SC, HQ, and ER in SEERs were generally higher than those in NEERs and SCK, with HQ and ER showing the most pronounced performance. This suggests that SEERs were more conducive to enhancing regulating, supporting, and cultural services, but provisioning services were not involved. Furthermore, regions with three or more S–ES spatial hotspots were superimposed as class I regions, while coldspots and regions with fewer than three S–ES hotspots were grouped as class II, and the remaining regions were not significant. Class II regions were distributed over a small area within the SCK, mainly in south-eastern Guizhou and north-western and north-eastern Guangxi, and the proportion of class II regions within the SEERs was higher than within NEERs (Figure 6b). Overall speaking, the supply capacity of provisioning services in SEERs was lower than in NEERs, but the S–ES hotspots were more concentrated within the SEERs.

3.3. Characteristics of S–ES Relationships

3.3.1. Change Trends in S–ES Trade-Offs/Synergies

Figure 7 illustrates the 15 groups of correlations across different scales. Blue ovals indicate positive correlations, while orange ovals indicate negative correlations. On the SCK scale, all S–ES relationships passed the significance tests. A trade-off relationship was observed between FP and the other five services. The strongest trade-off was between FP and HQ, with an average annual correlation coefficient of −0.466. All other services demonstrated varying degrees of synergy, with WR-HQ exhibiting the strongest synergistic interaction with an average annual correlation coefficient of approximately 0.416.

The FP-HQ trade-off was the most significant in both SEERs and NEERs, with average correlation coefficients of −0.419 and −0.336, respectively. The most synergistic relationship were different in SEERs and NEERs: CS&OR-WR was the most pronounced in SEERs, with an average correlation coefficient of 0.410, while WR-HQ was the most in NEERs (0.319). Some pairs performed extremely differently, even developing in an inverse direction. For FP-ER, we observed continuous synergistic trends in SEERs and NEERs compared to the SCK. The FP-CS&OR within the two boundaries gradually shifted from a trade-off to a synergistic relationship, increasing in significance. In SEERs, the HQ-ER exhibited a trade-off relationship initially but subsequently evolved into a synergistic relationship, exhibiting increased significance and synergy coefficients. The CS&OR-SC in NEERs consistently showed a trade-off relationship that intensified over time. From 2000 to 2020, most S–ES pairs within the SEERs boundaries exhibited an increasing trend towards synergy, and more than half of the S–ES pairs within the NEERs moved towards trade-offs.

In order to better articulate the changes in trade-offs/synergies between S–ESs, as well as to further explore the trade-off/synergy intensity, we quantitatively explored the S–ES pair relationships within different regions over the years 2000–2020 by taking the average S–ES values as study targets. Specific results are displayed in Figure 8. In addition to these pairs of FP-SC, WR-SC, SC-ER, and HQ-ER, there was a greater degree of synergy in S–ESs in SEERs than in NEERs.

3.3.2. Trade-Off/Synergy Intensity

In this study, the greater the RMSD value between S–ES pairs in a synergistic relationship was, the greater was the degree of synergy. The RMSD distribution of the S–ES pairs in SEERs and NEERs is illustrated in Figure 9. The highest synergy intensity emerged in WR-HQ, with RMSD values of 0.440 in SEERs and 0.389 in NEERs. The FP-HQ exhibited a notable degree of trade-off, with RMSD values of 0.334 in SEERs and 0.383 in NEERs. The RMSD in these two pairs were spatially similar, although their correlation coefficients were opposite, with high RMSD values in SEERs concentrated in Guizhou, Guangxi, and eastern Hunan, and in the NEERs concentrated in western Sichuan. The degree of synergy for the WR-SC was the lowest, with RMSD values of 0.041 for SEERs and 0.051 for NEERs. The differences in the trade-off or synergy intensity for the remaining S–ES pairs were minimal, with RMSD values ranging from 0.109 to 0.266 for SEERs and 0.120 to 0.271 for NEERs. In a word, most S–ES pairs within SEERs had lower trade-off intensity and higher synergy intensity compared to NEERs.

3.3.3. Identification of Social-Ecological Drivers Based on OPGD

• Single-factor detection

The q-values within SEERs and NEERs are shown in Figure 10. Generally, the S–ES pairs within the SEERs were more responsive to the natural and socio-economic factors. In the SEERs group, all factors were significantly correlated (F < 0.001). Within the NEERs, the night light index did not significantly affect the CS-ER relationship (F = 0.518) and was therefore excluded from further analyses.

Table S7 presents the top three q-statistic values for the impact of the S–ES pairs. Some S–ES pairs were affected by the same single factor within both SEERs and NEERs boundaries. The NDVI significantly contributed to the FP-CS&OR, with q-values of 0.424 in SEERs and 0.220 in NEERs. The land use change exerted the greatest explanatory power on the FP-WR, FP-SC, FP-HQ, WR-HQ, and SC-HQ relationships, with q-values of 0.439, 0.385, 0.098, 0.333, and 0.230 in SEERs, and 0.377, 0.265, 0.169, 0.378, and 0.214 in NEERs, respectively. The DEM was the dominant factor influencing the FP-ER, with q-values of 0.396 in SEERs and 0.250 in NEERs. The annual average temperature had the largest effect on the CS&OR-WR and SC-ER, with q-values of 0.436 and 0.432 in SEERs, and 0.350 and 0.186 in NEERs. The effect of the slope on the WR-SC was notable, with the q-value of 0.104 in SEERs and 0.132 in NEERs, respectively.

The remaining S–ES pairs were influenced by different factors. Within the SEERs, the land use change also exerted the greatest explanatory power on the CS&OR-HQ (q = 0.221) and HQ-ER (q = 0.190). GDP played a significant role in shaping the relationship of CS&OR-ER, with a q-value of 0.113. The annual average temperature significantly affected the CS&OR-SC (q = 0.439) and WR-ER (q = 0.366). Within the NEERs, the DEM played a dominant role in shaping the WR-ER, with a q-value of 0.132. The annual average temperature had notable effects on the CS&OR-HQ, with a q-value of 0.142. The CS&OR-SC was primarily influenced by annual precipitation (q = 0.413). The landform type significantly affected the CS&OR-ER, with a q-value of 0.183. The GDP exhibited the greatest explanatory power in the HQ-ER, with a q-value of 0.061.

2.

Interaction factor detection

The interaction analysis suggested no mutual weakening of the interaction combinations within either SEERs or NEERs boundaries. The interactions were predominantly identified as bi-enhanced and non-linear enhancements. Table S8 illustrates that the influence of the factors increased under interaction, indicating that no individual factor in the study area singularly governed the regional heterogeneity of each S–ES relationship.

Within the SEERs boundary (Figure 11), nearly half of the S–ES pairs were influenced by the interactions involving the land use changes and other natural factors, primarily of the bi-enhanced type. Specifically, the interaction between land use change and DEM significantly affected the FP-SC, FP-HQ, CS&OR-HQ, WR-HQ, and SC-HQ. Some natural factors presented significantly less influence when considered independently but gained significance when interacting with other factors. For instance, the interaction of annual average temperature with the DEM increased the q-value to 0.581 in the CS&OR-WR (compared with 0.436 for annual average temperature alone). Similarly, the interaction of DEM with NDVI in the FP-ER resulted in a q-value of 0.434 (compared with 0.396 for the DEM alone). The interaction between the annual average temperature and landform type had the greatest impact on the CS&OR-SC and SC-ER (q = 0.529, 0.477). Some natural factors such as slope, landform type, and lithology, while individually exhibiting lower impact values, enhanced their effects when interacting with each other. For instance, the interaction factor of landform type and lithology in the CS&OR-ER reached 0.184 (compared to 0.111 for landform type alone).

Within the NEERs boundary (Figure 12), the interaction of the land use change with other factors significantly affected most S–ES pairs. Particularly noteworthy was the interaction of NDVI with other factors, which exerted considerable impacts on the FP-CS&OR, FP-ER, and CS&OR-WR relationships, notably in the CS&OR-WR (NDVI × annual precipitation: 0.606). For CS&OR-SC and SC-ER, the dominant interaction factors in this area mirrored those in the SEERs, specifically the annual average temperature × landform type. Other natural factors, which were initially less impactful when considered individually, demonstrated a strong impact after interacting with other factors.

4. Discussion

4.1. Spatiotemporal Variations in S–ESs and the Necessity for Continuous Ecological Engineering Programs

4.1.1. S–ESs Supply Capacity

The consistent enhancement of CS&OR and SC can indicate the effectiveness of forest conservation and restoration efforts in karst areas [92,93]. However, favorable vegetation cover can limit the increase in WY under similar weather conditions because of its water-trapping and evapotranspiration capabilities [61]. Consequently, our study suggested a declining WR trend during the study period. The increase in FP can be attributable to the ongoing agricultural reforms and the adjustments in cultivation and farming practices. Statistical analysis indicated that the total grain production in the study area increased from 1.54 × 10⁷ t in 2000 to 1.59 × 10⁷ t in 2020. The HQ was closely tied to the land suitability for habitation and was threatened by human activities and landscape fragmentation. The built-up land expanded by 104.17% between 2000 and 2020 in SCK (Table S9 in the Supplementary Materials), indicating increased human disturbances and their impacts on the HQ. The opposite trend in the spatial distribution of FP and HQ hotspots also demonstrates the important impact of land fragmentation and human activities on HQ. The capacity of the ER to respond comprehensively to various indices has seen an observed increase, likely due to some key factors. For instance, vegetation restoration contributed to the enhanced landscape recreational potential, evident in the 7.3% increase in NDVI from 2000 to 2020. Moreover, urban developments improved the transportation infrastructure, thereby facilitating easier access to recreational opportunities.

During the study period, the implementation of ecological projects was more conducive to enhancing the supply of regulating, supporting, and cultural services, and S–ES hotspots areas had higher proportions in the SEERs than in NEERs. This revealed the effectiveness of ecological management in enhancing the S–ES provision capacity. Taking KDC as an example, by the end of 2015, the KDC had received over 1.3 × 10¹¹ CNY in total investment, with approximately 1.2 × 10¹⁰ CNY sourced from the central government, covering an area of approximately 7.0 × 10⁴ km² [94]. The implementation of various ecological restoration initiatives led to notable shifts in the land-use patterns, influencing the dynamics of S–ESs across the study area.

4.1.2. S–ES Trade-Offs/Synergies

Previous studies have generally assumed a trade-off relationship between the provisioning services and regulating or supporting services [95]. This pattern is evident across the entire SCK but does not hold uniformly within the SEERs and NEERs. The spatial heterogeneity in the ecosystem distribution can be attributed to the variations in ecological context and resource conditions at different geographic scales, thereby influencing the ecosystem relationships [96]. There was a general trend towards synergistic interactions between the S–ES in SEERs, and the synergy intensity in SEERs tends to be larger than that in NEERs. Studies have confirmed that the synergy among ESs can be typically greater in ecological engineering areas than in the non-engineering areas, with the afforestation programs affecting the vegetation cover changes and subsequent ES functions [91]. Therefore, continuous ecological engineering efforts can play a significant role in enhancing coordination and achieving ecological management goals.

A detailed comparison of the FP-SC and CS&OR-SC pairs within the SEERs and NEERs boundaries revealed significant differences. Within the SEERs, the FP-SC exhibited a trade-off trend, with the highest RMSD values concentrated in north-central Guizhou, western Hunan, and southeastern Hubei, characterized by the mountainous terrain and poor soil texture. Practices such as steep-slope plowing and sloppy cultivation can reduce the capacity of the soil to store rainfall, resulting in decreased SC and increased FP [97]. The trade-off trend between CS&OR and SC was larger in NEERs. Since the implementation of the ecological construction, the vegetation cover in karst areas has generally increased, with higher NPP and effective suppression of soil erosion, thereby increasing the SC; moreover, the NPP and SC can exhibit a synergistic relationship [98]. This finding confirmed the importance of ecological restoration projects in karst areas for regional SC improvement.

4.2. Single/Interaction Factor Characteristics and the Necessities of Zoning Restoration and Influencing Factor Management

The results of single-factor detection highlighted that natural factors and land use changes were pivotal in explaining the spatial patterns of S–ES trade-offs. The natural environmental factors, such as rainfall, vegetation coverage, and elevation, could form the foundational elements to affect the spatial trade-offs and synergies [35,99]. Land use change can be another critical factor directly affecting the S–ES trade-offs by altering the ecosystem composition and structure to affect the ES supply dynamics [100]. Chen et al. emphasized that the influence of various factors on the trade-off of ESs could exhibit significant spatial imbalance, with forestland having the most prominent impact in a karst region [31]. Furthermore, the study indicated that the S–ES relationships within the SEERs could be more responsive to both the natural and socio-economic factors than the NEERs under similar single-factor impacts. This sensitivity suggests that the environmental and anthropogenic changes within the SEERs could exert a more pronounced impact on the S–ES relationships.

The interaction detection indicated that the S–ES relationships could be affected by bi- or non-linear enhancement effects, which could be supported by research in karst areas. Previous studies have highlighted that the interactions between different drivers can enhance the explanatory power of the spatial variations in the trade-off strength [101,102]. For the S–ES pairs related to the FP, the intensity of trade-offs was notably affected by the interactions involving land use and other natural factors. Research has indicated that the impact of FP on various ESs can result largely from the changes in land-use patterns [103]. Regarding the S–ES pairs linked to the HQ, the interactions between the land use change and other factors such as DEM, annual average temperature, and population density can exhibit stronger effects than the interactions within the natural or socio-economic factors alone. Therefore, the results of the interactions between multiple factors should be considered in the SES management. The development and land use should be selected to align with natural resources, and the economic and social development. For certain S–ES pairs, such as CS&OR-WR, CS&OR-SC, WR-ER, and SC-ER, the interactions between meteorological factors (annual average temperature, potential evaporation, and annual precipitation) and topographic factors (landform types, lithology, DEM, and slope) exerted significant effects. This was particularly pronounced within the SEERs boundary, where these interactions were consistently synergistically enhanced. The karst areas were characterized by complex and spatially heterogeneous structures, where the landforms, climate conditions, vegetation, and soil properties varied significantly. These factors could affect the ESs and their spatial relationships, whereas macroscopically, the landform types could directly shape the natural ecological environment, affecting the supply and maintenance of ESs [104]. Given the substantial topographic relief within the SEERs, it could be crucial to differentiate and strengthen management strategies tailored to the different landform zones.

Notably, NDVI could be a common factor influencing the trade-offs/synergies in the FP-CS&OR and FP-ER within both the SEERs and NEERs boundaries, as well as in the CS&OR-WR and CS&OR-HQ, specifically within the NEERs. Research has demonstrated the significance of NDVI in these S–ES relationships, particularly in shaping the trade-off patterns between the FP and NPP [105]. The increase in NDVI driven by the ecological programs notably enhanced the NPP and HQ [32]. As mentioned earlier, NDVI serves as a crucial indicator for assessing the ER potential. These findings have indicated the importance of enhancing the NDVI to enhance the synergistic trend among FP, CS&OR, ER, and HQ. In the future, we should pay more attention to the effects of factor interactions on S–ES relationships when analyzing the influencing factors, and implement the factor effects as much as possible, so as to better serve ecological management and decision-making. We believe that it is difficult to control climate-type impacts of natural factors at the anthropogenic level. Natural factors such as topography and vegetation cover can usually be improved through land use.

4.3. Suggestions for SES Sustainable Management in SCK

(1)

It is necessary to consolidate the results of previous ecological engineering construction and seek more integrated protection-restoration projects. In some regions where ecological engineering is concentrated, some measures could involve large-scale transformation and the reconstruction of inefficient plantation forests, upgrading and transforming the shrub forests on the slopes and foothills of KD mountains, afforestation beneath the forest canopies, the selection of ecologically and economically viable forestry and grass species, and the development of derivative ecological industries, such as agroforestry and specialty economic crops and medicinal plants [94]. In terms of integrated protection–restoration aspects, the “mountains–rivers–forests–farmlands–lakes–grasslands” program provides a suitable development direction, and it highlights the ecosystem connectivity and diversity through vegetation restoration, effectively enhancing ecosystem stability and quality [106]. Furthermore, in smaller karst ecological engineering implementation areas (watersheds, administrative units, or grid units), it is necessary to manage zoning according to typical ES supplies and relationships, so as to better carry out ecological effect assessment and later monitoring. At the same time, the concept of ecological priority should be strongly advocated in those areas where the impact of ecological engineering is small or where there is no ecological engineering implementation.

(2)

It is necessary to strictly control land development planning at all levels of units, especially those ecological engineering regions that are more intensely subject to anthropogenic factors. On the one hand, the development strategy of ecological red line and arable land red line should be strictly enforced. For instance, as the foundation of the national and regional ecological security, the ecological protection red line preserves the essential functions of key ecosystems and enhances their ecological support capacity for social and economic development [107]. Preserving the red line of the cultivated land is fundamentally vital for ensuring the food security and enhancing the productivity and efficiency of agricultural spaces [105]; this is especially important in areas with high FP values. On the other hand, attention should be paid to the problem of ecological land compression brought about by the expansion of construction, and increasing green patch sites around existing built-up areas and introducing green space systems are necessary [108].

(3)

There is a need to promote the study of ES trade-offs/synergies under different scenarios in order to adjust SES development strategies at different scales. Scenario modeling is an effective way to integrate land use, climate change, and factor management to provide decisions for SES development. Future scenarios can not only predict the evolution of trade-offs/synergies on the supply side of ESs [62,109], but can also simulate the evolution of trade-offs/synergies on the supply–demand side of ESs [110]. This is a scientific approach and a favorable reference for optimal SES development at different scales, and decision makers can be helped to develop targeted and differentiated strategies [111].

(4)

There is a need to play market mechanisms and improve diversified ecological protection and restoration paths. Firstly, enterprises, private individuals and non-governmental organizations can be introduced to reduce the government’s financial and management burden. Secondly, county governments can give full play to their autonomy, vigorously broaden the channels for realizing the value of ecological products, and actively strive for the introduction of ecological industry-type projects, so as to transform rich natural resources and high-quality ecological environments into ecological products and enhance market competitiveness. Thirdly, certain policy preferences and financial incentives will be given to regions with good results in the implementation of ecological projects, so as to balance the ecological restoration needs of different regions.

4.4. Research Limitation and Outlooks

Our study had several limitations. First, although we aimed to simulate typical S–ESs in the study area as accurately as possible, uncertainties remained owing to the utilization of proxies and models. For example, the WR and SC were evaluated in our study, but due to the special “dual” structure above and below ground in the karst region, surface water and soil seepage occurs along rock cracks, funnels, and sinkholes, and migrate to other areas [112]. Therefore, the assessed values of these two types of S–ES may be higher than the actual values. Second, the karst and non-karst landscapes exhibited distinct behaviors concerning the ES trade-offs and synergies [113], suggesting the need for a flexible methodology to explore these relationships across different landscape types. Third, although we incorporated many influencing factors based on existing studies, some of the factors that have a significant impact on karst areas, such as mining, tourism, and urbanization, were not considered, and current research has not spatially explored positive or negative changes in natural and socio-economic factors. In the future, more scientific impact factors should be considered, and spatial characteristics of the positive or negative impacts on the S–ES relationships can be explored by spatial models (e.g., geographically weighted regression models), with a view to guiding development policies at different scales (administrative or grid units). Fourth, the effects of various influencing factors on S–ES relationships are not linear; thus, we should focus on exploring the thresholds at which various factors affect S–ES relationships to better predict the trends of trade-offs/synergistic transitions and make better decisions.

5. Conclusions

By overlaying three sub-regions, we identified two ecological engineering impact areas (SEERs and NEERs). We then assessed the key S–ESs (FP, CS&OR, WR, SC, HQ, and ER) in 2000, 2010, and 2020 with various assessment models. Further, we compared the S–ES trade-offs/synergies relationships, and their influencing drivers within SEERs and NEERs.

Our study revealed several key findings. First, SEERs were predominantly concentrated in the central and southern SCK regions, with the most concentrated distribution observed in Guizhou and Guangxi. Second, the mean values of the S–ESs within the SEERs, except for the FP, exceeded those in the NEERs and the entire SCK, and the S–ES hotspots were more concentrated within the SEERs. Third, the trade-off and synergistic relationships varied across different boundary scales, with most S–ES pairs within the SEERs exhibiting synergistic trends, with lower trade-off intensity and higher synergy intensity compared to NEERs. Finally, the two-factor interactions exerted a stronger impact on the S–ES trade-off intensity than the single-factor interactions, with the land use change playing a particularly significant role. Over the past two decades, human activities have emerged as an important driver affecting the changes in the S–ES relationships in SCK. Our study provides theoretical support for the optimization of SES functions and the consolidation of ecological construction in karst areas, even in vulnerable areas.