Introduction

Metacoupling, first proposed by Liu in 2017, encompasses human-nature interactions within a particular system (intracoupling), between adjacent systems (pericoupling), and between distant systems (telecoupling), providing a new approach to analyze socio-economic and environmental interactions within human-nature coupled systems across spatial scales1. Traditional approaches (e.g., intracoupling2, telecoupling3, etc.) usually focus on a single dimension and are limited to the impacts of human-nature interactions within a system or across a region, which the metacoupling framework compensates for by simultaneously considering interactions and relationships within a region and with neighboring and distant regions4. The framework has been widely applied in environmental-economic studies and provides a solid basis for analyzing physical, non-physical, and virtual flows across different distances5. For example, Du et al.6 revealed the impact of global virtual water trade on local water stress through a metacoupling; Carlson et al.7 revealed the interrelationships of marine fishing at local, regional, and global scales through metacoupling. Introducing a metacoupling framework for research has become an important direction for promoting synergistic development of the environment and economy and optimizing ecosystem service assessment globally. The study of carbon sink service flows (CSSFs) based on the metacoupling framework can provide new theoretical perspectives for the cross-regional assessment of ecosystem services.

CSSFs, driven by anthropogenic and natural factors, are transferred from supply regions to demand regions via specific pathways and media8. In recent years, rapid population growth, the expansion of construction land, and an increase in consumer demand9,10have exacerbated an imbalance between the supply and demand of the carbon sink service (CSS)11,12. CSSF can effectively connect the supply and demand relationship of spatial changes, quantify and draw the spatial flow of CSS, and become the focus of future research13. Luo et al.14 analyzed how land use affects the flow of carbon sinks(CSs); Wu et al.15 identified the transfer path between supply and demand areas. We study the whole process of CSS from generation to flow and use, which helps to identify how CSS moves and provides theoretical support for regional ecological protection and sustainable development.

Research on CSSF is actively evolving, with incomplete quantitative assessment methods16. Current CSSF assessment methods are divided into physical quantity-based (e.g., ARIES17and EcoMetrix models18, etc.) and value-based (e.g., ecological footprint method19, MRIO model20, breakpoint model21, and field strength model22, etc.) approaches. The breakpoint model begins with economic impacts between cities and can be used for various ecosystem service flow measurements because of its similarity to ecosystem service flow theory. Now have been a large number of studies based on the fracture point model for ecosystem service measurements, such as Yu et al.23, Ji et al.24, Guan et al.25, Zhai et al.26, and others. The physical quantity-based method is limited by multi-source data and scale applicability issues27. In contrast, the value quantity-based method has broader applicability. Most current studies focus on a single perspective and scale, overlooking the spatial coupling relationship between nonneighboring areas and the subsystems of composite ecosystems. The reach of CSS can extend to regional, national, and global levels. Zhang et al.28 used the equivalent factor method and the breakpoint model to analyze the ecosystem service flow inside and outside the Huangshui Basin under the framework of metacoupling. Vergara et al.29 investigated the spatial connections among cities, regions, and neighboring regions in southern Chile’s marine space via a metacoupling framework. Da Silva et al.30 evaluated the impacts of intracoupling and intercoupling on the environmental and socioeconomic factors of soybean production in Mato Grosso, Brazil, via a metacoupling framework. Most of these studies focus on the water resource supply and grain economic trade. However, there is little research on CSSF at home and abroad, and there is a lack of in-depth studies on the coupling mechanism of supply, flow, and demand of CSSs under the metacoupling framework.

The shortage of CSS supply and demand is increasingly critical, and conducting an assessment and pathway modeling of CSSF based on “supply‒flow‒demand” has become a vital method to solve this issue. Most current pathway modeling methods for ecosystem service flows focus on water provisioning services31,32,33, and methods based on CSSF include ARIES8, Service Path Attribute Networks (SPANs)17, and RI34, etc. SPANs model can simulate the supply, use, and flow processes of a CSS under ideal and difficult-to-access conditions, becoming an essential focus in spatial CSSF research35. However, this model is currently concentrated in North America and Europe due to the large number of required parameters and complex data36,37. By simplifying the SPAN model and adopting its conceptual framework of “source‒sink‒use”, many researchers have quantified the transfer process of CSSF by combining it with Bayesian belief networks (BBNs) and ecological models. There are limited studies on uncertainty and sensitivity analysis of SPANs‒BBNs models38,39, and it is possible to combine SPANs with BBNs models and utilize the sensitivity analysis of BBNs to assess the performance of SPANs models in different uncertainties and stability in ecosystem service assessment40.Liu et al.41 simulated water resources in Inner Mongolia via the SPANs model, classified variables affecting water supply services via the BBNs model, and calculated the state probability. Yang et al.42 drew the direction and trend of regional CSS spatial flow by using the SPANs model and GIS tools through wind direction elements. The “source‒sink‒use” process of the CSSF can be effectively simulated by the SPANs‒BBNs model.

The composite ecosystem of mountains, rivers, forests, farmlands, lakes, and grasslands (CEMRFFLG) shoulders the task of establishing an ecological protection compensation mechanism that reflects the value of the CS. Establishing a robust ecological compensation system for the CSS is crucial for achieving the “dual carbon” goal. The CEMRFFLG is a crucial CS, and protecting and restoring it enhances the CSS. Current research on CSS focuses mainly on forest CS43, grassland CS44, urban CS45, and ocean CS46, which cover a wide range of fields. Xu et al.47 assessed the CS capacity of China’s forest, grassland, farmland, wetland, and shrub ecosystems. Lou et al.48 assessed the carbon emission (CE) and sink capacities of China’s forests, grasslands, and farmlands. Existing studies have ignored the differentiation and spatial and temporal change characteristics of CSSs among different ecosystem subsystems of the CEMRFFLG, and there are fewer studies on the flow of CSSs in the CEMRFFLG. The research perspective is informative for the study of the metacoupling process of supply-flow-demand of global composite ecosystem services.

Researchers are increasingly focusing on CSSF. The study of the CSSF in the CEMRFFLG raises two key questions: (1) How can we identify the spatial dislocation of the CSS supply and demand in the composite ecosystem and analyze the coupling relationships among the CSS supply, flow, and demand in the CEMRFFLG? (2) How can we create a comprehensive multiscenario evolutionary model of a CEMRFFLG and analyze the mechanisms governing the spatial flow of the CSS in different scenarios?

Based on the quantitative assessment of the supply and demand of CSSs in the CEMRFFLG, we will reveal the metacoupling mechanism of “supply-flow-demand” in the composite ecosystem, and explores the distribution law of CSSF among different subsystems; and based on the spatial mismatch between the supply and demand of CSSs in CEMRFFLG, multi-scenario simulations were conducted to reveal the flow mechanism of CSSs, and provide a scientific basis for optimizing the management of CSSs. Based on the above research objectives, this paper takes Chongqing as the study case area, and quantitatively assesses the supply and demand of CSSs in the CEMRFFLG through the Carnegie‒Arms‒Stanford approach (CASA) model and the IPCC inventory method; explores the coupling process between different subsystems through the breakpoint model combined with the metacoupling framework; and simulates the flow process of the CSSs based on the SPANs-BBNs model, to reveal the flow mechanism is revealed. The study results can improve the understanding of the objective law regarding the spatial flow of the CSS in Chongqing and facilitate the green development of the ecosystem composite pattern.

Materials and methods

Research framework

In this paper, we select four subsystems, namely, water, forest, farmland, and grassland, to measure the total supply and demand of the CSS and analyze the CSSF between the systems by establishing a focal‒adjacent‒distant coupling framework. This paper’s focal system is “one district‒two groups”; the adjacent system is Chongqing, and the distant system is the Yangtze River Basin (YRB). Zhang et al.28 used the abundance of water resources in the Yellow River as a criterion for judging the supply and demand zones; Liang et al.13 took into account the area of different regions as well as their locations in the watershed when determining the supply and demand zones of carbon sequestration services; Yan et al.49 determined the spatial distribution of ecosystem service supply and demand by examining different area differences in land use types to determine the spatial distribution of ecosystem service supply and demand. Depending on the subsystem’s size, the supply area for the adjacent and distant systems is selected in “one district‒two groups”. The Chongqing metropolitan area (CM) is the focal system of the farmland and water subsystem, the southeast Chongqing urban agglomeration (SC) is the focal system of the forest subsystem, and the northeast Chongqing urban agglomeration (NC) is the focal system of the grassland subsystem.

The research framework of this paper is shown in Fig. 1. In the first part, the CASA model and the IPCC inventory method are used to obtain the total CSS supply and demand of Chongqing and the YRB, respectively. The second part uses the breakpoint model to quantify the CSSF between the four subsystems in focal‒adjacent‒distant systems. The third part uses the SPANs‒BBNs model to simulate the space flow path of the CSS. By dividing the four subsystems of the CSSF into service clusters, some suggestions are proposed to optimize the spatial pattern of the CSSF in Chongqing.

Fig. 1
figure 1

Research framework.

Full size image

Research area

Chongqing is located in Southwest China, at the heart of the Three Gorges Reservoir area (Fig. 2a). According to the Third National Land Survey, Chongqing spans approximately 8.24 million hectares, with farmland and forest comprising 22.70% and 56.91% of the total area, respectively. The region has 38 districts and counties, forming a development pattern known as “one district–two groups”. This pattern includes the CM, NC, and SC. The terrain is complex, with hills and mountain ranges covering approximately 90% of the area. To achieve the “dual carbon” goal, Chongqing reduced its carbon energy intensity by 19.4% during the “13th Five‒Year Plan period”. Owing to its rapidly developing economy, the CE remains high, and the imbalance between the supply and demand for the CSS needs to be addressed.

This paper divides the CEMRFFLG into four subsystems (Fig. 2b): water, forest, farmland, and grassland. These four subsystems are not individual elements but are interconnected, influencing each other to create a complex ecosystem. The water subsystem consists of river, lake, and beach elements. The forest subsystem includes mountain and forest elements. Chongqing is a typical mountain city, with mountains and hills covering approximately 76% and 18% of the area, respectively. Most forests grow on mountains, so mountains and forests are placed in the same system. The farmland subsystem consists of paddy fields and dryland elements. The grassland subsystem includes grassland elements.

Fig. 2
figure 2

Map of the study area. (a) Study area (This map depicted in the figure was generated by the authors using ArcGIS 10.4 (Esri, URL: https://www.esri.com/en-us/arcgis/products/arcgis-desktop/overview) and does not require any license). (b) The four subsystems of this paper.

Full size image

Data source

This paper uses three main data types for 2000, 2010, and 2020. The first type of data is vector data, including administrative boundary data; the second type of data is raster data, including DEM, NDVI, LUCC, POP, and GDP data. When the resolution of the original data is 1 km, it is necessary to resample the data to 30 m by the resampling tool in ArcMap software when evaluating Chongqing; and the third type of data is statistical data, including energy consumption data (from the Chongqing Statistical Yearbook (https://tjj.cq.gov.cn/), and China Energy Statistical Yearbook (https://data.cnki.net/)), and meteorological data (solar radiation, precipitation, temperature, evapotranspiration) (from the National Meteorological Information Center‒China Meteorological Data Network (https://data.cma.cn/)). The first and second data types are from the Resource and Environmental Science Platform (https://www.resdc.cn/). In this paper, the research accuracies of Chongqing and the YRB are not consistent, so they are calculated separately. The resolution of the data in Chongqing is unified as 30 m, and the data in the YRB are unified as 1 km.

CSS supply model

Net primary productivity (NPP) represents the chemical energy converted from the energy of the sun by photosynthesis. The CASA model uses light use efficiency to estimate NPP at large or small scales. This paper expresses the supply of the CSS through the NPP value. NPP is represented by establishing a relationship between absorbed photosynthetically active radiation (APAR) and light use efficiency (Ɛ). The basic principle of the model is as follows50:

$$\:NPP\left(x,t\right)=APAR(x,t)\times\:\epsilon\:\left(x,t\right)$$
(1)

where NPP(x, t) denotes the NPP (gC·m−2·a−1) of image element x at time t; APAR(x, t) denotes the APAR(gC·m−2/month) absorbed by image element x at time t; and Ɛ(x, t) denotes the actual light energy utilization of image element x at time t. The algorithm of APAR(x, t) is51:

$$\:APAR\left(x,t\right)=SOL(x,t)\times\:FPAR(x,t)\times\:0.5$$
(2)

where FPAR(x, t) denotes the proportion of photosynthetically active radiation absorbed by the plant at time t by image x, which can be determined from the NDVI; SOL(x, t) denotes the total solar radiation at time t by image x; and 0.5 is the proportion of total solar radiation available to the vegetation. The algorithm for Ɛ(x, t) is as follows:

$$\:\epsilon\:\left(x,t\right)={T}_{\epsilon\:1}(x,t)\times\:{T}_{\epsilon\:2}(x,t)\:\times\:{W}_{\epsilon\:}(x,t)\times\:{\epsilon\:}_{max}$$
(3)

where TƐ1(x, t) and TƐ2(x, t) are temperature stress factors reflecting the effect of temperature on light energy use efficiency; WƐ(x, t) is a moisture stress factor indicating the impact of moisture on light energy use efficiency; and Ɛmaxis the maximum light energy use efficiency under ideal conditions. Based on existing studies52,53, this paper assigns different Ɛmax values to different site types, as shown in Table 1. Zhu et al.54 describe the model in more detail.

Table 1 Ɛmax for different site types.
Full size table

CSS demand model

This paper’s CE includes direct and indirect CEs55, representing the need for a CSS. Direct CE comes from farmland, forest, grassland, water, and unused land56. Farmland has a positive CE factor and is a net CE source, whereas other land use types play a role in CS57. Indirect CEs are emissions from urban energy production and living on self‒built land58. This paper selects coke, raw coal, gasoline, kerosene, diesel oil, natural gas, and power for calculation. The formula for calculating the CE is as follows:

$$\:E={E}_{D}+{E}_{I}$$
(4)

where ED is direct CE and EI is indirect CE. ED is calculated as follows:

$$\:{E}_{D}=\sum\:_{i=1}^{n}{A}_{i}\times\:{\epsilon\:}_{i}$$
(5)

where Ai denotes the area of land use type i and Ɛi denotes the CE coefficient (g·m−2·a−1)59,60,61 of land use type i, as shown in Table 2. The EI is calculated as follows:

$$\:{E}_{I}={C}_{e}\times\:{\alpha\:}_{e}\times\:{\beta\:}_{e}$$
(6)

where Ce denotes the total energy consumption of the lower e energy sources; αe denotes the conversion coefficient between the e energy source and standard coal, and βedenotes the energy CE coefficient. Referring to existing studies62,63 and the IPCC guidelines for calculating CEs (2006), the details are shown in Table 3.

Table 2 CE factors for land use types.
Full size table
Table 3 Energy standard coal conversion factor and CE factor.
Full size table

CSSF coupling process model

The flow of the CSS from the supply area to the demand area generally follows the principle of distance attenuation, which is affected by various natural and social factors, such as water flow, wind direction, land use type, population, and economy64. This paper uses the breaking point model and the metacoupling framework to measure the CSSF.

The model is divided into three steps. First, the radial distance from the supply area to the demand area is calculated via the following equation:

$$\:{d}_{ij}={D}_{ij}/(1+\sqrt{{m}_{j}/{m}_{i}})$$
(7)

where dij denotes the radial distance from supply zone i to demand zone j; Dij denotes the distance from supply zone i to demand zone j; and mi and mjdenote the flow difference in the supply zone and the flow difference in the demand zone, respectively21. The formulas for mi and mj are the same in this paper:

$$\:{m}_{i}=NPP\left(x,t\right)-E$$
(8)

Then, the field model is used to calculate the average radiation value between regions, and the formula is as follows65:

$$\:{F}_{ij}={m}_{i}/{D}_{ij}^{2}$$
(9)

where Fij denotes the average radiation value between regions.

Finally, the CSSF is calculated according to the average radiation value and radiation range, and the formula is as follows:

$$\:{E}_{ij}={F}_{ij}\times\:{A}_{i}\times\:\alpha\:$$
(10)

where Eij denotes the flow from supply area i to demand area j; Ai denotes the radial area of supply area i, which is a circular area geometrically centered on i; and α denotes the spatial conversion coefficient. α is generally thought to be 0.6, Ma et al.66, Wang et al.67 (2021) and Wang et al.68 (2024) set it as 0.6 when studying ecosystem service flow in Loess Plateau, Yangtze River Delta region and Guangdong region respectively. Therefore, α of CSSF in Chongqing is set as 0.6 in this paper.

To assess the CSSF from the districts and counties in the “one district‒two groups”, the area ratio method was employed to determine the supply area of each district and county. In the breakpoint model, we calculate the movement of the CSS from the supply area to other regions. Then, we use the supply area as the demand area to calculate the flow of services from other regions into this area.

To quantify the flow of the CSS between the “one district‒two groups” and the adjacent and distant systems, the “one district‒two groups” is considered the supply area. The breakpoint model calculates the CSSF for adjacent and distant systems. The natural breakpoint method was used to divide districts and counties into three levels: low outflow, medium outflow, and high outflow. The remaining CSS outflows are the main CSS outflows, with the exception of the low outflow districts and counties.

Then, considering the “one district‒two groups” as the demand area, the breakpoint model was used to calculate the CSSF from the districts and counties of the adjacent and distant system to the demand area. We divided these districts and counties into three levels: low inflow, medium inflow, and high inflow. The rest are the main CSS inflows, with the exception of the low‒inflow districts and counties. The major CSS outflow and inflow areas of the focal system were divided according to the same principle.

CSSF multiscenario model

CSSF is characterized by three attributes: flow direction, flow rate, and flow volume. Based on the characteristics of the spatial distribution of the CSS supply and demand in the study area, the flow path of the CSS in the study area was determined via the SPANs model. The logical relationships among the supply, flow, and demand of the CSS were established via the BBNs, and the results were spatialized via integration with ArcGIS software. We analyzed the area and evolution trend of different regions under the carbon peak (CP) and carbon neutrality (CN) scenarios.

SPANs model

The SPANs model is a spatial framework for quantifying ecosystem service flows on the basis of the “source‒sink‒use” idea69. This paper’s service flow parameters are screened according to the model framework and localized on the basis of the research content, as shown in Table 4.

Table 4 SPANs model conceptual framework localization.
Full size table

BBNs model

The BBNs model is a method for representing and reasoning about uncertain data and consists of two components: the directed acyclic graph and the conditional probability table (CPT)71. The nodes of the model represent the probability distributions of the variables, and the directed edges between the nodes represent the CPT‒based causal relationships. We combine the prior data with the actual data to find the conditional probability distribution for each node72. For a detailed description of the BBNs model, please refer to Pham73. This paper uses Netica32 (https://norsys.com/netica.html) and GeNIe 4.0 (https://www.bayesfusion.com/) software to model Chongqing’s CSSF. This includes identifying the areas where CS supply services originate (“source” areas), where CS supply services are utilized (“use” areas), and where the CSS is depleted (“sink” areas).

SPANs‒BBNs model parameterization

Based on the model objectives and the CSSF’s supply and demand results, the model contains both discrete and continuous data. To simplify the calculation, we used ArcGIS to discretize the data and selected a sample of 36,143 data points. Of these, 80% were used for model learning, and 20% were used for model testing. The model comprises 21 nodes (see Table B.1) and is illustrated in Fig. A.1 (Fig. A.1 in Supplementary Information).

SPANs‒BBNs model evaluation

In this paper, we assess the performance of SPANs‒BBNs from two perspectives. First, we conducted a sensitivity analysis using mutual information (MI) and variance of belief (VB) to gauge how well the target nodes align with the input nodes. The larger the MI and VB are, the stronger the node relationship becomes. Afterward, the accuracy of the model is evaluated based on the overall accuracy of the confusion matrix72. The sensitivity results for the main nodes of the model are presented in Table B.2. The model’s overall accuracy exceeds 85%, demonstrating its strong applicability and the high scientific validity and accuracy of CSSF prediction under various scenarios. The detailed confusion matrices can be found in Tables 5, 6 and 7.

Table 5 Confusion matrix of the CSS “sink” area.
Full size table
Table 6 Confusion matrix of the CSS “use” area.
Full size table
Table 7 Confusion matrix of the CSS “source” area.
Full size table

SPANs‒BBNs model multiscenario simulation

The paper uses the CSS of Chongqing in 2020 as a base scenario and then modifies key data points while keeping others constant to simulate different scenarios. Based on the sensitivity analysis, this paper focuses on changes in parameters such as key nodes A, APAR, POP, GDP, Coke, and Rawcoal to predict the CSS trends in Chongqing for 2030 under the CP and CN scenarios. According to the Chinese policy and existing research results, the specific setting parameters are shown in Table 8. According to the parameter settings of the key nodes under the two scenarios, the simulation results of the key nodes of the CSSF in Chongqing under different scenarios were obtained by substituting different conditional probabilities into the SPANs-BBNs model. According to the node data of the simulation results, the table data were converted into spatial distribution data by ArcMap software.

CP scenario assumes that carbon emissions peak just before 2030. During this period, economic growth decreases, the share of coal consumption in energy consumption decreases, population growth decreases, and the capacity of ecosystems to sequester carbon increases. Specific indicators are set as follows: The China Peak Carbon by 2030 Study predicts that China’s GDP will grow at an average annual rate of more than 5% during the Fourteenth and Fifteenth Five-Year Plan periods. According to the National Population Development Plan (2016–2030), the country’s total population will reach its peak in 2030, with a growth rate of 2.1%. According to the data description of CMIP6 by the Center for Resource and Environmental Science and Data of the Chinese Academy of Sciences, through the data description of SSP1-2.6 and SSP2-4.5, the temperature in these two scenarios increases by 36.8% compared with the 2020 data of this paper. The key nodes of this paper do not have temperatures, so we are realized through “APAR” and “A” (the principle of CN scenario is the same). The World and China Energy Outlook 2060 suggests that China’s share of coal consumption in 2025 will be 8% lower than in 2020, with the addition of a 15% reduction in coal consumption in 2030 compared to 2025. There is no coal in the key node of this paper, so it is realized by “Coke” and “Rawcoal” (the principle of CN scenario is the same).

CN scenario assumes that carbon neutrality will be achieved by 2060, in order to project 2030 data. During this period, the economy continues to develop at a high rate, energy consumption continues to decrease, the population grows moderately, and the CSS supply service is deployed more slowly, so the carbon sequestration capacity still shows an upward trend. The specific indicators are set as follows: According to Deng et al.74, CN can be achieved around 2060 under the SSP1-2.6 scenario, under which precipitation rises by 10.2% in the mean value of temperature and precipitation compared to the 2020 data of this paper. According to the setup in the Study on China’s Long-Term Low-Carbon Development Strategies and Transition Paths, China’s achievement of CN in 2060 can be equivalent to the situation studied in the 1.5 °C temperature control target scenario of the study, with a 40% coal share in 2030, which is a 16.8% reduction compared to 2020. Based on the historical data from the statistical yearbook, the population data will grow by 3.8% and GDP by 6.6% in 2030.

Table 8 Scenario parameter setting (unit: %).
Full size table

Results

Supply and demand of the CSS in CEMRFFLG assessment

CSS supply

From 2000 to 2020 (Fig. A.2‒A.5; Fig. A.2‒A.5 in Supplementary Information), there was little change in the spatial class distribution of NPP in the water, forest, and grassland subsystems, with significant differences in the farmland subsystem. The contribution rate of NPP value of water, forest, farmland, and grassland subsystems to the total NPP value was 0.54%, 44.83%, 36.78%, and 7.03%, respectively. The primary sources of NPP are forest and farmland.

From 2000 to 2020 (Fig. A.2), the water subsystem recorded relatively high NPP levels, which were located mainly along the YRB. Wanzhou, Fuling, Changshou, Hechuan, Jiangjin, and Yunyang had relatively high NPP values, averaging 1.88 × 1010gC/m2. Yuzhong, Dadukou, Shapingba, Qijiang, Nanchuan, and Chengkou had lower NPP values, averaging 5.05 × 109gC/m2.

From 2000 to 2020 (Fig. A.3), the forest subsystem recorded relatively high NPP levels, mainly in SC and NC. Chengkou, Wuxi, Fengjie, Youyang, Pengshui, and Fengjie presented relatively high NPP values, with a mean value of 8.24 × 1012gC/m2. Moreover, Yuzhong, Dadukou, Jiangbei, Nan’an, and Tongnan had lower NPP values, with a mean value of 6.84 × 1010gC/m2. Among them, Yuzhong had the lowest NPP value, at less than 6.40 × 107gC/m2.

From 2000 to 2020 (Fig. A.4), the farmland subsystem recorded higher NPP levels, mainly in the NC; the NPP value decreased by 37.53%. The areas with decreasing NPP values are mainly in CM and SC, possibly because urban expansion has led to reduced farmland. Hechuan, Kaizhou, and Wanzhou consistently had high NPP values, with levels above 5 × 1011gC/m2, whereas Yuzhong had the lowest value, averaging 3.29 × 107gC/m2.

From 2000 to 2020 (Fig. A.5), the grassland subsystem recorded relatively high NPP levels, mainly in the NC. Fengjie, Yunyang, Chengkou, Kaizhou, Wuxi, and Wanzhou presented relatively high NPP values, with a mean value of 2.82 × 1011gC/m2. Notably, Yuzhong, Dadukou, Jiangbei, and Nan’an had no grasslands, so they were not analyzed. Tongliang, Rongchang, Jiulongpo, Bishan, Shapingba, Beibei, and Yongchuan had lower NPP values, with a mean value of 4.69 × 108gC/m2.

CSS demand

From 2000 to 2020 (Fig. A.6‒A.9; Fig. A.6‒A.9 in Supplementary Information), the spatial distributions of the CEs of the four subsystems were the same, and the values changed little. The average CE values for the water, forest, farmland, and grassland subsystems are −1.51 × 1010gC/m2, −2.03 × 1012gC/m2, 1.79 × 1012gC/m2, and −1.66 × 1010gC/m2, respectively. Farmland is the primary source of CEs.

In the water subsystem (Fig. A.6), areas with greater carbon sequestration were concentrated in the SC. The lowest level of carbon sequestration was in Chengkou, with an average value of 2.25 × 107gC/m2, whereas the highest level of carbon sequestration was in Jiangjin, with an average value of 1.24 × 109gC/m2. From 2000 to 2020, Wanzhou, Fuling, Yunyang, Fengjie, Wushan, and Wanzhou significantly increased, with an increasing trend of approximately 112% in carbon sequestration.

In the forest subsystem (Fig. A.7), areas with greater carbon sequestration were concentrated in CM. From 2000 to 2020, the values fluctuated and were relatively small. The lowest level of carbon sequestration was in Yuzhong, with a mean value of 1.18 × 107gC/m2, whereas the highest level of carbon sequestration was in Youyang, with a mean value of 1.96 × 1011gC/m2.

In the farmland subsystem (Fig. A.8), the distribution of CEs was relatively uniform, and the values showed a general downward trend. The lowest CE is in Yuzhong, with an average value of 1.41 × 107gC/m2. Nanchuan and Kaizhou show clear increasing trends, with increases of 6% and 18%, respectively.

In the grassland subsystem (Fig. A.9), areas with greater carbon sequestration were concentrated in CM, and the values fluctuated. The Yuzhong, Dadukou, Jiangbei, and Nan’an regions do not have grasslands, so they were not included in the analysis. The amount of carbon sequestered in Kaizhou decreased significantly from 1.87 × 109gC/m2 in 2000 to 1.08 × 109gC/m2 in 2020, a decrease of approximately 42%. Carbon sequestration in Pengshui generally increased from 6.60 × 108gC/m2 in 2000 to 1.05 × 109gC/m2 in 2020, a rise of approximately 60%.

Metacoupling of the CSSF in CEMRFFLG

Intracoupling of the CSSF in CEMRFFLG

The rank distribution for intracoupling is determined by adding up the outflows and inflows.

From 2000 to 2020 (Figs. 3 and 4, Fig. A.10‒A.15; Fig. A.10‒A.15 in Supplementary Information), the total flow first increased but then decreased. In the outflow plot, the total flow in the supply area changes as 4.94 × 1013g→6.73 × 1013g→4.96 × 1013g. In the inflow plot, the total flow in the demand area changes as 1.54 × 1013g→2.10 × 1013g→1.53 × 1013g. The difference in the outflow plot is obvious in the spatial rank distribution. The rank distribution of the inflow plot is generally consistent.

In the water subsystem, the outflow plot (Fig. 3) shows that the total flow from the supply area has a mean value of 8.34 × 1011g. The maximum flow value in the supply area is in Hechuan, with a mean value of 3.82 × 1011g, accounting for 45.87% of the total supply area. From 2000 to 2020, the high‒level phases accounted for 29.41%, 47.06%, and 5.88%, respectively. In the inflow plot (Fig. 4), the total flow in the demand area has an average value of 3.36 × 1010g. The largest flow in the supply area is in Tongnan, with an average value of 1.04 × 1010g, representing 30.82% of the entire supply area. The high levels all accounted for approximately 30%.

In the forest subsystem, in the outflow plot (Fig. A.10), the total mean flow value provided by the supply area is 7.83 × 1012g. Among them, the maximum flow value in the supply area is Youyang, with a mean value of 4.78 × 1012g, which accounts for 61.06% of the total flow in the supply area. All high‒level outflows occurred in 2010, and no high‒level outflows occurred in 2020. In the inflow plot (Fig. A.11), the total flow in the demand area has an average value of 2.75 × 1012g. The smallest flow in the supply area is Qianjiang, with an average value of 5.74 × 1011g, accounting for 20.84% of the total flow in the supply area. All were high levels of inflows in 2010.

In the farmland subsystem, the outflow plot (Fig. A.12) shows that the total flow provided by the supply area has a mean value of 3.22 × 1013g. There were no low‒level outflows in 2010. In the inflow plot (Fig. A.13), the demand area has a total flow mean of 4.78 × 1012g; the smallest flow in the supply area is Yuzhong, with a mean of 2.97 × 104g, which is less than 1% of the total flow in the supply area. The high levels all accounted for approximately 24%.

In the grassland subsystem, in the outflow plot (Fig. A.14), the mean value of the total flow provided by the supply area was 5.92 × 1012g. Among the supply areas, Kaizhou has the smallest flow, with a mean value of 1.61 × 1012g, accounting for 27.17% of the total flow. From 2000 to 2020, the high‒level periods accounted for 25%, 50%, and 12.5%, respectively. In the inflow plot (Fig. A.15), the demand area has an average flow of 1.33 × 1012g. The highest flow in the supply area is Chengkou, with an average of 4.39 × 1011g, accounting for 33.06% of the total flow. The proportion of high levels is consistently approximately 50%.

Fig. 3
figure 3

Outflow plot of intracoupling of water subsystem. (a) Hierarchical distribution of CSSFs in 2000. (b) Hierarchical distribution of CSSFs in 2010. (c) Hierarchical distribution of CSSFs in 2020. (d) CSSF outflow values in 2000. (e) CSSF outflow values in 2010. (f) CSSF outflow values in 2020.

Full size image
Fig. 4
figure 4

Inflow plot of intracoupling of water subsystem. (a) Hierarchical distribution of CSSFs in 2000. (b) Hierarchical distribution of CSSFs in 2010. (c) Hierarchical distribution of CSSFs in 2020. (d) CSSF inflow values in 2000. (e) CSSF inflow values in 2010. (f) CSSF inflow values in 2020.

Full size image

Pericoupling of the CSSF in CEMRFFLG

The spatial distribution and total flow characteristics (Figs. 5 and 6, Fig. A.16‒A.21; Fig. A.16‒A.21 in Supplementary Information) during the coupling of the focal and adjacent systems are generally consistent with the intracoupling. The total flow in the outflow area varies from 4.35 × 1014g→5.81 × 1014g→4.07 × 1014g; the total flow in the inflow area varies from 2.13 × 1012g→3.16 × 1012g→2.45 × 1012g.

In the water subsystem (Fig. 5), the largest outflow area is Chengkou, with an average flow of 2.08 × 1011g, accounting for 7.34% of the total flow. The smallest outflow area is Yunyang, accounting for 4.42% of the total area. From 2000 to 2020, high‒level outflows accounted for 35.29%, 58.82%, and 5.89%, respectively. As shown in Fig. 6, the largest inflow area is Yunyang, with a flow increasing from 3.67 × 108g in 2000 to 3.90 × 109g in 2020, an increase of 408.66%, accounting for 26.75% of the total flow. Qianjiang, Liangping, Chengkou, and Youyang all have smaller flows, each accounting for no more than 1% of all flows. The high‒level distribution increased by 1.5 times.

In the forest subsystem (Fig. A.16), Yuzhong had the largest outflow, with an average value of 1.32 × 1013g, accounting for 4.41% of the total outflow. Wuxi had the smallest outflow, representing only 1.91% of the total. Only Wuxi had a low grade in 2010. As shown in Fig. A.17, the inflow into Wuxi was significantly greater than that into other areas, increasing from 3.31 × 1011g in 2000 to 4.21 × 1011g in 2020, a 33.56% increase. The flows in Yuzhong, Dadukou, Jiangbei, and Shapingba are relatively small, each accounting for less than 1% of the total.

In the farmland subsystem (Fig. A.18), the outflow in all areas showed a significant upward trend in 2010, with a high‒level distribution and an overall increase of approximately 30%. The differences in the flow between different districts and counties were relatively small. In the inflow area (Fig. A.19), Kaizhou had significantly higher flow than did the other areas, increasing from 5.43 × 1010g in 2000 to 1.03 × 1011g in 2020, an 89.96% increase. Xiushan had the smallest flow, with an average value of 6.49 × 109g, representing only 1.11% of the total inflow.

In the grassland subsystem, the absence of grassland was removed. In the outflow area (Fig. A.20), the lowest flow was in Qijiang, with an average value of 2.39 × 1012g, accounting for 3.03% of the total flow. No low level was recorded in 2010. Among the inflow areas (Fig. A.21), the largest flow was Shizhu, which decreased from 2.25 × 1010g in 2000 to 1.30 × 1010g in 2010, a decrease of 42.39%, accounting for 24.81% of the total flow. The flows in Jiulongpo, Shapingba, Tongliang, and Bishan were relatively small, each accounting for less than 1% of the total flow.

Fig. 5
figure 5

Outflow plot of pericoupling of water subsystem. (a) Hierarchical distribution of CSSFs in 2000. (b) Hierarchical distribution of CSSFs in 2010. (c) Hierarchical distribution of CSSFs in 2020. (d) CSSF outflow values from 2000 to 2020.

Full size image
Fig. 6
figure 6

Inflow plot of pericoupling of water subsystem. (a) Hierarchical distribution of CSSFs in 2000. (b) Hierarchical distribution of CSSFs in 2010. (c) Hierarchical distribution of CSSFs in 2020. (d) CSSF inflow values from 2000 to 2020.

Full size image

Telecoupling of the CSSF in CEMRFFLG

From 2000 to 2020 (Figs. 7 and 8, Fig. A.22‒A.27; Fig. A.22‒A.27 in Supplementary Information), during the coupling process between the focal and distant systems, the flow in the outflow area first increased but then decreased. The total flow rate in the outflow area changes to 1.67 × 1015g→2.40 × 1015g→1.51 × 1015g. In the inflow area, the flow rate decreases and then increases, with the total flow rate in the inflow area changing to 2.43 × 1014g→2.09 × 1014g→3.37 × 1014g. The differences in spatial class distributions are more noticeable in the outflow plot, with a significant increase in class distribution in 2010, and in the inflow plot, with a significant increase in class distribution in 2020.

In the water subsystem, the spatial class distribution is more uniform and less varied. In the outflow area (Fig. 7), most low‒level areas are distributed in the middle reaches of the YRB. Naqu and Hezhou have relatively high flows, both exceeding 2 × 1011g. Among the inflow areas (Fig. 8), Haixi, Jingzhou, and Yushu have relatively high flows, and the largest flow is in Jingzhou, with an average value of 2.92 × 1011g, accounting for 8.90% of the total flow. The main inflow areas are located in the middle reaches of the YRB, with high grades accounting for approximately 8% of the total area.

In the forest subsystem (Fig. A.22), the outflow area had a large proportion of low levels, averaging 31.08%, with a clear distribution of high levels in the lower reaches of the YRB. In 2010, the outflow area experienced a significant increase in flow, from 8.16 × 1014g in 2000 to 1.18 × 1015g in 2010, indicating a 44.84% increase. In the inflow area (Fig. A.23), the flow in 2020 also increased significantly, from 1.59 × 1014g in 2000 to 2.17 × 1014g in 2020, an increase of 36.73%. Among them, the regions of Ganzi, Liangshan, Aba, and Ganzhou experienced relatively high flows. Ganzi has the highest flow, with a change in flow of 1.38 × 1013g→1.29 × 1013g→2.17 × 1013g, reflecting an increase of approximately 57.60% compared with 2000 and accounting for 9.50% of the total flow.

In the farmland subsystem (Fig. A.24), in the outflow area, the total flow increased from 2.19 × 1013g in 2000 to 3.38 × 1013g in 2020, a 54.01% increase. There is no high‒level distribution in the lower reaches of the YRB. The flows in Nanyang, Liangshan, and Hanzhong are significantly greater than those in other areas, with the highest flow occurring in Nanyang, averaging 1.16 × 1012g, accounting for 4.85% of the total flow. As shown in Fig. A.25, the inflow area mainly consists of low‒level flows, averaging 33.08%. In 2010, only Nanyang and Liangshan had low‒level flows below 4 × 1012g.

In the grassland subsystem (Fig. A.26), the main outflow areas are located mainly in the middle and lower reaches of the YRB. In 2010, there was a significant increase in the outflow area, from 2.77 × 1014g in 2000 to 3.73 × 1014g in 2010, representing a 34.60% increase. As shown in Fig. A.27, the main inflow areas were located primarily in the upper YRB, with Ganzi, Aba, and Yushu having significantly higher inflows than other areas. Ganzi had the largest inflow, with an average value of 2.20 × 1013g, accounting for 33.30% of the total inflow.

Fig. 7
figure 7

Outflow plot of telecoupling of water subsystem. (a) Hierarchical distribution of CSSFs in 2000. (b) Hierarchical distribution of CSSFs in 2010. (c) Hierarchical distribution of CSSFs in 2020. (d) CSSF outflow values from 2000 to 2020.

Full size image
Fig. 8
figure 8

Inflow plot of telecoupling of water subsystem. (a) Hierarchical distribution of CSSFs in 2000. (b) Hierarchical distribution of CSSFs in 2010. (c) Hierarchical distribution of CSSFs in 2020. (d) CSSF inflow values from 2000 to 2020.

Full size image

Multiscenario results of CSSF in CEMRFFLG

Figure 9 shows that changes in the “sink” area are insignificant; the “source” area shows an increasing trend, with changes mainly concentrated in NC and SC; the “use” area shows a decreasing trend, with changes mainly concentrated in the CM. In the CP scenario, the “source” area of high‒level CSS increases by 188.72%, and the “use” area of high‒level CSS decreases by 63.83%. The supply and demand for the CSS are most moderate under the CP scenario.

In the three scenarios (Fig. 9a, b, and 9c), the higher “source” areas for CSS are mainly along the Daba Mountains, the Dalou Mountains, and the Wuling Mountains. The CSS “sink” areas are primarily determined by the soil respiration rate, with high “sink” areas mainly concentrated in forests and grasslands, resulting in minimal differences among the three scenarios. High “use” areas for the CSS are predominantly concentrated in urban districts and counties, indicating that regions with greater human activity and faster economic development are the main demand areas for the CSS. Overall, the flow of CSS primarily moves from areas with significant forest cover, such as Wuxi, Chengkou, and Fengjie, etc., to areas with more human activity, such as Dadukou, Jiulongpo, and Yuchong, etc.

Among the “source” areas (Fig. 9a, b, and 9c), under the base scenario, only Wuxi, Chengkou, Fengjie, Nanchuan, and Kaizhou have a high‒level CSS “source” area share of more than 10%. The “source” areas of low CSS account for more, with 32 areas, including Tongliang, Tongnan, and Dadukou, having more than 50% of the low level. Compared with the base scenario, under the CP and CN scenarios, the districts and counties with a high CSS “source” area account for more than 10% of the total area increase from 5 to 12 and 7, respectively. The regions with a greater reduction in the share of “source” areas served by low CSS include the Chengkou, Yunyang, and Wushan districts and counties, all of which have declined by more than 50%.

Among the “sink” areas (Fig. 9a, b, and 9c), in the base scenario, CP scenario, and CN scenario, the low CSS sink area is the main area, accounting for 86.84%, 76.32%, and 71.05%, respectively. The average proportion of regions with high‒level status is approximately 29%, and the proportion of high CSS “sink” is the largest in Wuxi, which is approximately 50% in all the scenarios.

Among the “use” areas (Fig. 9a, b, and 9c), in the base scenario, more than 50% of the high‒CSS areas are Dadukou, Jiulongpo, Jiangbei, etc. The percentage of areas with more than 50% low‒CSS areas is 44.74%. Compared with the base scenario, the percentages of areas with more than a 50% reduction in high‒CSS “use” areas under the CP and CN scenarios are 83.87% and 45.16%, respectively. The percentages of regions with more than a 50% increase in the “use” area for low CSS are 48.65% and 27.03%, respectively.

Fig. 9
figure 9

CSSF paths in Chongqing under different scenarios. (a) Base scenario. (b) CP scenario. (c) CN scenario. This map depicted in the figure was generated by the authors using ArcGIS 10.4 (Esri, URL: https://www.esri.com/en-us/arcgis/products/arcgis-desktop/overview) and does not require any license.

Full size image

Discussion

Global implications of the metacoupling framework

The metacoupling framework contains multiple systems of focal, adjacent, and distant systems, which not only increases the complexity between systems also allows to deal with complex features (e.g., problems such as nonlinearities75, time lags76, and legacy effects77), and has become one of the world’s most complex frameworks one of the most complex1. With the application of the metacoupling framework, scholars have conducted a large number of studies from the provincial scale78, national scale79to the global scale1, which involve a wide range of areas such as ecosystem services, goods, environmental footprints, virtual water, virtual carbon, and other fields80,81, enriching the existing theoretical system of environmental economy. Due to the omnidirectionality of CSSFs82, the CSSFs of the CEMRFFLG in Chongqing can affect the neighboring provinces and municipalities of Chongqing, China, and even have a certain impact on CSSFs on a global scale. Based on the differences in supply and demand of CSSs, the metacoupling framework allows for the balance of CSSFs within systems as well as across spatial distances, allowing for a more complete inter-system study, with significant potential in the integrated assessment of CSSFs. With rapid population and economic growth, it will increase the imbalance between the supply and demand of CSSs in the CEMRFFLG across space and time, which will change the strength of the metacoupling. In order to fully realize the potential of the framework, more time and resources are needed to conduct more systematic and integrated research on the CEMRFFLG to accelerate the further application of the framework worldwide and to achieve global sustainable development.

Differences in the focal‒adjacent‒distant coupling of the CSS in different subsystems

Based on the results of the focal–adjacent–distant coupling from 2000 to 2020, it is evident that the CSS outflow generally exhibits an inverse relationship with the CSS inflow. To further analyze this relationship, we discuss the CSS outflow in the outflow plot and the CSS inflow in the inflow plot for the four subsystems in 2020 (Fig. 10).

In intracoupling (Fig. 10a), an inverse linear relationship between CSS outflow and inflow is observed in the water and farmland subsystems. This relationship is weaker in the forest and grassland subsystems, being most notable in the forest subsystem. The main factors contributing to this are the distance between supply and demand areas and the level of CSS demand in the demand areas. When the distance difference is insignificant, areas with lower CSS demand experience greater inflows; when the difference in demand for CSSs is insignificant, more distant areas have larger inflows. Pengshui, Youyang, Kaizhou, Fengjie, and Yunyang in the outflow plot are supply areas. Given the lower CSS demand in the supply area, the outflow from it is less than that from the demand area, leading to a significant supply‒demand gap. Compared to the water and farmland subsystems, the supply‒demand distance in the grassland and forest subsystem is shorter, and the more conspicuous decrease in inflow causes a departure from the linear relationship.

In pericoupling (Fig. 10b), an inverse linear relationship exists between the outflows and inflows of CSSs in the four subsystems, with a weaker linear relationship in the forest subsystem. This is mainly due to the geomorphology of Chongqing. The eastern, southern, and southeastern parts of Chongqing are mostly mountain and hilly areas, such as the Wushan, Wuling, and Daloushan mountain ranges. The special geomorphology makes the mountainous areas of SC and NC rich in forest resources. Forest areas are more evenly distributed in SC and NC, narrowing the difference between supply and demand areas. The MC has a smaller distribution of forest areas, high demand, and low supply. As a result, the linear relationship between supply and demand areas is weakened. The distribution trend of resources in the supply and demand areas of other subsystems is the same, the resources in the supply area are more abundant, and the resources in the demand area are more scarce. Ample resources confer a higher supply capacity; when resources are relatively scarce, the demand for resources is high, which is more evident in the CSS inflows.

In telecoupling (Fig. 10c), the outflows and inflows of CSSs in the four subsystems exhibit a negative correlation, yet the linear relationship is attenuated, and in some regions, the outflows of CSSs are less than the inflows. Telecoupling encompasses larger ecosystems, and its flow mechanism is influenced by multiple factors. On one hand, it may result from internal transformation within the distant system, which diminishes the coupling effect between the focal and distant systems via the interaction of CSSs in the neighboring region. On the other hand, various factors such as climate, topography, land use, policies, and human activities in different regions need to be taken into account. Consequently, the linear relationship fails to precisely depict the relationship between inflows and outflows. Under metacoupling, overall CSS outflows exceed inflows, but in telecoupling, certain regions have CSS outflows smaller than inflows. These regions are mainly located in Yunnan, Sichuan, Shaanxi, Qinghai, Jiangxi, Hunan, and Hubei, primarily because these regions possess relatively more abundant resources compared to the “one district‒two groups”.

Fig. 10
figure 10

Correlation between the outflows and inflows of CSSs. (a) The situation of four subsystems in intracoupling. (b) The situation of four subsystems in pericoupling. (c) The situation of four subsystems in telecoupling.

Full size image

Differences between the CSSFs are analyzed based on their share of the four subsystems in each region from 2000 to 2020. Due to the regional inconsistency of the focal system in this paper and the small amount of data for intracoupling and pericoupling, only the regions in the distant system are analyzed here. According to Fig. 11, the region’s most significant change from 2000 to 2020 was a decrease in CSSFs of grassland, farmland, and water subsystems to compensate for the increase in CSSFs of the forest subsystem. The forest subsystem and the farmland, grassland and water subsystems present trade-offs (when an increase in the CSSFs from the forest subsystem leads to a decrease in the CSSFs from the rest of the subsystems83), and the farmland, grassland and water subsystems present synergy (simultaneous decreases in the CSSFs from the subsystems). Land use change and urbanization may lead to a decrease in CSSFs in the farmland and grassland subsystems; the decrease in CSSFs in the water subsystem is related to changes in hydrological conditions due to water management policies and climate change. According to the Guiding Opinions on Strengthening Afforestation and Greening of the Yangtze River Economic Belt issued in 2016, the YRB has vigorously carried out the return of farmland to forests and grasslands to improve forest quality; during the 13th Five-Year Plan period, it is clear that it will be necessary to increase the carbon sink of forests and to reduce the intensity of carbon dioxide emissions. Based on the above reasons, it is possible that the CSSF of the grassland, water, and farmland subsystems may be sacrificed in order to enhance the CSSF of the forest subsystem.

Fig. 11
figure 11

Transformation rates of the share of subsystem CSSFs in the distant system by region from 2000 to 2020.

Full size image

Optimal design of the CSSF in CEMRFFLG

This paper uses the self‒organizing mapping method84 to identify the carbon sink service clusters (CSSCs) in the districts and counties of Chongqing on the basis of the similarity of the CSSFs of the four subsystems in each district and county. These clusters are then mapped in space and assigned to the CSSC. As shown in Figs. 12 and 13, in terms of spatial distribution, the supply side of the SC is B1, and the demand side is Bi. SC is the ecological barrier area of the Wuling Mountains. This area has a high percentage of forest cover and favorable conditions for afforestation. The area of farmland, water, and grassland areas are all small, and the ecological environment quality is high85. The supply side of NC is dominated by the CSS in the grassland subsystem, with Chengkou, Wuxi, and Yunyang designated B5; Wanzhou, Kaizhou, Fengjie, and Wushan designated B4; the remaining areas designated B6; and the demand side is all Biv. The NC contains the ecological barrier area of the Dabashan Mountains, with limited forested areas and the largest proportion of grassland. Farmland accounts for approximately 38% of the city area and is an essential main production area for agricultural products in Chongqing. Water resources come from the main stream and many tributaries of the YRB, but the region still faces a shortage of water resources due to karstic landforms86. The land use structure of the CM is quite complex, with the farmland subsystem being the dominant feature. Hechuan is classified as B2 and requires forest services. Tongnan, Changshou, Fuling, Ba’nan, and Jiangjin are classified as B3 and require forest services, whereas the remaining areas are classified as B6 and need more green services. The farmland area of the CM accounts for about 47% of the city, and the area of grassland and forest land is relatively scarce87. It contains runoff from many rivers such as the main stream of the YRB, Jialing River, and Wuriver, and is rich in water resources.

This paper proposes optimization suggestions on the basis of the spatial distribution pattern of supply and demand service clusters related to the CSSF within the CEMRFFLG in Chongqing. These suggestions consider natural ecological features and relevant policies.

(1) SC should increase the forest area and expand farmland and water bodies appropriately. SC’s forest ecosystem provides ample space, making it a key area for building the “Green Mountains on Both Sides of the Taiwan Strait–Thousand Miles of Forest Belt” and the National Reserve Forest. This expansion aligns with national policy goals and helps address rocky desertification issues in SC, thereby increasing the CSS of the forest ecosystem and compensating for the CSSF in CM. SC water and farmland ecosystems have relatively low CSSs, so measures should be taken to build lakes and wetlands in suitable areas to increase their CSS while protecting the original ecological environment.

(2) NC should protect blue and green spaces and reasonably increase the area of farmland. Protecting the blue–green space also helps maintain the CSS of the forest, grassland, and water subsystems and compensates for the CSS flows in the CM. Interregional water transfer, including the Sichuan‒Chongqing‒Northeast Integrated Water Resources Allocation Project and the Changzheng Canal Diversion Project, has also increased water resources. The Daba Mountain ecological barrier area and the Wushan‒Qiyao Mountain ecological barrier area in NC contain critical CSSs in the forest and grassland subsystems of Chengkou, Wuxi, and Yunyang, followed by Wanzhou, Kaizhou, Fengjie, and Wushan. The Liangping, Dianjiang, Zhongxian, and Fengdu areas have the largest farmlands, but their CSSs are insufficient. Therefore, expanding the amount of arable land is necessary to enhance its CSS.

(3) CM should respect the red line for farmland and establish compensation mechanisms for the CSSF. The CM is dominated by farmland, which is an essential ecological conservation area. It is appropriate to “return farmland to forests and grasslands” to protect the farmland red line and increase the CSS of the forest and grassland subsystems. Hechuan is a major grain‒producing county with three river confluences; Tongnan, Changshou, Fuling, Banan, and Jiangjin are the next most important areas. The remaining regions provide fewer CSSs in the farmland subsystem, the land is tight, and the area available for afforestation is relatively small. Relying on the ecological pattern of “three belts, four screens, multiple corridors, and multiple points”, a compensation mechanism for the CSSF can be constructed.

Fig. 12
figure 12

Supply‒side CSSC in the CSSF of the EC in Chongqing. (a) Spatial-temporal patterns of CSSC. (b) Composition and relative magnitude of CSSFs in CSSC. (B1: forest service cluster; B2: high food‒water service cluster; B3: low food‒water service cluster; B4: low green service cluster; B5: high green service cluster; B6: low food service cluster).

Full size image
Fig. 13
figure 13

Demand‒side CSSC in the CSSF of the EC in Chongqing. (a) Spatial-temporal patterns of CSSC. (b) Composition and relative magnitude of CSSFs in CSSC. (Bi: high food‒water‒grass synergistic services demand cluster (SDC); Bii: high ecological SDC; Biii: low ecological SDC; Biv: food‒water‒forest synergistic SDC; Bv: low green SDC; Bvi: forest‒field‒water synergistic SDC).

Full size image

Deficiencies and prospects

This paper offers a reference for Chongqing to enhance understanding of the relationships among CSS supply, flow, and demand within the CEMRFFLG. It lays a foundation for ecosystem management and restoration, yet has certain limitations. Firstly, the application of different land use data in various models was overlooked. Here, CSS supply was measured with a single model for expediency, with the CASA model used for forest CS assessment and all land use data based thereon, potentially increasing data error. Secondly, the interaction between adjacent and distant systems was not explored. In the focal‒adjacent‒distant coupling process, only the interactions between the focal system and adjacent/distant systems were considered, while the CSSF and mutual influence between adjacent and distant systems were neglected. Thirdly, insufficient multiscenario modeling was conducted. Three scenarios were explored based on national policies, with multiscenario parameters set subjectively, necessitating further exploration to combine subjective and objective data for simulation.

Given these shortcomings, future research is warranted. Firstly, different LUCC data and more suitable models should be considered. The CASA model may not be able to accurately assess the characteristics of the ecosystems of the water, farmland and grassland subsystems, but it has become the norm for most studies to calculate the CSs of each land use type through this model. As different models are different from each other due to different principles and calculation methods, they also bring errors between models. Which specific way brings less error needs to be further explored in the future. Secondly, the connection between adjacent and distant systems should be incorporated. Future work could improve this connection within the metacoupling framework. Thirdly, a combination of shared socioeconomic pathways (SSPs) and representative concentration paths (RCPs) should be employed for multiscenario setup. The SSP focuses on global social, demographic, and economic futures without considering climate change impacts, while the RCP accounts for greenhouse gas and radiative forcing effects. Future research could involve further discussions in the context of SSPs‒RCPs.

Conclusions

(1) From 2000 to 2020, the average NPP values for the water, forest, farmland, and grassland subsystems were 0.54%, 44.83%, 36.78%, and 7.03%, respectively. From 2000 to 2020, the distribution of CEs across the four subsystems remained consistent, with only minor changes in values.

(2) The focal system of this paper is “one district‒two groups”, with Chongqing as the adjacent system and the YRB as the distant system. During intracoupling, the farmland subsystem has the highest total value in the outflow and inflow plots, accounting for 1.77% and 2.64% of the total, respectively. For pericoupling, the forest subsystem has the highest total value in the outflow and inflow plots, accounting for 12.50% and 0.52% of the total, respectively. Finally, in telecoupling, the forest subsystem again has the highest total value in the outflow and inflow plots, accounting for 35.51% and 61.24% of the total, respectively.

(3) In this paper, the “source” area refers to the location of the supply of the CSS, the “sink” area refers to the potential demand area and potential supply area of the CSSF, and the “use” area is where the demand for the CSS is located. The base scenario used in the paper is the CSS in Chongqing in 2020, with only the data of key nodes being changed to create the CP and CN scenarios. Under the CP scenario, the “source” area for high CSS increases by 188.72%, the changes in the “sink” area are not very noticeable, whereas the “use” area for high CSS decreases by 63.83%. The overall trend observed was that the CSSF mainly originated from areas with large forested areas, such as Wuxi, Chengkou, and Fengjie, to areas with more human activities, such as Dadukou, Jiulongpo, and Yuchong.

(4) In focal–adjacent–distant coupling, the amount of outflow flow from the supply area of the outflow plot is related to the supply of resources in the area itself, and the more resources there are the higher the outflow. The amount of inflow flow in the demand area of the outflow plot is related to the proximity of the supply area and the difference in the flow of the demand area of its CSS. Generally, the inflow plot’s flow is inversely proportional to the outflow plot’s flow. The forest subsystem and the farmland, grassland and water subsystems present trade-offs, and the farmland, grassland and water subsystems present synergy.

(5) To optimize the spatial distribution of the CSS in Chongqing’s CEMRFFLG, we put forward corresponding suggestions based on the ecological pattern of “one district‒two groups”.