Spatiotemporal Evolution and Multi-Scenario Prediction of Carbon Storage in the GBA Based on PLUS–InVEST Models

Abstract

In the context of sustainable development and dual-carbon construction, in order to clarify the future changes in land use and carbon storage in the GBA, this study used the PLUS and InVEST models as well as Geoda software to simulate and predict the spatial development pattern of land use as well as the changes in carbon stocks in the GBA in 2030 under multiple scenarios. The results show that (1) From 1990 to 2020, carbon storage decreased year by year. (2) In 2030, except for the EPS, the future carbon storage prediction values of the remaining scenarios are lower than those in 2020, especially the carbon storage prediction value under the EDS, which is the lowest at 8.65 × 10⁸ t. (3) The spatial distribution of carbon storage in the GBA has significant spatial heterogeneity. The high-value areas of carbon storage are distributed in the east and west wings as well as southwest of the GBA, while the low-value areas are concentrated in the middle and east. The research results can provide a reasonable scientific basis for the territorial space resource planning of the GBA under the goal of “dual carbon”.

1. Introduction

The terrestrial ecosystem serves as a substantial carbon sink on Earth, and its carbon storage plays a significant role in the global carbon cycle and climate change [1,2,3]. Land use/cover change (LUCC) is a concrete embodiment of the interaction between human society, the economy, and nature. The carbon sequestration capacity of different land use types is also different [4,5]. For example, changes in land use brought on by urbanization and deforestation will directly reduce the amount of carbon stored in local terrestrial ecosystems [6]. The Intergovernmental Panel on Climate Change (IPCC) report points out that the impact of land use change on the carbon cycle of terrestrial ecosystems has become a major part of the field of climate change, and it is the second largest source of carbon emissions after the burning of fossil fuels [7]. Both quantitative analysis and prediction of the impact of regional land use change on terrestrial ecosystem carbon stocks are critical, as human activities are regulated and driven by policy [8,9].

Nowadays, the quantitative methods of terrestrial ecosystem carbon storage mainly include field investigation, model analysis, and remote sensing interpretation [10,11,12]. Among them, field surveys can struggle to quantify carbon storage at large scales, and remote sensing interpretations are weak in terms of mechanism interpretation. CA-Markov, CLUE-S, GeoSOS-FLUS, and PLUS combined with the InVEST model have been widely used to predict land use and carbon storage changes at different spatiotemporal scales due to their strong explanatory power and superiority at spatiotemporal scales [13,14,15,16]. At present, many scholars have carried out research on space–time evolution and the influencing factors of ecosystem carbon storage from different scopes and angles. For example, Zhu, Gong, and Tian et al. [17,18,19] analyzed carbon stocks in different regions based on a combination of multiple models. However, most studies have focused on analyzing the impact mechanisms of land use and ecosystem carbon stocks, while few studies have investigated spatial zoning management of ecosystem carbon stocks in urban agglomerations driven by policy factors. Urban agglomeration is an important part of the terrestrial ecosystem, and its land use change is mainly driven by government policies [20,21]. Especially in recent years, the ecological policies implemented by the Chinese government in urban agglomerations in China have affected the carbon sequestration capacity of terrestrial ecosystems. Therefore, in the context of ecological policies implemented by the Chinese government, it is of great significance to clarify the change of carbon storage in urban agglomerations in China for the formulation of future carbon emission policies in China.

With the rapid development of global urbanization, the urban area continues to expand, the population explodes, and the ecological environment further deteriorates [20,21]. As one of the four major bay areas in the world, the Guangdong–Hong Kong–Macao Greater Bay Area (GBA) is not only a crucial component of the economic development pattern, but also a typical research area for the growth of urban agglomeration [22]. Especially in recent years, under the guidance of China’s “dual carbon” goal, the GBA has needed to not only overcome the problems of energy structure adjustment and carbon emission reduction while developing the city economy, but also to solve the balance between urban development and ecological environment [23]. However, while current studies on the spatial and temporal characteristics of concentrated carbon density in the GBA are generally available, studies on the spatial zoning management of carbon stocks are relatively lacking, especially in future prediction [19,24]. Therefore, the objective of this paper is to clarify the temporal and spatial characteristics and future trends of carbon storage in the GBA against a background of ecological policies, and also to provide data for the future carbon emission policies of urban agglomerations in China.

Based on land use change data from 1990 to 2020, coupled with the PLUS–InVEST model, this study predicts and analyzes changes in carbon stocks in the GBA under different development scenarios, and uses the spatial autocorrelation analysis method to study the zoning management of carbon stocks in the GBA. Additionally, we explore the effects of land use change on ecosystem carbon stocks under multiple scenarios in the future, revealing the response mechanism of land use to carbon stocks in various scenarios. Finally, the research results of this paper suggest that the GBA’s carbon storage shows an annual trend of shrinking and that the decline is steadily slowing down. This indicates that more attention should be paid to urban greening and land management in future urban development. The research results of this paper can provide reference significance for adjusting the policies of relevant departments, as well as a reasonable scientific basis for the territorial space resource planning of the GBA under the goal of “dual carbon”.

2. Materials and Methods

2.1. Study Area

The GBA is located on the southeast coast of China (111°21′ E–115°25′ E, 21°34′ N–24°23′ N), with a total land area of 5.6 × 10⁴ km² (Figure 1). It has a humid subtropical monsoon climate, with an average annual temperature of 21.4–22.4 °C and an annual rainfall of 1300–2500 mm. The GBA’s economy is growing rapidly, with its GDP accounting for about 12% of China’s total economic development. The GBA is not only the economic hinterland of South China but also the central link of the domestic economic cycle and an important hub of the international cycle. However, with the accelerating urbanization process of the GBA and the increasing demand for development and construction land, the contradiction between social and economic development and environmental and ecological protection continues to intensify [25].

2.2. Data Resources

The land use data are from the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (http://www.resdc.cn/ (accessed on 18 January 2023)) for the years 1990, 2000, 2010, and 2020, with a resolution of 30 m. According to the research requirements, 25 secondary classifications of land use data were reclassified and combined into 6 primary classifications, such as cropland, forestland, grassland, water area, construction land, and unused land.

In this paper, according to the existing research results and the specific situation of the study area, 11 impact factors were selected from the natural factors and social factors (

Table 1. Data source and description.

Data Type	Data Name	Data Accuracy	Source
Natural factor	Digital elevation model	30 m	The United States Geological Survey (https://earthexplorer.usgs.gov/ (accessed on 10 November 2022))
	Slope	30 m
	Mean annual temperature	1 km	Resources and Environmental Sciences and Data Center, Chinese Academy of Sciences (http://www.resdc.cn/ (accessed on 12 November 2022))
	Mean annual precipitation	1 km
	Soil type	1 km
	Distance to open water	30 m	National Center for Basic Geographic Information (http://www.ngcc.cn/ngcc (accessed on 15 November 2022))
	Distance to river	30 m
Social factor	Population	1 km	Resources and Environmental Sciences and Data Center, Chinese Academy of Sciences (http://www.resdc.cn/ (accessed on 13 November 2022))
	Gross domestic product	1 km
	Distance to national highway	30 m	National Center for Basic Geographic Information (http://www.ngcc.cn/ngcc (accessed on 15 November 2022))
	Distance to highway	30 m
	Distance to railway	30 m
	Distance to government	30 m

) [26,27]. In addition, the data were uniformly preprocessed based on the PLUS model’s input requirements (including clipping, projection transformation, and resampling, among other things). The spatial resolution was 30 m raster data, and the geographic coordinate system was GCS_WGS_1984.

2.3. Research Framework

The research framework in this paper consists of three parts: the PLUS model for simulating land use/cover (LUCC) data, the InVEST model for estimating ecosystem carbon stocks, and the Geoda model for spatial autocorrelation analysis (Figure 2). The specific process is as follows: Firstly, based on the GBA land use data from 1990 to 2020, this paper uses the PLUS model to simulate land use change in 2030 under four different urban development scenarios (such as the Natural Development Scenario (NDS), Economic Development Scenario (EDS), Cropland Protection Scenario (CPS), and Ecological Protection Scenario (EPS)). Secondly, according to the GBA climate conditions, the carbon density in this region is modified by the carbon density model, and the carbon storage changes under four scenarios from 1990 to 2020 and 2030 are simulated by the InVEST model. Finally, the Geoda model is used to analyze the correlation of the GBA carbon storage spatial distribution.

2.4. Methods

2.4.1. Multi-Scenario Setting

Multiple factors influence future urban development and land use change; thus, when simulating and forecasting land use change, various environmental factors must be fully considered. In accordance with the optimization policy of “Urban, Agricultural, and Ecological Spatial Layout of Guangdong Province” and related literature [28,29], this paper combined the actual situation of the GBA with the formulation of conversion rules and conversion rates between increasing and decreasing land classes to set up four scenarios (

Table 2. Scenarios and principles of future land use change.

Scenario Type	Principles and Restricted Development Zones
NDS	Based on the 2010–2020 GBA land use conversion law simulation (the binding effects of planning policies on land use change are not considered).
EDS	Based on the needs of the rapid economic and urban development of the GBA, while following the laws of nature (the conversion probability of cropland, forestland, water area, and unused land to construction land will be increased by 20%, respectively, and construction land will not be transferred to other land types).
CPS	Strictly implement cropland protection tasks based on the GBA Development Program Outline (reduce the conversion rate of cropland to construction land by 60% and set the area where cropland is located as a restricted development zone).
EPS	Based on the “green development and ecological protection” requirements in the GBA Development Program Outline (the conversion probability of forestland and grassland to construction land will be reduced by 50%, and the conversion probability of cropland and water area to construction land will be reduced by 30%. At the same time, the reduced part is included in the probability of the transfer of cropland to forestland, and the forestland and water areas are set as restricted development zones).

).

In addition, in order to reflect the change of land use type in a certain period of time, we adopted the dynamic attitude method of land use for calculation. According to the research object, it can be divided into single land use type dynamic attitude (K) and comprehensive land use dynamic attitude (LC) [11,19]. The calculation method is as follows:

$$K=\frac{U_b-U_a}{U_a}\times \frac{1}{T}\times 100\%$$

(1)

where U_a and U_b are the area of this land use type at the beginning and end of the study, respectively, (km²). T is the interval between the end of the study and the beginning of the study (year).

$$LC=\left(\frac{{{\displaystyle \sum }}_{i=1}^n\Delta L{U}_{i-j}}{2{{\displaystyle \sum }}_{i=1}^{n}L{U}_i}\right)\times \frac{1}{T}\times 100\%$$

(2)

where LU_i is the area of the i land use type in this region at the beginning of the study period (km²). ΔLU_i−j is the absolute logarithm of the area where the i land type changed into other land types from the beginning to the end of the study (km²). n is the total number of land types in this area. T is the interval between the end of the study and the beginning of the study (year).

2.4.2. PLUS Model

The patch-level Land Use Simulation model uses the rules of land expansion analysis strategy to mine the Land Expansion Analysis Strategy (LEAS) framework and the CA model based on multi-type random seeds (CARS). In this way, the landscape pattern of different scenarios in the future can be simulated with higher precision [30,31].

Firstly, the accuracy of the model needs to be verified. The land use data of the 2000 and 2010 periods are loaded into the extract land expansion module of the PLUS model to obtain the 2000–2010 land use expansion map. Secondly, the land use expansion map and various influence factors are loaded into the LEAS module of the model to obtain the development probability of each region during this period, and the land use of each region in 2020 is simulated by Markov chains. Finally, the land use data of 2010 and the development potential of each region are loaded into the CARS module, and various parameters are input. The simulation of land use in 2020 was compared with the actual land use data in 2020 to verify its accuracy. The results showed that the kappa coefficient was 0.89, the Fom was 0.553, and the Overall Accuracy (OA) was 0.92, with a good simulation effect and high reliability.

Therefore, the specific parameter settings of this study are determined as follows: In the LEAS module, the number of regression trees is set to 50, and mTry is set to 15, indicating the number of features used to train the random forest regression model. The sampling rate is set to 0.09 by default, meaning that about 9% of the pixels are selected for training. In the CARS module, the neighborhood weight is obtained by calculating the expansion law in the historical process of each land use type. The domain factor can directly reflect the strength of the expansion ability of each category, and the value range is 0–1. The closer a value is to 0, the weaker the expansion capacity of this type of land; thus, the more easily the land type is transformed into other land types. The total area of each land use type patch can quantitatively describe the degree of expansion [32]. In this paper, the land use data from 1990 to 2020 are used to calculate the area of each land type by transfer matrix, and the dimensionless processing is combined with Formula (1) to determine the neighborhood weight. Finally, based on the neighborhood weight under the natural development scenario, combined with the relevant literature and after many tests [33], the neighborhood weights of each land type in the remaining scenarios are adjusted. The specific results are shown in

Table 3. Domain weight parameters.

Scenario Use	Cropland	Forestland	Grassland	Water Area	Construction Land	Unused Land
NDS	0	0.12	0.43	0.32	1	0.37
EDS	0	0.12	0.43	0.32	1	0.27
CPS	0.2	0.22	0.53	0.42	0.8	0.37
EPS	0.2	0.12	0.43	0.32	0.8	0.37

. The specific calculation formula is:

$$W_i=\frac{T{A}_i-T{A}_{min}}{T{A}_{max}-T{A}_{min}}$$

(3)

where W_i is the domain weight of land type i, TA_i is the expansion area of land use of type i, TA_min is the minimum expansion area of each type of land use, and TA_max is the maximum expansion area of each type of land use.

2.4.3. InVEST Model

The InVEST model’s carbon storage module categorizes ecosystem carbon storage into four categories: aboveground carbon storage, underground carbon storage, soil organic carbon storage, and dead organic matter carbon storage. Following land use classification, the average carbon density of four basic types of different land classes was calculated to obtain the total carbon storage in the study area [34]. The specific calculation formula is:

$$C_{total}=\left(C_{above}+C_{below}+C_{soil}+C_{dead}\right)\times A_i$$

(4)

where i is the average carbon density of the earth class, A_i is the area of the earth class, C_toal is the total carbon storage of all land types (t/hm²), C_above is the aboveground carbon storage (t/hm²), C_below is the underground carbon storage (t/hm²), C_soil is the soil organic carbon storage (t/hm²), and C_dead is the dead organic matter carbon storage (t/hm²).

2.4.4. Carbon Density Correction

The carbon density data come from the National Data Center for Ecological Sciences (http://www.cnern.org.cn/ (accessed on 20 January 2023)). According to the climatic conditions of the GBA and related research results, the biomass carbon density and soil organic carbon density were modified by temperature and precipitation factors [35,36]. The modified formula refers to the research results of Alam et al., Giardina et al., and Chen et al. [37,38]. The specific calculation is:

$$C_{SP}=3.968\times MAP+3996.1\left(R^2=0.11\right)$$

(5)

$$C_{BP}=6.798\times e^{0.0054MAP}\left(R^2=0.70\right)$$

(6)

$$C_{BT}=28\times MAT+398(R^2=0.47,P<0.01)$$

(7)

where C_SP is soil carbon density based on annual precipitation (kg/m²), C_BP and C_BT are biological carbon densities obtained according to annual precipitation and average annual temperature, respectively, (kg/m²), MAP is the average annual precipitation (mm), and MAT is the average annual temperature (°C). In this paper, the average annual temperature and precipitation of Guangdong Province and the GBA region are substituted into the above formula, respectively, (the average annual temperature and precipitation in the two regions were 22.1 °C and 22.54 °C, and 1799.2 mm and 1910.5 mm, respectively), and the ratio of the two is the correction factor. The carbon density data of Guangdong Province multiplied by the correction factor are used as the GBA’s carbon density data.

$$K_{BP}=\frac{C_{BPGBA}}{C_{BPGDP}}$$

(8)

$$K_{BT}=\frac{C_{BTGBA}}{C_{BTGDP}}$$

(9)

$$K_B=K_{BT}\times K_{BP}$$

(10)

$$K_S=\frac{C_{SPGBA}}{C_{SPGDP}}$$

(11)

where C_BPGBA and C_BPGDP are the GBA and Guangdong Province’s biomass carbon density, C_SPGBA and C_SPGDP are the GBA and Guangdong Province’s soil carbon density, K_BP and K_BT are the correction coefficients of the precipitation factor and temperature factor of biological carbon density, respectively, and K_B, K_S are the correction coefficients of biomass carbon density and soil carbon density, respectively. Revised carbon density values for different land use types for the GBA are shown in

Table 4. Carbon density values of various land use types (unit: t/hm²).

Land Use Type	Aboveground Carbon Storage	Belowground Carbon Storage	Soil Organic Carbon Storage	Dead Organic Matter Carbon Storage
Cropland	26.77	1.85	82.47	1
Forestland	42.09	15.32	164.32	6.5
Grassland	2.03	9.6	65.15	1.9
Water area	0.04	0	0	0
Construction land	0.22	1.72	59.13	0
Unused land	0.37	0	60.17	0

.

2.4.5. Spatial Autocorrelation Analysis

Spatial autocorrelation is a common method to test whether the attribute values of a feature are spatially related, as well as the degree of spatial correlation. It is used to measure the degree of aggregation or dispersion between attributes. It can be divided into Tocal Indicators of Spatial Association (TISA) and Local Indicators of Spatial Association (LISA) [39,40]. The calculation formula of TISA is:

$$MI=\frac{{{\displaystyle \sum }}_{i=1}^n{{\displaystyle \sum }}_{j=1}^n{w}_{ij}\left(x_i-\overline{x}\right)\left(x_j-\overline{x}\right)}{S^2{{\displaystyle \sum }}_{i=1}^n{{\displaystyle \sum }}_{j=1}^n{w}_{ij}}$$

(12)

n represents the number of subjects; x_i and x_j are the values of the i and j observation objects in space, respectively; w_ij is the spatial weight matrix; and the value of MI is between −1 and 1. A positive value means that the research object has a positive spatial correlation. The larger the value is, the closer the relationship between the units is; a negative value indicates that the distribution of the research objects has a negative spatial correlation and is discrete. If the value is 0, it means that there is no spatial autocorrelation in the distribution of research objects, which is a random distribution.

Since global spatial autocorrelation only reflects the overall distribution characteristics of spatial elements, the GBA is divided into a 1000 m × 1000 m grid and adopts local spatial autocorrelation (LISA) to show the spatial concentration of carbon reserves in local areas. The formula for LISA is as follows:

$$M{I}_i=\frac{\left(x_i-\overline{x}\right)}{S^2}{\displaystyle \sum }_i{w}_{ij}\left(x_j-\overline{x}\right)$$

(13)

The value of MI is between −1 and 1, and its variable meaning is the same as that of Formula (11). In this paper, the local concentration of carbon reserves is divided into five types: High-High Cluster, Low-Low Cluster, High-Low Cluster, Low-High Cluster, and Not Significant.

3. Results and Analysis

3.1. Analysis of Land Use Change and Simulation Analysis of Different Scenarios

3.1.1. Analysis of Land Use Change

From 1990 to 2020, the primary land type in the GBA will be forestland, which will account for 54% of the total area. Cropland accounts for 22% of the total area; the construction land area continues to grow, accounting for 15% of the total area by 2020. Grassland, water, and unused land account for less than 10% of the total area (

Table 5. Area of each year from 1990 to 2020 (unit/km²).

Land Use Type	1990	2000	2010	2020
Cropland	15,921.8	14,445.5	12,643.5951	12,096.8613
Forestland	30,920.5	30,647.4	30,062.1636	29,706.6879
Grassland	1273.36	1225.1	1101.8556	1187.3151
Water area	3852.64	4374.28	4095.2277	4035.3561
Construction land	3137.65	4466.62	7381.5984	8342.5905
Unused land	23.8311	23.607	11.4246	6.9876

). In the past 30 years, only the area of construction land continued to expand, but the increasing scale showed a decreasing trend. From 2000 to 2010, the comprehensive and dynamic attitude toward land use was at its highest. A large amount of cropland was transferred to construction land, resulting in a rapid increase of 2842.81 km² (

Table 6. Changes in the area of each land type from 1990 to 2020.

Time Scale	Property	Cropland	Forestland	Grassland	Water Area	Construction Land	Unused Land
1990–2000	Variable-area (km2)	−1478.49	−282.09	−49.91	517.89	1292.82	−0.22
	Dynamic Index (%)	−0.929	−0.091	−0.392	1.344	4.121	−0.094
	Comprehensive dynamic Index (%)	0.181
2000–2010	Variable-area (km2)	−1807.50	−587.00	−124.11	−404.35	2842.81	−12.24
	Dynamic Index (%)	−1.251	−0.192	−1.013	−0.713	6.365	−5.186
	Comprehensive dynamic Index (%)	0.515
2010–2020	Variable-area (km2)	−546.41	−370.22	83.96	−75.70	912.90	−4.53
	Dynamic Index (%)	−0.432	−0.123	0.763	−0.185	1.238	−3.974
	Comprehensive dynamic Index (%)	0.180

). The area of cropland, forestland, and unused land decreased continuously, among which the cropland area decreased the most, and cropland was mainly converted into construction land and water area. The main conversion of cropland to construction is mainly due to the rapid development of the economy and urbanization in the Greater Bay Area, resulting in a large amount of construction land occupying cropland. The main reason for the conversion of cropland into waters is mainly due to the implementation of the policy of “returning farmland to lakes and wet“ and the revival of the mulberry-based fish pond industry. The fluctuation range of grassland and water area is small, with grassland exhibiting “decrease–decrease–increase” characteristics and water area exhibiting “increase–decrease–decrease” characteristics (Figure 3).

3.1.2. Simulation Analysis of Different Scenarios

Under the four scenarios, the primary land types for the GBA in 2030 are forestland and cropland (Figure 4). Under the NDS, the GBA’s cropland, forestland, water area, and unused land are projected to decrease by 28.8249 km², 366.864 km², 81.5814 km², and 1.728 km², respectively, compared with 2020. The construction land and grassland area will continue to increase, with an increase of 79.0818 km² and 399.9164 km², respectively, of which the grassland has the fastest growth rate, of 0.67% (

Table 7. Area changes of different land types in different scenarios from 2020 to 2030.

Scenario Use	Property	Cropland	Forestland	Grassland	Water Area	Construction Land	Unused Land
NDS	Area (km2)	−28.8249	−366.864	79.0818	−81.5814	399.9164	−1.728
	Dynamic Index (%)	−0.024	−0.123	0.666	−0.202	0.479	−2.473
EDS	Area (km2)	−547.3449	−434.2275	71.2854	−131.8968	1042.3701	−0.1863
	Dynamic Index (%)	−0.452	−0.146	0.600	−0.327	1.249	−0.267
CPS	Area (km2)	13.9244	−255.3212	79.8911	−73.1946	235.2106	−0.513
	Dynamic Index (%)	0.012	−0.086	0.673	−0.181	0.282	−0.734
EPS	Area (km2)	38.2622	53.6386	−27.5361	0	−63.734	−0.63
	Dynamic Index (%)	0.032	0.018	−0.232	0	−0.076	−0.902

). From the perspective of transfer direction (Figure 5), cropland and construction land were mainly converted to grassland, and the transfer areas were 28.82 km² and 26.26 km², respectively. The forestland was converted to grassland and construction land, and the transfer area was 66.14 km² and 331.02 km², respectively. The water area was mainly converted into construction land, with a transfer area of 64.45 km². The total amount of unused land is small, and the area reduction is relatively small. It is mainly transferred to construction land, with an area of 1.54 km². In summary, although both construction land and grassland have been transferred out, the increase in construction land is the largest among all land types, mainly in forest land, so it shows an increasing trend. The unused land is the most intense change, mainly because only the land type is transferred out but no other land type is transferred in (Figure 5).

Under the EDS, the GBA’s cropland, forestland, water area, and unused land are projected to decrease by 547.3449 km², 434.2275 km², 131.8968 km², and 0.1863 km², respectively, compared with 2020. The grassland and construction land showed an increasing trend, with increments of 71.2854 km² and 1042.3701 km², respectively, among which construction land had the strongest change, with an increase of 1.25% (Table 7). From the perspective of transfer direction (Figure 4 and Figure 5), the cropland, forestland, water area, and unused land of the GBA are mainly converted into construction land, and the transfer area is 547.34 km², 346.16 km², 134.32 km², and 0.19 km², respectively. It can be seen that the development of construction land against the backdrop of rapid economic development will destroy the balance between land types and cause natural environmental problems.

Under the CPS, forestland, water area, and unused land in the GBA are projected to decrease by 255.3212 km², 73.1946 km², and 0.513 km², respectively, compared with 2020. Cropland, grassland, and construction land showed an increasing trend, with increments of 13.9244 km², 79.8911 km², and 235.2106 km², respectively, among which grassland had the highest growth rate, with an increase of 0.67% (Table 7). From the perspective of transfer direction (Figure 4 and Figure 5), forestland, water area, and unused land are mainly converted into construction land, but the transfer area is relatively small. The increase in cropland area was mainly due to the transfer of grassland and construction land, of which the transfer of construction land was greater, and the transfer area was 11.02 km².

Under the EPS, the cropland and forestland of the GBA are projected to increase by 38.2622 km² and 53.6386 km², respectively, over 2020. Grassland, construction land, and unused land showed a decreasing trend, decreasing by 27.5361 km², 63.734 km², and 0.63 km², respectively, of which unused land decreased to 0.9% (Table 7). From the perspective of transfer direction (Figure 4 and Figure 5), the area of construction land converted into cropland is 38 km², while the area of grassland and construction land converted into forestland is 35.33 km² and 18.1 km², respectively. The transfer of the remaining land types was not obvious.

3.2. Analysis of Carbon Stock Change and Simulation Analysis of Different Scenarios

In this paper, the InVEST model is used to analyze the temporal and spatial changes of carbon stocks in the GBA from 1990 to 2020, and the temporal and spatial changes of carbon stocks in the GBA in 2030 are predicted under different scenarios (Figure 6 and Figure 7). In terms of time change, the GBA’s carbon reserves from 1990 to 2020 were 9.14 × 10⁸ t, 8.98 × 10⁸ t, 8.82 × 10⁸ t, and 8.74 × 10⁸ t, with a gradually decreasing trend. In the past 30 years, carbon storage has been reduced by 41.3 × 10⁶ t, an average annual reduction of 1.38 × 10⁶ t. Among them, the carbon loss from 2000 to 2010 was the largest, at 0.17 × 10⁸ t. The GBA’s economic and urbanization growth accelerated during this period, which necessitated the development of a large amount of land for construction. After 2010, although construction land continued to expand, the expansion speed was relatively slow, land changes tended to be moderate, and the carbon loss rate was gradually alleviated.

In terms of the results of spatial variation, the spatial distribution of GBA carbon stocks has significant spatial heterogeneity, and the overall distribution is increasing from the middle to the periphery (Figure 6). The GBA’s carbon reserves are concentrated primarily in the northwest, southwest, and northeast regions, and the highest carbon storage density was 20.5407 t/hm². These regions have a strong carbon sequestration capacity due to topographic factors that limit the sprawl of construction land, and human activities cause relatively minor damage to the natural environment. The regions with low carbon reserves are distributed in the central area of the GBA. Due to the high population density, rapid urbanization, and land use type (mainly construction land), the carbon sequestration capacity of this region is weak.

Based on the analysis results of the GBA carbon storage in 2030 under four scenarios (Figure 7), under the NDS, the carbon storage is 8.68 × 10⁸ t, which is 5.64 × 10⁶ t less than that in 2020, and the average annual reduction is 0.56 × 10⁶ t. Under the EDS, the carbon storage is 8.65 × 10⁸ t, which is 9.1 × 10⁶ t less than that in 2020, and the average annual reduction is 0.91 × 10⁶ t. Under the CPS, the carbon storage is 8.70 × 10⁸ t, 3.6 × 10⁶ t less than that in 2020, and the average annual reduction is 0.36 × 10⁶ t. Under the EPS, the carbon storage is 8.75 × 10⁸ t, 1 × 10⁶ t more than that in 2020, and the average annual increase is 1 × 10⁵ t. The results show that the GBA’s carbon sequestration effect can be achieved by ecological protection measures, and the effect is significant.

3.3. Autocorrelation Analysis of Carbon Storage Space

Geoda software was used to conduct a TISA analysis on the spatial and temporal differences in carbon storage in all regions of the GBA in 2030 under four future scenarios. The results show that (Figure 8) in the NDS, Moran’s I value is 0.552. In the EDS, Moran’s I value is 0.555. In the CPS, Moran’s I value is 0.548. Moran’s I value in the EPS is 0.551. It means that there is a spatial clustering correlation rather than a random distribution condition.

The spatial correlation characteristics of the GBA are as follows: The regions with higher carbon reserves tend to be adjacent to regions with higher carbon reserves. The regions with lower carbon stocks tend to be adjacent to regions with lower carbon stocks. Since most of the regions are located in the first quadrant (hot spot area) and the third quadrant (cold spot area), the hot spot area belongs to high and high aggregation, while the cold spot area belongs to low and low aggregation.

ArcGIS was used to conduct LISA analyses on the carbon storage in the GBA under various scenarios (Figure 9). The results showed that the carbon storage in the four scenarios had similar spatial distributions. The regions with a high carbon storage value were concentrated in Zhaoqing, the south of Jiangmen, the north of Guangzhou City, and the north of Huizhou. The regions with low carbon storage values were concentrated in Shenzhen, the central part of Guangzhou, Zhongshan, and Zhuhai. In different scenarios, the high-value area gathered the most under the EPS, and the low-value area gathered the most under the EDS.

4. Discussion

4.1. Effects of Land Use Change on Carbon Stocks in the GBA

Results from GBA carbon stocks under different scenarios in 2030 (Figure 10) show that, under the NDS, GBA carbon storage is projected to be 8.68 × 10⁸ t, which is 5.64 × 10⁶ t less than that in 2020. It is caused by the transfer of cropland, forestland, and water areas to construction land. Under the EDS, a large number of forestlands are transferred to construction land, which leads to a further reduction of the predicted carbon storage value. The predicted carbon storage value was 8.65 × 10⁸ t. Under the CPS, the predicted value of carbon storage was 8.70 × 10⁸ t. Although cropland is no longer transferred, part of the forestland is still transferred, so the carbon storage will continue to decrease by 3.6 × 10⁶ t compared with 2020. Under the EPS, the predicted value of carbon storage was 8.75 × 10⁸ t. Carbon storage showed an increasing trend, increasing by 1 × 10⁵ t compared with 2020. It is mainly caused by restricting the transfer of cropland, forestland, and water areas, while some construction land is converted to forestland.

Under the NDS and EDS, based on the above conclusions, except for grassland and construction land, the significant reduction of ecological areas such as cropland, forestland, and water areas is a critical reason for reducing ecosystem carbon storage. Under the CPS, reducing forestland and water areas has a great impact on carbon storage. Under the EPS, the increase in cropland and forestland area leads to an increase in total ecosystem carbon storage. The main reason for this is that this scenario limits the expansion of construction land, and reduces the conversion of forestland, cropland, grassland, and water to construction land. Under the four scenarios, forestland and water area have the greatest impact on ecosystem carbon storage. Unused land has little effect on changes in carbon stocks [41,42]. Thus, the change in the GBA’s carbon stocks is mainly due to the expansion of construction land, which encroaches on other ecological lands. Although the future construction land expansion rate decreased, its progress is still increasing. Without effective environmental protection measures, there will still be carbon loss in parts of the GBA in the future.

4.2. Partition Management of Carbon Stocks in the GBA

This paper elucidates the spatial pattern of carbon stocks by dividing the spatial change values of carbon stocks in the GBA during 1990–2020 into three categories: decrease, unchanged, and increase. The results indicate that the spatial variation of GBA carbon stocks presents the characteristics of patch aggregation and sporadic distribution (Figure 11). The carbon storage of most GBA regions is unchanged and of the same type, and the area accounts for more than 70%. A quarter of the regional carbon storage showed a decreasing trend, mainly distributed in the central and eastern regions, namely Shenzhen, central Guangzhou, Zhongshan, and Zhuhai. Mainly due to the expansion of the GBA cities, there is an overall integration and a closer relationship with the surrounding areas. The construction land between the cities and the junction around the waters continues to expand, and the cropland in the marginal area is converted into construction land. The patch aggregation degree is increased and the carbon storage is reduced. The carbon storage that was in the central part of the GBA as well as sporadically distributed in coastal areas, that is, in southern Guangzhou and western Foshan, showed an increasing trend mainly due to the establishment of a forest ecological network system, of which improved the protection of forests and wetlands and park green space construction, enhancing its carbon sequestration capacity.

Based on the analysis of this paper and related studies, the promotion of ecological policies can protect forestland and grasslands to a certain extent, which will help to improve regional carbon storage, especially in the economically developed urban agglomeration [43,44,45]. As an economically developed city cluster, with the continuous progress of urbanization carbon storage in the central and eastern parts of the GBA will continue to decline in the future. Therefore, the carbon reserves of different regions of the GBA need to be managed in specific zones. Appropriate policy guidance should be established for land planning in different regions to limit the conversion of grassland and forestland with high carbon densities to construction land with low carbon densities, which will help improve the carbon reserves of the GBA in the future.

5. Conclusions

Based on land use change data from 1990 to 2020, this study, coupled with the PLUS–InVEST model, analyzed changes in carbon stocks of the GBA under different development scenarios. Spatial autocorrelation analysis was conducted by Geoda to study the zoning management of carbon storage. The following three conclusions can be drawn:

(1)

From 1990 to 2020, the GBA’s carbon storage showed a trend of decreasing year by year. From 2000 to 2010, the GBA’s carbon storage changed dramatically, and carbon loss was at its most severe. However, after 2010, the speed of carbon loss slowed down, indicating that the conversion rate of high carbon density land to low density land slowed down.

(2)

The prediction results for 2030 show that, except for the EPS, the future carbon storage values of all scenarios are lower than those of 2020. Especially in the EDS, the GBA has the lowest predicted value of carbon storage, which is only 8.65 × 10⁸ t, which is 9.1 × 10⁶ t less than that in 2020. This indicates that urbanization will be the critical factor affecting GBA carbon storage in the future.

(3)

There is a strong positive spatial correlation between GBA’s carbon storage and spatial distribution. Regions with higher carbon storage tend to be adjacent to regions with higher carbon storage, while regions with lower carbon storage tend to be adjacent to regions with lower carbon storage. Moreover, the carbon storage values of the four scenarios have certain similarities in their spatial distribution. The regions with high carbon storage values were distributed in the eastern, western, and southwestern parts of the GBA. The regions with low carbon storage values are clustered in the central and eastern regions. This suggests that Chinese government policies play a major role in the spatial distribution of carbon stocks.