Study area
With an area of about 2.33 × 105 km2, the Western Sichuan Plateau (27.11°–34.31°N and 97.36°–104.62°E) is located in the transition zone between the Qinghai-Tibet Plateau and the Sichuan Basin, including all of Garze Prefecture and Aba Prefecture, and parts of Liangshan Yi Autonomous Prefecture28 (Fig. 1). With an altitude of 780–7556 m, this area is dominated by mountain and ravine areas and high mountain and plateau areas, and the terrain is high in the west and low in the east. The climate belongs to the subtropical plateau monsoon climate, with large temperature difference between day and night and abundant sunshine. The annual average temperature is about 9.01–10.5 °C, and the precipitation is about 556.8–730 mm28. The study area is rich in water resources, including the Yalong River, Minjiang River and other important river systems in the upper reaches of the Yangtze River, and the Baihe River, Heihe river and other river systems of the Yellow River. The main types of soil are plateau meadow soil dark brown soil, brown soil, cold frozen soil and cinnamon soil, and the main vegetation types are alpine meadow and scrub. With rich and diverse soil vegetation types and distinctive vertical zonal distribution characteristics, it is one of the global biodiversity conservation hotspots29.
Location of the study area. The map is created in the support of ArcGIS 10.2 (ESRI). The China map and Western Sichuan Plateau boundary data were collected from Resources and Environmental Science and Data Center (http://www.resdc.cn/). The Qinghai-Tibetan Plateau boundary data were collected from the Global Change Research Data Publishing & Repository (http://www.geodoi.ac.cn/WebCn/Default.aspx).
Data source and processing
Multisource archival data were used in this study (Table 1). The land use remote sensing monitoring data, administrative boundary data and geological disaster vector data were obtained from Resources and Environmental Science and Data Center. The spatial resolution of land use remote sensing monitoring data is 30 × 30 m, including 6 first-level classification and 26s-level classification. The first-level classification includes cropland, woodland, grassland, water body, built-up land, and unused land. The accuracy of remote sensing classification is not less than 95% for cropland and built-up land, not less than 90% for grassland, woodland, and water body, and not less than 85% for unused land, which meets the need of the research. Landsat remote sensing monitoring data is used as the main information resources, among which Landsat-TM/ETM remote sensing monitoring data is used in 2000, 2005, 2010 and Landsat 8 remote sensing monitoring data is used in 2015 and 2020. In light of actual conditions and the implementation of policies and philosophies including the natural forest protection project, return of farmland to forest, land remediation, ecological civilization, the period from 2000 to 2020 is selected as the study period, and the land use data of each period is cropped using ArcGIS 10.2 to reclassify the 26 secondary classifications into cropland, woodland, grassland, water body, built-up land and unused land.
The DEM data were obtained from SRTM (Shuttle Radar Topography Mission) of Resources and Environmental Science and Data Center, the spatial resolution of 30 × 30 m, absolute horizontal accuracy ± 20 m, absolute elevation accuracy ± 16 m, elevation and slope are extracted from the downloaded DEM. The Qinghai-Tibetan Plateau boundary data were collected from the Global Change Research Data Publishing & Repository. Data of carbon density of different land types were obtained from Chinese Ecosystem Research Network Data Center (http://www.nesdc.org.cn/).
A total of 29,284 evaluation units were collected for spatial grid processing of the Western Sichuan Plateau according to 3 km × 3 km by ArcGIS 10.2. The impact factors obtained in this study include grid data per kilometer of GDP spatial distribution, grid data per kilometer of population spatial distribution, annual mean temperature spatial interpolation data, annual mean rainfall spatial interpolation data, long-term normalized difference vegetation index (NDVI) comes from Resources and Environmental Science and Data Center with a resolution of 1 km × 1 km. The Human Active Index (HAI), with a resolution of 30 m × 30 m, can be calculated by formula30,31, and the factors are discretized into the data type required for the geodetector by the natural breakpoint method.
Methods
The InVEST model
The InVEST model was developed by Stanford University, the University of Minnesota, the Nature Conservancy and the World Wide Fund for Nature (WWF). The model’s terrestrial ecosystem services assessment includes four modules: soil conservation, water retention, carbon storage and biodiversity assessment, and provides an overall measurement of regional ecosystem services32. The carbon storage model of the InVEST model divides the carbon storage of the ecosystem into 4 basic carbon pools, namely above-ground carbon, underground carbon, soil carbon, dead organic matter carbon7.
The calculation formula of total carbon storage in the Western Sichuan Plateau is as follows7:
$$C_{total} = C_{above} + C_{below} + C_{soil} + C_{dead}$$
(1)
In formula (1), Ctotal is the total carbon storage; Cabove is the above-ground carbon storage; Cbelow is the underground carbon storage; Csoil is the soil carbon storage, and Cdead is the dead organic matter carbon storage.
Based on the carbon density and land use data of different land use type, the carbon storage of each land use type in the Western Sichuan Plateau is calculated by the formula7:
$$C_{{text{total}}i} = (C_{{text{above}}i} + C_{{text{below}}i} + C_{{text{soil}}i} + C_{{text{dead}}i}) times A_{i}$$
(2)
In formula (2), i is the average carbon density of each land use, and Ai is the area of this land used.
The carbon density data of different land use types in this study were obtained from the shared date of the National Ecological Science Data Center and some documents33,34,35,36,37. Since the carbon density data were collected from the results of studies in different parts of China, the selected documents should be close to or similar to the study area as far as possible to avoid excessive data gap. At the same time, the carbon density varies with climate, soil properties and land use38, so the carbon density should be modified according to the climate characteristics and land use types of the Western Sichuan Plateau. Existing research results show that the carbon density is positively correlated with annual precipitation and weakly correlated with annual average temperature. The quantitative expression of the relationship between carbon density and temperature and precipitation is as follows39,40,41,42:
$$C_{SP} = 3.3968 times P + 3996.1;;left( {{text{R}}^{{2}} = 0.{11}} right)$$
(3)
$$C_{BP} = 6.7981e^{0.00541p};;;left( {{text{R}}^{{2}} = 0.{7}0} right)$$
(4)
$$C_{BT} = 28 times {text{T}} + 398;;left( {{text{R}}^{{2}} = 0.{47,};{text{P}} < 0.0{1}} right)$$
(5)
In these formula, CSP is the soil carbon density (kg m−2) based on the annual precipitation; CBP is the biomass carbon density (kg m−2) based on the annual precipitation; CBT is the biomass carbon density (kg m−2) based on annual average temperature; P is the average annual precipitation (mm), and T is the annual average temperature (°C). According to the data of China Meteorological Data Service Centre (http://data.cma.cn/), in the past 20 years, the average annual temperature of China and the Western Sichuan Plateau was 9.0 °C and 6.3 °C, and the average annual precipitation was 643.50 mm and 812.65 mm respectively.
The modified formula of carbon density in the Western Sichuan Plateau is as follows7:
$$K_{BP} = frac{C^{prime}{_{BP}}}{{C^{primeprime}{_{BP}}}}$$
(6)
$$K_{BT} = frac{C^{prime}{_{BT}}}{{C^{primeprime}{_{BT}}}}$$
(7)
$$C_{BT} = 28 times T + 398;;left( {{text{R}}^{{2}} = 0.{47,};{text{P}} < 0.0{1}} right)$$
(8)
$$K_{S} = frac{C^{prime}{_{SP}}}{{C^{primeprime}{_{SP}}}}$$
(9)
In these formula, KBP is the modified indices of precipitation factor in biomass carbon density; KBT is the modified indices of temperature factor; C’BP and C”BP are the biomass carbon density obtained from annual precipitation in the Western Sichuan Plateau and the whole country respectively. C’BT and C”BT are the biomass carbon density obtained from annual average temperature; C’SP and C”SP are the soil carbon density data obtained from annual average temperature; KB and KS are the biomass carbon density modified indices and soil carbon density modified indices respectively. The carbon density values of each land use type after modified in the Western Sichuan Plateau are shown in Table 2.
Exploratory spatial analysis method
Global spatial autocorrelation
Global Moran’s I was used to describe the spatial differentiation characteristics of carbon storage in the study area, and the expression formula is as follows43:
$$I = frac{{nsumnolimits_{i = 1}^{n} {sumnolimits_{j = 1}^{n} {w_{i,j} left( {x_{i} – overline{x} } right)left( {x_{j} – overline{x} } right)} } }}{{sumnolimits_{i = 1}^{n} {sumnolimits_{j = 1}^{n} {omega_{ij} } } sumnolimits_{i = 1}^{n} {left( {x_{i} – overline{x} } right)^{2} } }}$$
(10)
wij is the spatial weight; x is the attribute mean; xi and xj are the attribute values of elements i, j, respectively; n is the number of cells, and the correlation is considered significant when |Z| > 1.96.
Local indications of spatial association (LISA)
LISA reveals the local cluster characteristics of spatial unit attributes by analyzing the difference and significance between spatial units and surrounding units, and the expression formula is as follows42:
$$I_{i} (d) = frac{{n(x_{i} – overline{x} )sumnolimits_{j = 1}^{n} {w_{ij} (x_{j} – overline{x} )} }}{{sumnolimits_{i = 1}^{n} {(x_{j} – overline{x} )^{2} } }}$$
(11)
Correlation analysis
In order to evaluate the influence of natural factors and socioeconomic factors on carbon storage in the study area, the correlation coefficients of temperature, rainfall, NDVI, GDP, population density (PD), HAI and carbon storage were calculated according to the Pearson correlation coefficient method. The calculation formula is as follows44:
$$r_{xy} = frac{{sumnolimits_{i = 1}^{n} {(M_{i} – overline{x} )(y_{i} – overline{y} )} }}{{sqrt {sumnolimits_{i = 1}^{n} {(M_{i} – overline{x} )^{2} sumnolimits_{i = 1}^{n} {(y_{i} – overline{y} )} } } }}$$
(12)
rxy represents the correlation coefficient between x and y; Mi represents the carbon storage in the ith year; yi represents the value of the impact factor Y in the ith year, and ({overline{text{x}}}) and ({overline{text{y}}}) respectively represents the average value of carbon storage and impact factor in the research period over several years.
Human influence index analysis method
Land use is significantly spatially clustered in the study area31, and LUCC changes will have a certain impact on the structure and process of the ecosystem. HAI has the characteristics of spatial variability, which can reflect the impact of human activities on land use and landscape composition changes. In this study, Human Influence Index Analysis Method (HAI) index was used to analyze the correlation between carbon storage and human interference intensity in the Western Sichuan Plateau. The calculation formula is as follows30,
$$HAI = sumlimits_{i = 1}^{n} {left( {A_{i} P_{i} /TA} right)}$$
(13)
HAI is Human Active Index; Ai is the total area of the ith land use type; Pi The intensity parameter of human impact reflected by type i land use type; TA is the total final surface area of land use type in evaluation unit; n is the number of land use types. Combined with the land use type of this study, Pi is assigned by Delphi method, in which cropland is 0.67, woodland is 0.13, grassland is 0.12, water body is 0.10, built-up land is 0.96, and unused land is 0.0530,45.
Geodetector
Geodetector is an algorithm that uses spatial heterogeneity principle to detect driving factors of carbon storage, which can quantitatively detect the influence of impact factors on carbon storage and explore the interaction between driving factors. Geodetector includes factor detection, risk detection, interaction detection and ecological detection46.
Differentiation and factor detection: the influence factors were discretized, and then the significance test of the difference in the mean values of the impact factors was conducted to detect the relative importance among the factors. The statistical quantity q is used to measure the explanatory power of impact factors on the carbon storage spatial differentiation and the value range of q is between 0 and 1. The larger the value, the stronger the explanatory power of the factor47.
$$q = 1{ – }frac{{sumnolimits_{h = 1}^{L} {N_{h} sigma_{h}^{2} } }}{{Nsigma^{2} }}$$
(14)
In this formula, h = 1, 2…, L is the classification or partition of variable (Y) or factor (X); Nh and N are layer h and regional number units respectively; and (sigma_{h}^{2}) and (sigma_{{}}^{2}) are the variance of the layer h and regional value Y respectively.
The variance of the regional value Y is calculated as follows,
$$sigma^{2} = frac{{sumnolimits_{i = 1}^{n} {(Y_{i} – overline{Y} )^{2} } }}{N – 1}$$
(15)
where, Yi and (overline{Y}) are the mean value of sample j and the region Y, respectively.
$$sigma^{2} = frac{{sumnolimits_{i = 1}^{{n_{h} }} {(Y_{h,i} – overline{{Y_{h} }} )^{2} } }}{{N_{h} – 1}}$$
(16)
where, Y and (overline{Y}) are the value and mean value of sample i in layer h, respectively.
Interaction detection: it is used to identify the interaction between different impact factors Xs, that is, to evaluate whether the combined action of X1 and X2 will increase or weaken the explanatory power of vegetation coverage Y, or the influence of these factors on Y is independent of each other. The evaluation method is to first calculate the value q of the two factors X1 and X2 for Y respectively: q(X1) and q(X2), and calculate the value q of their interaction (the new polygon distribution formed by the tangent of the two layers of the superimposed variables X1 and X2) : q(X1 ∩ X2) and compare q(X1) and q(X2) with q(X1 ∩ X2)46.
Source: Ecology - nature.com
