The Temporal Variation and Spatial Scale Dependence of the Trade-Offs and Synergies among Multiple Ecosystem Services in the World Heritage Site of South China Karst

Abstract

The World Natural Heritage Sites (WNHSs), which have unique ecosystems, ecological landscapes, and biodiversity, are the common heritage of all of humanity. The evolutionary pathway of ecosystem services (ESs) in the Karst WNHS between the years 2000 and 2020 has been examined, and the trade-offs and synergies among the ESs have been quantified. This research serves as a point of reference for the rational utilization of natural resources and for the protection of the ecological environment within the Karst WNHS. This research aims to assess the following ESs in the Karst WNHS, for the period 2000–2020: water conservation (WC), soil conservation (SC), carbon sequestration (CS), and habitat quality (HQ). Additionally, the objective of this study is to examine the space–time evolution of these ESs. Spearman’s correlation coefficient and spatial auto-correlation analyses were utilized to ascertain the temporal and spatial trade-offs and synergies for each ES. The results of this study indicate the following: (1) Between the years 2000 and 2020, the inter-annual changes in HQ and CS at the Shibing and Libo-Huanjiang WNHS exhibited a relatively stable pattern, with a gradual increasing trend, and in contrast, the inter-annual changes in WC and SC exhibited greater fluctuations; (2) the distribution of ESs is closely linked to land use patterns, and woodland is the most critical land type of the ESs, contributing the most to each ES; (3) population density is negatively correlated with various ecosystem services, while GDP is positively correlated with each ecosystem service; and (4) the Karst WNHS in Southern China demonstrates a clear and significant synergistic relationship between WC and CS. This relationship is primarily characterised by a strong synergistic effect. The synergistic relationship between HQ and WC exhibited a weakening trend within both study areas. Furthermore, the SC demonstrated a spatial trade-off relation with HQ, CS, and WC.

1. Introduction

The term “ecosystem services” is used to describe the natural resources and all of the goods and services of life which support long-term human survival [1]. Concurrently, all of the advantages that people derive from the Es are evident [2]. Ecosystem services are grouped into four categories: provisioning, regulating, supporting, and cultural [3]. ESs interact with each other at a given time and space, exhibiting different levels of trade-offs or synergies [4]. ES trade-offs are cases where an enhancement in one type of ES results in a decrease in another type of ES; ES synergies are characterized by a simultaneous rise or fall of both ESs [5]. Nevertheless, the pursuit of excessive trade-offs or synergies among ESs may give rise challenges in achieving the optimal utility of each type of service [6]. Should human activities persist in their alteration of ESs with the objective of obtaining more freshwater resources, when this will have an effect on other types of services to varying degrees [7]. Given the joint impact of global environmental change and intensifying human activities, the identification of effective measures to maximize the benefits of different types of ES has a significant meaning for enhancing the overall benefits of ESs, promoting regional sustainable development and achieving a mutually beneficial outcome for both human well-being and ESs.

The clarification of the intricate interconnections between different ESs is a crucial preliminary step in the process of determining the most appropriate approach to the sustainability management of several ESs [8]. Howe et al. [9] analyzed the results of a large number of cases in which global ESs have been used for human well-being and found that trade-offs have a record that is almost three times as strong as synergies. Furthermore, they concluded that assuming a mutually beneficial outcome situation is not universally relevant, and that trade-offs are more likely to lead to desirable outcomes. Qiao et al. [10] and Lang and Song [11] conducted ESs studies utilizing the extreme value method and the production possibility boundary method, respectively, with the objective of exploring the trade-offs or synergies among different ESs. Qiao et al. [12] employed the Taihu Lake basin in China to examine the time–space dependence of trade-offs and synergy among the various ESs at different periods and scales of urbanization. The objective was to facilitate decision-making for sustainable development. Li and Luo [13] employed InVEST to evaluate ESs and a geoprobe to investigate the trade-offs and synergistic influences on ESs in southern China. Their findings revealed that the main influencing factors were lithology, altitude, and precipitation. This study enriches the research content in the field of ecosystem restoration in karst areas. The research of trade-offs and synergies in environmental services not only elucidates the linkages between different ESs, but also guides humanity to make more rational use of natural resources [14].

Global karst covers 12% of the continental surface area, with China’s karst reaching 3.44 million km² [15]. China is home to the largest distribution area of karst landforms globally. Concurrently, it is also the region of the globe with the greatest concentration of karst landforms, the most advanced karst development, the most diverse landscape types, and the richest biodiversity [16]. The karsts of Southwest China span four latitudinal zones, from the hilly terrain to the highland valleys. They are home to the most unique types of karst primeval forests in the same latitudinal zone of the world, covering the largest area and with the most typical continuous distribution [17]. World Heritage represents a distinctive and invaluable asset of exceptional universal significance, as acknowledged by UNESCO and the World Heritage Committee [18]. The Karst WNHS fulfils the criteria set out in the World Heritage List, specifically criteria vii and viii [19]. It is a significant part of World Heritage, demonstrating unique geological formations, remarkable aesthetic landscapes, rare ecological environments, and a profound humanistic heritage [20]. Meanwhile, Karst WNHS represents a natural ecosystem with the richest biodiversity and the most complex structure. As an important part of the global ecosystem, it can effectively create more social goods and provide important ESs [21]. It is therefore crucial to investigate the space–time variations of ESs and talk over trade-offs and synergistic states within the Karst WNHS.

In recent years, research on ESs for world heritage sites has attracted considerable attention. Most studies have concentrated on the development of assessment methods and the analysis of influencing factors. However, there is a lack of quantitative tools that can provide a unified description of the ESs in question. Furthermore, there has been little research conducted on the relationships among ESs. The ecological environment monitoring network of the Karst WNHS is not yet fully developed, and the results are currently fragmented. Consequently, the implementation of a comprehensive research programme on ESs is of paramount importance. This research employs InVEST to assess the ESs of Karst WNHS. It also analyses the space–time changes in ESs between 2000 and 2020 in Karst WNHS. Furthermore, it analyses the numerical trade-offs and synergistic relations among the ESs of Karst WNHS, using Spearman’s correlation. In parallel, the bivariate spatial auto-correlation method is applied to examine the spatial variability among the ESs of the Karst WNHS. This approach enables researchers and decision-makers to gain a more profound comprehension of the ecosystem structure and the lasting management of the Karst WNHS.

2. Materials and Methods

2.1. Study Area

The Shibing Karst WNHS is located at Shibing County, Qiandongnan Prefecture, Guizhou Province (Figure 1), with an area of 282.98 km² [22]. It is located within the humid climate zone with subtropical monsoon, characterized by a mild and humid climate. The average temperature throughout the year is around 16 °C, with annual precipitation of 1220 mm, and the average altitude of the area is 912 m [23]. The mean annual population density from 2000 to 2020 was 87.6852/km², while the mean annual GDP is 13,523 ¥ per capita, with the majority of this figure derived from the tourism sector. he Shi Bing Karst is a superlative example of tropical-subtropical dolomite karst on the global scale [24]. Dolomite is a carbonate rock containing calcium and magnesium, which is less soluble than limestone. The distinctive peaks and gorges that have developed under such lithological conditions are a defining feature of the Shibing Karst WNHS.

The Libo-Huanjiang Karst WNHS is located at the confluence of Libo County in Guizhou Province and Huanjiang County in Guangxi (Figure 1). It encompasses an area of 730.16 km² [22], within the subtropical humid monsoon climate zone, with an average year-round temperature of approximately 15 °C and an average annual rainfall of 1752 mm, and the average elevation of the area is 758 m above sea level [25]. The mean annual population density from 2000 to 2020 was 62.0958/km², while the mean annual GDP was 19,517 ¥ per capita, with the majority of this figure derived from the tourism sector. The Libo-Huanjiang Karst is a typical subtropical limestone cone-shaped karst landscape. It is the most typical example of the Peak Tree-Peak Forest Karst, which is an illustration of the evolution of similar karst landscapes across the globe. Concurrently, the Libo-Huanjian and Huanjiang karsts collectively constitute a vast and comprehensive subtropical cone-shaped karst World Heritage ecosystem, which offers an optimal environment for biological growth and reproduction [26].

2.2. Data Sources

This research examines the space–time changes in ESs and the trade-offs and synergistic relations at the Shibing Karst WNHS and the Libo-Huanjian Karst WNHS between the years 2000, 2010, and 2020.

In order to conduct research, it is necessary to have access to meteorological data, topographic and geomorphological data, and land use data (

Table 1. Data source and processing methods.

Data Type	Data Name	Data Sources
Meteorological data	Precipitation	The National Weather Science Data Center (https://data.cma.cn/metadata/#/layerType, accessed on 20 April 2024)
	Temperatures	The National Weather Science Data Center (https://data.cma.cn/metadata/#/layerType, accessed on 20 April 2024)
	Potential Evapotranspiration	The National Earth System Science Data Center (https://www.geodata.cn/main/, accessed on 20 April 2024)
Topography	The World Soil Database	The China Soil Dataset (v1.1) is derived from the World Soil Database (HWSD) (https://www.fao.org/soils-portal, accessed on 8 May 2024).
	DEM	NASA (https://www.nasa.gov/, accessed on 26 April 2024)
Land use	Land Cover	The Geospatial Data Cloud (https://www.gscloud.cn/#page1/1, accessed on 10 April 2024)
Socio-economic Statistics Data	Population Density	Statistical Yearbook Sharing Platform (https://www.guizhou.gov.cn/, accessed on 24 June 2024)
	GDP

). Meteorological data encompass precipitation, temperature, and potential evapotranspiration, derived from spatial station records from the Meteorological Data Center of the China Meteorological Administration (CMA) and spatially interpolated using the Kriging method. Topographic and geomorphological data encompass soil data, digital elevation models (DEMs), and other relevant information. As raster data exhibit varying resolutions, inconsistent coordinates, or offsets, the study was standardized to raster data in order to accurately quantify the trade-offs and synergies among the various ESs and to ensure the scientific validity and accuracy of the study. In order to address the issue of differing resolutions, the resampling tool of ArcGIS10.2 was employed to standardize the resolution to 10 m by 10 m. In order to address the issue of the non-uniformity of coordinates and of coordinate shifts, a variety of geographic proofreading, projection, and other tools for coordinate standardization were employed. Due to the limited spatial extent of the study area, remote sensing images from 2000, 2010 and 2020 were employed as the principal data source. The analysis was conducted using ArcMap 10.2 and ENVI 5.3 software, employing supervised classification and a confidence interval of 95%. The results of the field validation, when considered alongside the findings of the land use interpretation, demonstrate that the latter is aligned with the current land use status of the study area.

2.3. Research Methodology

2.3.1. Cellular Automata (CA)-Markov Model

The CA-Markov model is a widely used tool for prediction of land use and land-cover change [25]. It combines the simulation advantages of the CA model with the temporal prediction advantages of the Markov model [27], enabling the evaluation and simulation of spatial and temporal changes in image elements [28]. Consequently, the CA-Markov model used in this paper to simulate the land-cover change of the Karst WNHS, thereby enabling an assessment of its ecological quality.

CA models are frequently employed to depict the space–time unfolding of intricate systems, and are hence designated as a local discrete-time network dynamics model [29]. The principle of calculation is as follows:

$$S_{t+1}=f\left(S_{\left(t\right)},n\right)$$

(1)

where S is the meta-cellular state, t and t + 1 are the corresponding meta-cellular moments, and n is the distance between the two instants [25].

The fundamental tenet of Markov modelling is the comparison of the subsequent state of affairs with that of the preceding state [30].

$$X_{\left(t+1\right)}=X_{\left(t\right)}{P}_{xj}$$

(2)

where $X_{\left(t+1\right)}$ and $X_{\left(t\right)}$ are the states at moments $t+1$ and $t$. $P_{xj}$ are the transfer matrix.

2.3.2. Water Conservation (WC)

Following the concept of the water balance, the use of InVEST simplifies the confluence process, without distinguishing between surface runoff, loamy runoff, and baseflow. It is assumed that the water produced by the grid reaches the water outlet through any of the aforementioned pathways [12].

$$Y_{xj}=\left(1-{\displaystyle \frac{AE{T}_{xj}}{P_x}}\right)\times P_x$$

(3)

where $Y_{xj}$ is the annual water production, $P_x$ is the annual average rainfall of grid cell x, and $AE{T}_{xj}$ is the actual annual average evapotranspiration of grid cell x on land use type j [31] (Table S4 in the Supplementary Materials).

Given that the selected study areas in this paper are all mountainous, it is essential to consider the influence of topographic factors such as elevation and slope on the water-holding function. In order to achieve this, it is necessary to utilize topographic indices, soil saturated hydraulic conductivity, and flow rate coefficients in order to correct the water yield [32]. The final calculations yielded a water conservation which is as follows:

$$Retention=min\left(1,{\displaystyle \frac{249}{Velocity}}\right)\times min\left(1,{\displaystyle \frac{0.9\times TI}{3}}\right)\times min\left(1,{\displaystyle \frac{Ks}{300}}\right)\times Yield$$

(4)

where Retention is the amount of WC, TI is the tectonic indicator, Ks is the soil saturated hydraulic conductivity, Velocity is the flow coefficient, and Yield is the annual water yield [33] (Figure S1 in the Supplementary Materials).

2.3.3. Soil Conservation (SC)

The sediment transport rate module of InVEST was used to assess conservation [34]. The calculation principle is as follows:

$$Sre{t}_x=\left(RKL{S}_x-USL{E}_x\right)+up{s}_{retai{n}_x}$$

(5)

where $Sre{t}_x$ is the amount of sediment retained by the grid (t) and $up{s}_{retai{n}_x}$ is the amount of sediment intercepted by the grid upstream of the grid (t) (Table S4 in the Supplementary Materials).

2.3.4. Habitat Quality (HQ)

The HQ component of InVEST was applied to evaluate the level of biodiversity, with the HQ Index serving to indicate the quality of habitats within a given area [35]. The index is a function of the intensity of the impact of the stressor (Table S1 in the Supplementary Materials), the sensitivity of the habitat to the stressor (Table S2 in the Supplementary Materials), the distance of the habitat from the stressor, and the degree of the legal protection of the land [36].

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

(6)

where $Q_{xj}$ and $D_{xj}$ are the HQ and the level of stress suffered by grid x in habitat type j, respectively; $H_j$ is the habitat suitability index; z is the normalized constant; and k is the semi-inclusive sum constant [37].

2.3.5. Carbon Storage (CS)

InVEST employs a specific methodology for calculating carbon stocks. This methodology involves the creation of a table of the four major carbon pools, with vegetation type serving as the statistical unit (Table S3 in the Supplementary Materials). The four major pools are then aggregated for different land cover types, and finally, the carbon stock in a given area is determined [38].

$$C_x={\displaystyle \sum }_{j=1}^J{A}_{xj}\left(C_{aj}+C_{bj}+C_{sj}+C_{dj}\right)$$

(7)

where $C_x$ is the carbon stock in area x and $A_{xj}$ is the range-of-land-use category j in zone x. $C_{aj}, C_{bj}, C_{sj}, C_{dj}$ represent the surface biomass, beneath the surface biomass, soil carbon sink, and dead organic matter carbon pool, respectively [39].

2.3.6. Spearman Correlation Coefficient

Correlation analyses are a common method of assessing interactions among ESs [40]. This paper investigates the synergies and trade-offs among ESs in the region using correlation analysis. Spearman’s correlation is a non-parametric method of statistical analysis, and geographic data are non-linear and non-normally distributed [40].

$$\mathsf{\rho }\left({\mathrm{F}}_{\mathrm{Xi}},{\mathrm{F}}_{\mathrm{Yi}}\right)=1-{\displaystyle \frac{6{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{F}}_{\mathrm{Xi}}-{\mathrm{F}}_{\mathrm{Yi}}\right)}^2}{\mathrm{n}\left({\mathrm{n}}^2-1\right)}}$$

(8)

where $\mathsf{\rho }\left({\mathrm{F}}_{\mathrm{Xi}},{\mathrm{F}}_{\mathrm{Yi}}\right)$ stands for the Spearman correlation coefficient between ${\mathrm{F}}_{\mathrm{Xi}}$ and ${\mathrm{F}}_{\mathrm{Yi}}$. ${\mathrm{F}}_{\mathrm{Xi}}$ and ${\mathrm{F}}_{\mathrm{Yi}}$ denote the land cover function. Xi and Yi in the space unit i. n denote the quantities of space units [41]. When $\mathsf{\rho }$ > 0, the two services are positively correlated and synergistic; conversely, when $\mathsf{\rho }$ < 0, the decrease in one type of ES is negatively correlated with an increase in the other type of ES and is a trade-off [42].

2.3.7. Bivariate Spatial Autocorrelation Model

Spatial correlation analysis is used to determine and evaluate patterns of spatial correlation, with the objective of detecting heterogeneity or convergence in spatial data. This encompasses the application of Global Moran’I and Local Moran’I [43].

(1) Global spatial autocorrelation is employed to estimate the overall level of linkage, significance, and geographical spread of ESs across all study units within the study area [4]. Moran’s I range if I > 0 indicates that the spatial distribution of ESs is positively correlated and that there is a ‘high-high’ or ‘low-low’ aggregation phenomenon; if I = 0, it indicates that the grid units are randomly distributed in space; if I < 0, it means that the ESs are negatively correlated in space, and that there is the phenomenon of a ‘high-and-low’ or a ‘low-and-high’ neighborhood [44].

$$I={\displaystyle \frac{n{{\displaystyle \sum }}_{i=1}^n{{\displaystyle \sum }}_{j\ne i}^n{\omega }_{ij}(y_i-\overline{y})·\left(y_i-\overline{y}\right)}{\left({{\displaystyle \sum }}_{i=1}^n{{\displaystyle \sum }}_{j\ne i}^n{\omega }_{ij}\right){{\displaystyle \sum }}_{i=1}^n{(y_i-\overline{y})}^2}}$$

(9)

where I is the global autocorrelation Moran index; $y_i$ and $y_j$ denote the ESV values of the ith and jth spatial grid cell, respectively; $\overline{y}$ is the average value of ESV of the evaluated grid cell; ${\omega }_{ij}$ is the weight matrix established based on the spatial grid adjacency; and n is the number of grids of the spatial cell.

(2) Local spatial autocorrelation is a statistical technique that can be employed to test the spatial correlation and spatial heterogeneity of ES functions across local regions [45]. Furthermore, the bivariate LISA (Local Indicators of Spatial Association) distribution map was used to characterize the correlation among ESs within the study area [41]. The calculation principle is as follows:

$$I_i={\displaystyle \frac{\left(x_i-\overline{x}\right){{\displaystyle \sum }}_{j=1}^n{\omega }_{ij}\times \left(x_j-\overline{x}\right)}{{{\displaystyle \sum }}_{j=1}^n{\left(x_i-\overline{x}\right)}^2/n}}$$

(10)

where $I_i$ denotes the local spatial autocorrelation Moran index; $x_i$ and $x_j$ are the ESV values of the $i$ and $j$ spatial cells, respectively; $\overline{x}$ denotes the average value of the ESV of the lattice grid; and ${\omega }_{ij}$ and n have the same meanings as in Equation (9).

3. Results

3.1. Characteristics of Land Use Conversion

The following figures illustrate the land use spatial distribution map and modifications of the Shibing and the Libo-Huanjian between the years 2000 and 2020 (Figure 2). The predominant land use type of the Shibing Karst WNHS is woodland (WL), with the proportion of WL area fluctuating between 60% and 80% across the three periods. This is characterized by a wide range of distribution and a considerable area. The cropland (CL) and building land (BL) are primarily situated within the buffer zone and are distributed in blocks. In respect to the amount of change from 2000 to 2020, the area of WL exhibited the greatest change, with a cumulative increase of 47.76 km². Unused land (UL) demonstrated a tendency to increase and then reduce, while the CL remained in decline. The range of grassland (GL), BL, and water (WB) areas did not change significantly.

The predominant land use type of the Libo-Huanjian Karst WNHS from 2000 to 2020 was WL. The proportion of WL area in all three periods ranged from 80% to 90%, with the distribution including the heritage site and the butterfly area. The remaining land types are primarily concentrated in the buffer zone in a linear distribution. In terms of inter-annual changes in the various types of areas from 2000 to 2020, the area of WL has undergone the most significant change, increasing by up to 54.92 km². This is followed by UL, which has undergone a cumulative decrease of 34.48 km². The majority of the UL has been transformed into GL and WL under the protection of localities, as evidenced by the spatial change map. The area of GL has exhibited a gradual increase, while the areas of CL and BL have shown a gradual decline. In contrast, the region of WB has exhibited minimal change.

The land use transfer matrix of the Shibing and the Libo-Huanjiang from 2000 to 2020 indicates that the majority of the transfer of the two study areas occurred between WL, followed by UL, BL and CL, with less transfer in and out of GL and WB areas. In summary, the Shibing and the Libo-Huanjiang are characterized by a high degree of WL resources.

3.2. Spatiotemporal Heterogeneity of ESs

3.2.1. Spatiotemporal Variations in ESs

From 2000 to 2020, the WC of the Shibing Karst WNHS exhibited a general trend of increasing and then decreasing. The WC decreased from 540.33 mm in 2000 to 476.186 mm in 2010, and then increased to 571.921 mm in 2020. The SC demonstrated a general trend of decreasing and then increasing, with the WC decreasing from 4487.7 t/(km²·a) in 2000 to 4590.93 t/(km²·a) in 2010, and then increasing to 6237.06 t/(km²·a) in 2020. The HQ index demonstrated an upward trajectory, expected to grow by an average annual 9.89%. The mean value of HQ increased from 0.711 to 0.778 between 2000 and 2010, expected to grow by an average annual 4.63%. The mean value of HQ in the 2010–2020 period exhibited a continuous upward trend, reaching an average value of 0.851 by 2020. The CS change demonstrates a fluctuating pattern, with a downward trend followed by an upward one. The average rate of change of CS from 2000 to 2010 was −0.07%, with minimal fluctuations in the average value. In contrast, the average rate of change of CS from 2010 to 2020 was 7.03%, accompanied by a notable increase in the average value, reaching 2.4953 by 2020.

The WC of the Libo-Huanjiang Karst WNHS shows an overall decreasing trend from 2000 to 2020. The WC decreased from 527.846 mm in 2000 to 476.186 mm in 2010, and then it reduced to 408.581 mm in 2020. The SC demonstrated an overall decreasing trend, with an upward trend and then a trend of decline. The WC decreased from 5856.83 t/(km²·a) to 4590.93 t/(km²·a) in 2010, and then it increased to 8399.63 t/(km²·a) in 2020. The HQ trend is comparable to that of Shibing, exhibiting a gradual increase in the trend of its average annual increase of 1.849%. From 2000 to 2010, the HQ increased from 0.798 to 0.813, representing a growth rate of 1.86%. Similarly, from 2010 to 2020, the average value of HQ increased from 0.786 to 0.808, with an average annual increase of 1.82%. From 2010 to 2020, the average value of HQ continued to increase, reaching 0.885 by 2020. The CS changes exhibited a tendency to fall and then to rise again, with an average rate of change in carbon storage of −0.11 per cent from 2000 to 2010 and an average rate of change in CS of 5.39 per cent from 2010 to 2020, with the average value reaching 2.5059 by 2020.

With regard to the spatial distribution of WC and CS in the Shibing and the Libo-Huanjiang, the overall pattern has changed little. The high values are focused on the WNHS, while the low values are found in the construction land and the water areas in the buffer zones. The high values of HQ in the two study areas are concentrated in the woodland and grassland areas within the heritage sites and buffer zones. In contrast, the low values are primarily focused on the CL and BL areas. The high values of SC in the Shibing are concentrated in the northern region, with a gradual decrease in the south. In contrast, the SC in the Libo-Huanjiang is concentrated in the southeast, with a corresponding decrease in the south. The negative values are distributed primarily in the buffer zone (Figure 3).

3.2.2. Changes in ESs Dynamics of Different Land Use Types

Over the past two decades, the WC, HQ, and CS of the Shibing Karst WNHS have demonstrated a general stability. The buffer zone has undergone a significant change and demonstrated a downward trend. In contrast, SC exhibited an upward trend, with a few areas displaying a linear downward trend. The land classes with the highest to lowest mean annual WC per unit area, in descending order, are as follows: WL, GL, CL, UL, BL, and WB. WC accounted for more than 70 per cent from 2000 to 2020, and for 89.59 per cent in 2020. WL, CL, GL, UL, BL, and WB are the land categories with the highest annual average SC per unit area, in descending order. The proportion of SC was above 80% from 2000 to 2020, and reached 93.16% in 2020. The ranking of the land categories by annual average HQ per unit area, from highest to lowest, are as follows: WL, GL, UL, CL, BL, and WB. WL accounted for more than 88 per cent from 2000 to 2020, with a peak of 91.43 per cent in 2020. In contrast, WB accounted for less than 1 per cent. The annual average CS per unit area is as follows, in descending order: WL, UL, CL, GL, BL, and WB. WL accounted for more than 70 per cent annually from 2000 to 2020, with all categories demonstrating a gradual growth trend.

Over the past two decades, the WC, HQ, SC, and CS of the Libo-Huanjiang Karst WNHS have demonstrated a consistent overall trend, with the SC exhibiting a linear decline in certain areas of the site. The different ESs exhibited a substantial decline in the buffer zone, demonstrating a block-like pattern of reduction. The land classes with the highest to lowest mean annual WC per unit area, in descending order, are as follows: WL, CL, GL, UL, BL, and WB. WC accounted for more than 80 per cent of the total from 2000 to 2020, and 90.16 per cent in 2020, while the other land categories accounted for less than 6 per cent of the total. The land categories, ranked from largest to smallest in terms of annual average SC per unit area, are as follows: WL, CL, GL, UL, BL, and WB. The proportion of SC from 2000 to 2020 was above 80%, reaching 93.16% and 92.25% in 2020. The annual average HQ per unit area, in descending order, are as follows: WL, UL, CL, GL, BL, and WB. WL accounted for more than 80% from 2000 to 2020, with a peak of 91.24% in 2020. This demonstrates a sustained growth trend over the past 20 years. The annual average CS per unit area is, in descending order, WL, UL, GL, CL, BL, and WB. WL accounted for more than 70% of the total annually from 2000 to 2020, with all categories demonstrating a gradual growth trend.

It can be observed that different land cover types are capable of providing different types of ESs. In this study, the four services WC, SC, HQ, and CS are further analyzed from the perspective of land cover types. Annual average WC and SC per unit area of different land cover types in the Shibing Karst WNHS and the Libo-Huanjiang Karst WNHS exhibited considerable variation, particularly between forest land and other land cover (Figure 4).

3.3. ESs Trade-Offs/Synergies

3.3.1. Spearman Correlation Coefficients between ESs

Among the six pairs of ESs in the Shibing Karst WNHS, WC and CS exhibited a highly correlated relationship (r ≥ 0.5), indicative of a strong synergistic effect, with a mean value of a correlation coefficient of 0.88. WC exhibited a weak synergistic relationship with HQ, while HQ demonstrated a similarly weak synergistic relationship with CS. The mean correlations are 0.38 and 0.37. Additionally, a weak trade-off was observed among WC and SC, HQ and SC, and CS and SC. The mean values of the correlation coefficients are as follows: −0.074, −0.076, and −0.15 (Figure 5).

Among the six pairs of ESs in the Libo-Huanjiang Karst WNHS, WC and CS in 2000 and WC and CS in 2020 were all in strong synergistic relationships; the coefficients are 0.85 and 0.8, respectively. WC and HQ in 2000 and HQ and CS were all in weak synergistic relationships, with correlation coefficients of 0.13 and 0.18. WC and SC, HQ and SC, and CS and SC are all weak compromises. The coefficients are −0.012, −0.048, and −0.14. In 2010, CS and WC, SC and WC, and HQ and CS were all weakly synergistic, with coefficients of 0.057, 0.037, and 0.11. WC and HQ, SC and HQ, and CS and SC are all weak trade-offs; the coefficients are −0.02, −0.13, and −0.058, respectively. In 2020, WC had a weak synergistic relationship with SC, and HQ had a weak trade-off relationship with SC, and SC had a weak trade-off relationship with HQ, with SC, and with CS; the coefficients are 0.0032, −0.14, and −0.14 (Figure 5).

3.3.2. Moran’s Index between ESs

The scatterplot is divided into four quadrants, and each is centred on the mean value (Figure 6). The statistical analysis of each quadrant revealed that the two study area points in the graph are primarily distributed in quadrants 1 and 2. Among the indices, the global autocorrelation indices of CS and WC, HQ and WC, and HQ and CS of the Shibing from 2000 to 2020 were significantly positive at the level of 0.05. This indicates a positive spatial correlation, indicating synergy among the ESs. Furthermore, the synergistic strengths of WC and CS and of WC and HQ are decreasing and then increasing, while the synergistic strengths of HQ and CS are continuously increasing. The global autocorrelation indices of SC with WC, SC with HQ, and SC with CS are significantly negative at the level of 0.05, indicating a negative spatial autocorrelation. This suggests that the ESs are in a trade-off. The strength of the trade-offs between SC and WC and SC and CS is weakened and then strengthened, while the strength of the trade-offs between SC and HQ is strengthened and then weakened. In 2020, the strength of the trade-offs between SC and HQ reached a synergistic state, as indicated by the Moran’s index, which is −0.005.

The global autocorrelation indices of WC and CS and of HQ and CS in the Libo-Huanjiang Karst WNHS exhibited a positive association at the 0.05 level during the period of 2000–2020. This indicates a positive spatial correlation, which suggests a synergistic relationship between the services. The synergistic relationship between WC and CS tended to decrease and then increase, while the relationship between HQ and CS gradually increased during the period of 2000–2020 In contrast, SC and HQ show a trade-off, with the level of trade-off gradually increasing. A positive spatial correlation has been identified between WC and CS, which suggests a synergistic relationship between the two services.

3.3.3. Spatial Representation of Trade-Offs and Synergistic Relationships for ESs

The results indicate that (Figure 7), among the six pairs of ESs in the Shibing Karst WNHS, the synergies between CS and WC, HQ and WC, and HQ and CS was more significant. Furthermore, the low-low clusters were mainly distributed within WNHS. The high-low aggregation area is focused on the cultivated land and the construction land regions. The spatial distribution of the high-first aggregation and low-high aggregation in the trade-off relationship among SC and WC, SC and HQ, and SC and CS is more dispersed. Among the three pairs of services, the spatial spread of the trade-off in 2020 exhibited the widest scope of distribution and the largest number of occurrences.

Among the six pairs of ESs in the Libo-Huanjiang Karst WNHS (Figure 8), the synergy between WC and HQ and between HQ and CS is more significant. The low-low clustering is mainly distributed in the part of the Huanjiang Protected Area in a block-like distribution, while the high-high clustering is mainly distributed in the buffer zone of arable land, construction land, and unused land area in a shape-like distribution. The SC and WC, SC and HQ, and SC and CS Weighing Guanzhong high and low clusters are concentrated in the core region of the Huanjiang WNHS, with the low and high clusters primarily situated in the Libo buffer zone.

4. Discussion

4.1. Spatial and Temporal Changes in Land Use and ESs

In recent years, the delicate karst ecosystem has undergone significant alterations to the karst landscapes, which have been caused by both climate change and human activities [15]. Nevertheless, land-cover change also has a direct influence on the temporal development and the spatial extent of ESs [46]. This paper employs a spatio-temporal analysis of ESs in the Shibing and Libo-Huanjiang Karst WNHS. The findings indicate that WC is related to precipitation, soil depth, and the available water content of vegetation. The heritage site is predominantly forest and grassland, with fertile soil and a high vegetation cover, resulting in greater soil depth and vegetation-available water content than the buffer zone. Consequently, the spatial distribution of WC is more prevalent in heritage sites than in buffer zones. Changes in the area of woodland and grassland also have a major impact on soil retention [47]. The interception of rainfall by woodland and grassland can be attributed to the forest canopy, the leaf litter, and to the grassland layer on the soil surface. This phenomenon reduces the direct scouring of the soil by rainwater, thereby increasing SC. These results are in accordance with those previously reported by Chen et al. [48].

A multitude of researchers have demonstrated that land-cover change has a profound impact on the structure, the characteristics, and the processes of ecosystems [49,50,51]. The CS in the region area was primarily focused on land classes with high vegetation cover, such as WL and GL. WL, as the land category with the highest CS among all land categories, demonstrated a significant increase in CS with the expansion of the WL area. These findings align with those of previous studies, including those by Qiao et al. [10] and Yuan et al. [52]. The HQ index is primarily indicative of the local ecological status and degree of biodiversity [53]. The buffer zone within the study area is predominantly composed of CL and BL, which are not conducive to the survival of organisms. Consequently, the HQ is relatively low. The HQ was higher in the heritage site due to the presence of WL, GL, and other ecological land types, which provided an optimal environment for the survival of organisms. These results are in accordance with those previously reported by Lin et al. [54] and Cui et al. [55]. In summary, the reduction in CL and BL and the increase in WL between the years 2000 and 2020 have resulted in an overall increase in CS and HQ. Consequently, alterations in land cover structure will predominantly give rise to analogous alterations in Ess’ functions.

4.2. Trade-Offs/Synergies of ESs

The innovation of this research lies in the application of the Spearman correlation coefficient to calculate ES trade-offs and synergies. This approach has been validated both in general and spatially using bivariate and local Moran’I indices, which have proven to be highly accurate. The correlation calculated by the Moran’I index was found to exhibit a certain lag, which reflected the ES’s trade-offs and synergies in a more realistic and comprehensive manner. The WC and the SC are significantly affected by precipitation. In contrast, the HQ and the CS are strongly influenced by land use change [56]. The extensive and continuous growth in the region of woodland in the research region, which has more than 75% vegetation cover, contributed to the improvement of HQ and CS. This resulted in a synergy between HQ and WC, CS and WC, and HQ and CS, which is more consistent with the results of existing studies [57]. The monoculture of vegetation in the core area, which is mainly forested, and the high distribution of BL and CL in the buffer zone, resulted in the uneven rain erosion resistance of the soil in the study area. This, in turn, indirectly affected the SC capacity, in accordance with the findings of Chang et al. [58] and Tian et al. [59]. Furthermore, ESs exhibit scale effects [60]. The outcomes of trade-offs and synergies can be influenced by a number of factors, including spatio-temporal heterogeneities, natural environmental factors, and socio-economic conditions.

An increase in population density and the clustering of economic activities results in a notable surge in demand for ecosystem services, which in turn precipitates ecosystem degradation and a reduction in the supply of ESs [61]. The implementation of appropriate regulatory measures to control the intensity of human activities is conducive to the synergistic development of ecosystem services. The Shibing Karst WNHS and the Libo Karst WNHS both demonstrate a negative correlation between population density and various ecosystem services, as well as a positive correlation between GDP and each ecosystem service (Figure 9). The results indicated that population density and GDP exhibited the most significant correlations with soil conservation, habitat quality, and carbon storage, with correlation coefficients exceeding 0.5. An increase in population density is associated with a reduction in ecosystem services. This illustrates the importance of considering the spatial variability of ecosystem services in response to changes in socio-economic factors. In the future, changes in different ecosystem services can be studied using longer-time-span data and regionally varying socio-economic parameters, thus providing a more comprehensive understanding of the factors influencing these services.

4.3. Limitations and Prospects

This research analyses the changes that have occurred in the Karst WNHS over a period of thirty years. The conclusions reached will have certain reference value for the more detailed examination of the relationship between ESs in the study area. However, there may be some differences compared with the long-term series of spatial and temporal dynamics studies. In the future, we will select continuous time scale data to conduct a multi-year continuous change study in order to reduce the uncertainty. Furthermore, the influence of precipitation, temperature, and topography on ESs is more complex. In the future, we will investigate the factors and mechanisms that drive the different ESs in depth, thus providing a foundation for the optimization of the distribution of land cover and the sustainable development of ecological and economic ecosystems in the Karst WNHS in the future.

As a globally significant ecological functional area, the objective of this paper is to propose a framework for promoting synergistic relationships among ESs within the site and to enhance the optimization of the ecological environment. This framework is based on the results of the research conducted on the Karst WNHS:

(1)

Rationalize the planning of land use types; continue to implement ecological conservation projects, for example restoring farmland to woodland and grasslands; and relocate some of the arable and building land to the townships around the buffer zone, with the objective of enhancing the overall ecological suitability of the heritage site and improving the sustainability of its ecosystem;

(2)

The objective is to optimize zoning management, to manage the existing heritage sites and buffer zones in a zoning pattern, and to integrate them with the space–time distribution of ESs and the spatial pattern of global trade-offs and synergies, as assessed in the present study. Furthermore, targeted zoning management and protection of the ecological service areas that require upgrading and the areas of the main trade-offs and relationships will be carried out. This will facilitate the ecological protection of the Karst Natural WNHS, while also contributing to global sustainable development.

5. Conclusions

In order to enhance the regional ESs’ provisioning capacity, to mitigate conflicts between Ess, and to achieve the sustainable development of WNHS, it is crucial to investigate the synergies and trade-offs among ESs in the Karst WNHS. In this research, the InVEST model was used to assess the four key ESs in the Shibing Karst WNHS and the Libo-Huanjiang Karst WNHS between the years 2000 and 2020. The changes in the development of ESs in different land use types were identified by histograms. The trade-offs and synergies among the ESs were analyzed with Sperman correlation analysis and a bivariate spatial autocorrelation model, and the following conclusions were reached:

(1) The inter-annual changes in WC and SC in the Shibing WNHS show considerable variability, displaying a declining and then increasing trend. At the Libo-Huanjiang WNHS, there is a gradual decrease in WC, a decrease followed by an increase in SC, and an uninterrupted increase in HQ and CS;

(2) The distribution of high-and low-value areas for ESs is closely linked to land use patterns. The WNHS is characterized by a high degree of vegetation and a rich biodiversity, which represents the high value area of each service. The buffer zone comprises CL, UL, and BL, which collectively represent the low value distribution area for each service. It has been demonstrated that WL represents the most critical land type in the development of ESs in the study area, with the greatest contribution to each service;

(3) A significant correlation was observed between changes in population density and GDP and changes in ecosystem services. The strongest correlations were found between population density and GDP and soil conservation, habitat quality, and carbon storage;

(4) The trade-offs and synergies among ESs were relatively stable, with the exception of WC and SC, which exhibited considerable fluctuations over time. There are significant synergistic relationships between CS and WC, HQ and WC, and HQ and CS in terms of their spatial distribution. Conversely, there are significant trade-offs among SC and WC, HQ and SC, and SC and CS in terms of their spatial distribution.

The clarification of the trade-offs and synergies among ESs and their spatial distribution represents the foundation for an examination of the ES’s changes in natural World Heritage sites. This investigation will facilitate the macro-control of natural resources in the sites and will serve as a reference point for the designation of a coordinated development strategy for resources and the environment in the natural World Heritage sites around the globe.