Early Warning of the Carbon-Neutral Pressure Caused by Urban Agglomeration Growth: Evidence from an Urban Network-Based Cellular Automata Model in the Greater Bay Area

Abstract

Carbon neutrality is becoming an important development goal for regions and countries around the world. Land-use cover/change (LUCC), especially urban growth, as a major source of carbon emissions, has been extensively studied to support carbon-neutral planning. However, studies have typically used methods of small-scale urban growth simulation to model urban agglomeration growth to assist in carbon-neutral planning, ignoring the significant characteristics of the process to achieve carbon neutrality: large-scale and long-term. This paper proposes a framework to model large-scale and long-term urban growth, which couples a quantity module and a spatial module to model the quantity and spatial allocation of urban land, respectively. This framework integrates the inertia of historical land-use change, the driving effects of the urbanization law (S-curve), and the traction of the urban agglomeration network to model the long-term quantity change of urban land. Moreover, it couples a partitioned modeling framework, spatially heterogeneous rules derived by geographically weighted regression (GWR), and quantified land-use planning orientations to build a cellular automata (CA) model to accurately allocate the urbanized cells in a large-scale spatial domain. Taking the Guangdong–Hong Kong–Macao Greater Bay Area (GHMGBA) as an example, the proposed framework is calibrated by the urban growth from 2000 to 2010 and validated by that from 2010 to 2020. The figure of merit (FoM) of the results simulated by the framework is 0.2926, and the simulated results are also assessed by some evidence, which both confirm the good performance of the framework to model large-scale and long-term urban growth. Coupling with the coefficients proposed by the Intergovernmental Panel on Climate Change (IPCC), this framework is used to project the carbon emissions caused by urban growth in the GHMGBA from 2020 to 2050. The results indicate that Guangzhou, Foshan, Huizhou, and Jiangmen are under great pressure to achieve the carbon-neutral targets in the future, while Hong Kong, Macao, Shenzhen, and Zhuhai are relatively easy to bring up to the standard. This research contributes to the ability of land-use models to simulate large-scale and long-term urban growth to predict carbon emissions and to support the carbon-neutral planning of the GHMGBA.

1. Introduction

In recent decades, global climate change has attracted wide attention from academic communities, the Intergovernmental Panel on Climate Change (IPCC), and countries [1,2]. In particular, the aggravation of the greenhouse effect has made climate protection an imperative topic for sustainable development [3,4]. Carbon emissions are a major contributor to the global greenhouse effect, primarily caused by fossil fuel consumption in urbanization and industrialization processes [5,6,7]. Many countries have formulated their carbon-neutral plans to realize the sustainable development of human society, which involves many aspects, including clean energy development, energy conservation, emission reduction, tree planting, and land-use planning [8,9,10,11,12,13]. Among them, land use is an important source of carbon emissions, and its planning directly affects the level of regional urbanization and industrialization development, thereby affecting the implementation of carbon-neutral plans [14,15,16].

Urban growth is the most active process in land-use/cover change (LUCC), and the carbon emissions associated with it are the primary source of land-use carbon emissions [17,18]. Therefore, most land-use planning focuses on the quantity and spatial allocation of urban land. This demand, combined with the wide use of 3S (Remote Sensing RS, Geographical Information Systems GIS, and Global Positioning System GPS) technologies, has led to the development of many data-driven approaches to model urban growth and support land-use planning, such as the agent-based model (ABM), cellular automata (CA), and their derivatives [19,20,21,22,23,24,25,26]. In particular, CA is the most widely used model to simulate urban growth because of its simple and open structure and advantages in easy combination with machine learning methods [27,28]. Numerous scholars have devoted themselves to improving the simulation capability of CA by optimizing its basic components that directly affect simulation performance, such as transition rules and neighborhood [29,30,31,32,33]. The neighborhood is the basic component of CA to represent local interactions between land patches and simulate the micro-self-organization process in urban growth [34]. Through optimizing the structure and local rules of the neighborhood, several kinds of neighborhoods have been designed to improve the simulation performance, including the distance-decay neighborhood, the spatial heterogeneity-weighted neighborhood, the dual-scale neighborhood, and so on [35,36,37]. Transition rules are at the core of CA and determine the quantity and spatial allocation of urban land in the simulation process [38]. The derivation of transition rules has undergone several stages, from subjective to objective and from external phenomenon analysis to internal mechanism exploration, which continually improves the simulation capability of CA [32,39]. Nowadays, machine learning methods such as the artificial neural network (ANN), random forest (RF), support vector machine (SVM), and genetic algorithm (GA) are commonly used in building transition rules, increasing the ability of CA to characterize complex urban dynamics significantly [40,41,42,43].

However, carbon neutrality is the long-term goal of large-scale regions, such as river basins and countries [44,45], while the widely used CA framework is usually developed for simulating urban growth in small-scale regions, like cities [27,28,46]. This means that even though the simulation capability of CA has been dramatically improved with the incorporation of machine learning, it still needs to be improved to accurately simulate large-scale urban growth and effectively support carbon-neutral plans [31,47,48,49]. The difference between small-scale and large-scale urban growth simulation is that large-scale urban growth simulation needs to consider the network structure and complex connections between cities [50]. Furthermore, the differences between urban growth characteristics and mechanisms in different cities will be more obvious in large-scale regions, which will significantly affect the simulation performance of urban growth. This means that the spatial heterogeneity of urban growth mechanisms should be incorporated into the derivation of CA transition rules for large-scale urban growth simulation [51,52]. The differences mentioned above can be represented in CA modeling as the differences between top-down macro quantity-control rules and bottom-up micro spatial-layout rules for small-scale and large-scale regions. In brief, to support carbon-neutral plans, urban growth modeling should consider the complex network interactions between cities when assigning long-term urban land quantity changes and should adopt spatially heterogeneous rules to allocate urbanized patches in large-scale spatial domains.

Some studies have explored novel solutions for modeling urban growth in a large-scale region. Accessibility, urban flows, and data field approaches such as gravity-based models and potential-based models are usually used to measure the spatial interactions between cities in an urban agglomeration [47,53,54,55,56]. However, these approaches are incapable of describing the network relationships between cities in an urban agglomeration. Moreover, they are usually quantified as spatial allocation rules to affect the distribution of newly urbanized cells, which ignores their direct effects on urban growth quantity at the macro level. On the other hand, the transition rules with spatial heterogeneity are not rare in CA modeling, but they are usually applied to single cities or small-scale regions, such as geographically weighted regression (GWR) [57,58] and machine learning-based search algorithms [59,60]. In addition, zoning or partitioning strategies are usually used to address the spatial heterogeneity of urban growth in some studies [31,47], but they can only characterize spatially heterogeneous urban growth features at the sub-region scale. More importantly, without considering spatial interactions between sub-regions, these strategies may result in the separation of research areas [54]. Furthermore, a recent study has proved that coupling partitioned quantity rules with spatially heterogeneous rules can significantly improve the simulation performance [61]. Therefore, integrating the quantity module based on network relationships and the spatial module based on spatial heterogeneous CA rules would be an appreciable solution to simulate large-scale and long-term urban growth and support carbon-neutral planning.

This paper proposes a framework to carry out the above solution, of which the core consists of two modules, namely, the network automata (NA)-based quantity module and the CA-based spatial module. The former module focuses on describing the network interactions between cities at the macro level, and the latter is aimed at deriving the spatial heterogeneous rules for the spatial allocation of urban cells. Taking the Guangdong–Hong Kong–Macao Greater Bay Area (GHMGBA) as an example, the proposed framework is applied to simulate and predict its urban growth to measure carbon-neutral pressure and provide early warning for regional development.

2. Study Area and Data Sources

2.1. Guangdong–Hong Kong–Macao Greater Bay Area

The Guangdong–Hong Kong–Macao Greater Bay Area (GHMGBA) is the region with the most complex administrative system in the world. The total land area of GHMGBA is about 5.6 × 104 km², and it covers 9 cities on the Chinese mainland, i.e., Guangzhou (GZ), Shenzhen (SZ), Zhongshan (ZS), Huizhou (HZ), Zhuhai (ZH), Foshan (FS), Dongguan (DG), Jiangmen (JM), and Zhaoqing (ZQ), and two special administrative regions, i.e., Hong Kong (HK) and Macao (MO). At the county administrative level, there are 50 counties in the mainland cities of GHMGBA, 18 sub-regions in Hong Kong, and 5 sub-regions in Macao. These regions not only differ in size but also population and socioeconomic development. To make different administrative units comparable, this paper refers to the Nomenclature of Territorial Units for Statistics (NUTS) of the European Union (EU) and the Tertiary Planning Unit (TPU) of Hong Kong and recombines the county-level regions into 36 normative units of spatial planning (NUSP) with the consideration of their land area and urbanization levels (Figure 1) [62,63].

2.2. Data Materials

2.2.1. Urban Land in 2000, 2010, and 2020

This study aims to explore the evolution trend of the GHMGBA urban agglomeration until the year 2050 (the target year of territorial spatial planning in China), which needs at least three periods of land-use data to conduct the calibration and validation processes. Therefore, the urban land-use changes from 2000 to 2010 are used to calibrate the transition rules, and the observed land-use data in 2020 are used for validation. From GlobeLand30 (https://www.globallandcover.com, accessed on 7 June 2021), the land-use datasets for the years 2000, 2010, and 2020 are retrieved. The GlobeLand30 land-use dataset consists of 10 different types of terrain, including artificial surfaces, cultivated land, forests, grasslands, water bodies, shrubland, wetland, tundra, perennial snow and ice, and bare land. The production of the GlobeLand30 dataset is supported by the National Geomatics Center of China (NGCC), and its source includes many 30 m resolution multispectral remote sensing images (TM5, ETM+, and OLI) of Landsat (USA), multispectral images of the China Environment and Disaster Reduction Satellite (HJ-1), and 16 m resolution multispectral images of the China High Resolution Satellite (GF-1). Chen et al. [64] created and developed the web-based system (GLCVal) (http://glcval.geo-compass.com, accessed on 9 June 2021) to support the global collaborative validation of GlobeLand30. With the support of GLCVal and the guidance of the technical specification, about 20 countries have completed the accuracy assessment of GlobeLand30. It turns out that the total accuracy of GlobeLand30 is 85.72% and the Kappa coefficient is 0.82, which indicates the reliability of GlobeLand30 [65]. This study reclassified GlobeLand30 into five categories, including farmland, forests, grassland, water bodies, and built-up land.

2.2.2. Ecological Conserved area Forbidden for Urban Growth

Spatial constraints are an important part of urban growth modeling data, especially when the modeling is aimed at supporting carbon-neutral planning. In urban planning, most spatial constraints originate from ecological protection policies, such as ecological redlines [66,67]. In combination with the geographical conditions of the GHMGBA, the Pearl River and its surrounding wetlands (ecological water), and forest areas with a relative elevation difference greater than 50 m (ecological forest) are considered as the ecological conserved area, in which the land cannot be developed into urban land (Figure 2). In addition, although permanent basic farmland is prohibited from being developed in China, its distribution is not immutable and can be adjusted according to development needs after the approval of the superlative administrative level. As a part of China’s national strategy, urban development in the GHMGBA will inevitably occupy some farmland, so the constraint effects of permanent basic farmland are not considered in the current work. On the contrary, once the regional ecological environment is damaged, it is difficult to recover.

3. Methodology

3.1. Modeling Framework

Figure 3 shows a diagram of the framework of the urban agglomeration growth (UAG) simulation model. As can be seen from the figure, the core of the modeling framework consists of two parts, namely, a quantity module and a spatial module. In the quantity module, each city in the urban agglomeration is regarded as a network node, and all the cities together compose the network structure of the urban agglomeration. The function of network nodes is to control the urban land quantity of corresponding cities, and their state of evolution is modeled by the Markov process and adjusted according to their interactions. In combination with urbanization laws and maximum scale constraints, the long-term spatiotemporal evolution of urban land quantity for each city in the urban agglomeration can be modeled by the network automata. In the spatial module, a partitioned CA model based on GWR is used to conduct the spatial allocation of newly urbanized cells, which can accurately control urban land quantity in each city and incorporate spatial heterogeneity into large-scale urban growth modeling. Through coupling the two modules, urban agglomeration growth can be effectively simulated from both quantitative and spatial aspects. It should be noted that the urban shrinkage caused by depopulation is not considered in this study, as the focus is to explore the trend of urban agglomeration development to get an early warning of carbon-neutral pressure. Meanwhile, the IPCC coefficient is used to estimate the increased carbon emissions caused by urban agglomeration growth and to provide early warnings of the carbon-neutral pressure [13,68].

3.2. UAG Simulation Model

3.2.1. Quantity Module: Network Automata

The network automata model is the main component of the quantity module and is used to predict the synergic development of urban land quantity in cities within the urban agglomeration. The establishment of the quantity module is described as follows. Firstly, the urban agglomeration is divided into several sub-regions, and each sub-region is regarded as a network node. All the network nodes are connected to form the urban agglomeration network. Since a county is the basic administrative unit for implementing land policies in China, it is used as the standard to divide the study area. Specifically, cities with a size similar to counties, such as Shenzhen, Dongguan, and Zhongshan, are directly identified as county-level regions, and Hong Kong and Macao are also identified as county-level regions due to their specialty. Secondly, the initial probability is obtained through urban growth from time t − 1 to time t, and then a Markov chain process is operated with the constraint of urbanization laws to predict the urban growth quantity from time t to time t + 1 [69,70]. Thirdly, the gravity field model is used to measure the radiant potential of each node in the urban agglomeration network at time t [55,56]. This radiant potential is then used as a weighting coefficient to make the urban agglomeration network vibrate and then balance. That is, the total urban growth quantity from time t to time t + 1 will be assigned to each network node according to the radiant potential. The number of vibrations in the simulation can be determined by practical needs. In theory, the more vibrations there are, the more balanced the relationship between network nodes should be. Fourthly, the maximum developed quantity of each network node is estimated with the constraint of the ecological conserved area and is used as the quantity constraint for the urban growth in the corresponding region. It should be noted that the change in ocean area is not considered in this study because sea reclamation is strictly restricted by policies. The network automata model can be mathematically represented as follows.

$$Q_n^{t+1}=Q_n^t+(w_{\alpha }{\alpha }_n^{t+1}+w_{\beta }{\beta }_n^{t+1})\times (u_n^t-u_n^{t-1})\times Q_n^{\max }, Q_n^{t+1}\le Q_n^{\max }, u_n^t=Q_n^t/Q_n^{\max }, w_{\alpha }+w_{\beta }=1$$

(1)

$${\alpha }_n^{t+1}=\frac{\left(\mathrm{g}\right(u_n^t,\lambda )-u_n^t)}{u_n^t-u_n^{t-1}}$$

(2)

$${\beta }_n^{t+1}=\frac{{\displaystyle \sum _m{\omega }_m\times (u_m^t-u_m^{t-1})}}{(u_n^t-u_n^{t-1})},\forall m\ne n$$

(3)

$${\omega }_m=\frac{(u_m^t-u_m^{t-1})}{{\displaystyle \sum _m(u_m^t-u_m^{t-1})}}$$

(4)

where $Q_n^t$ and $Q_n^{t+1}$ denote the urban land quantity of node n in years t and t + 1, respectively. ${\alpha }_n^{t+1}$ and ${\beta }_n^{t+1}$ are, respectively, the urbanization acceleration of node n itself and the effect it receives from the urban agglomeration network. $w_{\alpha }$ and $w_{\beta }$ are the weighting coefficients to balance the above two effects. Their default values are set to 0.5 and can be adjusted according to the testing data. $u_n^t$ is the urban land development ratio of node n at year t and equals the ratio of current urban land quantity $Q_n^t$ to the maximum developed quantity $Q_n^{\max }$. $\mathrm{g}$ is the predicted urban land development ratio of node n in the next year based on the S-curve that describes the urbanization process. $\lambda$ is the moving step of time nodes on the S-curve, with ten steps corresponding to ten years. The value of $\lambda$ should reduce with the iterations of the network automata model to represent the controlling effects of China’s Land Intensive Use Control Policy. This configuration will lead to the hierarchical convergence of the urban land quantity of all network nodes. This means that not all the urban land development ratios will approach 1, and some ratios may converge to other values according to their situations. Overall, if the urbanization of node n is in an accelerated stage, ${\alpha }_n^{t+1}$ is larger than 1, otherwise, it is smaller than 1. ${\beta }_n^{t+1}$ denotes the radiation effect of neighboring regions $m(m\ne n)$ on the focal region n. If region n is the radiation center, ${\beta }_n^{t+1}>1$, which means the weighted urbanization speed of its neighboring regions is larger than its urbanization speed, otherwise ${\beta }_n^{t+1}<1$. ${\omega }_m$ is the weighting coefficient of the radiation effect, which is determined by the inverse distance weighting of m neighboring regions of the focal region n.

As can be seen from the construction of the network automata, nearly all the parameters can be automatically obtained from historical land-use change, except for the effects of the S-curve and the urban agglomeration network [71]. For the former, the urbanization S-curve of the GHMGBA is fitted based on its urbanization process. That is, the urban land development ratios of all the nodes in the years 2000 and 2010 are sorted into a sequence and then extrapolated to 1 according to their tendency (Figure 4). With the help of Figure 4, the application of the g function can be explained, as if the current urban land development ratio of a node is close to the S-curve value at time t (0.1792), its urban land development ratio at the next stage should be adjusted to be close to the S-curve value at time t + 1 (0.3320). Based on that, the urbanization acceleration ${\alpha }_n^{t+1}$ can be obtained through the driving effect of S-curve laws, which can avoid linear urban growth.

In terms of the effect of the urban agglomeration network, the Generate Spatial Weights Matrix tool in ArcGIS is used to calculate the inverse distance weights of neighboring regions on the network node n, which is used to calculate ${\beta }_n^{t+1}$. The interaction network between the nodes in the GHMGBA in 2050 is obtained after 4 iterations and shown in Figure 5. In Figure 5, the yellow circle represents the urban land quantity in 2010, while the red circle represents that in 2050. All the circles have been proportionally processed to be comparable. The larger the circle radius is, the larger the urban land quantity should be, and the greater the difference between red and yellow circles is, the larger the urban growth quantity should be. On the contrary, minor differences between red and yellow circles denote that urban growth is reaching its limit or has no more potential. The blue lines in Figure 5 represent the interaction intensity between network nodes; the thicker the line is, the stronger the synchronization between two nodes and the stronger the effect of the urban agglomeration network on the node should be.

3.2.2. Spatial Module: Cellular Automata

The network automata can well predict the urban land quantity of each node (sub-region). However, the urban land quantity needs to be accurately allocated in the spatial domain to further complete the urban growth simulation process. Therefore, the partitioned CA framework is adopted to build the spatial module to reproduce and predict land-use patterns in this study, which can be described as:

$$p_{ij}^t=S_{ij}·{\mathsf{\Omega }}_{ij}^{t-1}(X_W^{t-1}\oplus \Delta x_W)·{\delta }_{ij}^{t-1}·{\mathsf{\Psi }}_{ij},{\mathsf{\Psi }}_{ij}\in \{0,1\}$$

(5)

$$X_{ij\in n}^{t+1}=\mathrm{CA}(p_{ij\in n}^t,Q_n^{t+1}),X\in \{0,1\}$$

(6)

where $p_{ij}^t$ is the transition probability of cell ij to be developed as an urban cell at time t, which consists of transition potential $S_{ij}$, neighborhood effect ${\mathsf{\Omega }}_{ij}^{t-1}(\varphi )$, the constraint of ecological conserved areas ${\mathsf{\Psi }}_{ij}$, and the stochastic perturbation ${\delta }_{ij}^{t-1}$. X denotes the cell state, with 1 meaning urban land and 0 meaning non-urban land. Among them, ${\mathsf{\Omega }}_{ij}^{t-1}$ is calculated by a moving-window function, $X_W^{t-1}$ denotes the state of cell ij in its Moore neighborhood W, and $\Delta x_W$ denotes the key development area. $\Delta x_W$ is a binary variable, and $\oplus$ is an integrative operator. Combined with them, the key development cells in the neighborhood will be labeled as 1 to strengthen the effect of spatial planning. ${\mathsf{\Psi }}_{ij}$ is the static constraint of ecological conserved areas, and ${\mathsf{\Psi }}_{ij}=0$ means that the objective cell is an ecological conserved cell that cannot be developed. ${\delta }_{ij}^{t-1}$ is a stochastic function to generate random numbers in iterations to represent the uncertainty in realistic urban growth. $X_{ij\in n}^{t+1}$ denotes the state of cell ij in the sub-region n at time t, which can be determined by $p_{ij}^t$ and the partitioned strategy $Q_n^{t+1}$.

The transition potential $S_{ij}$ is the core parameter of CA modeling [25,32,38]. It represents the probability of a cell being an urban cell and is the foundation of grid-based urban growth simulation. Generally speaking, some natural and anthropogenic factors need to be considered in the derivation of the transition potential distribution, such as topographic features, transport networks, population density, and so on [33]. However, several anthropogenic factors, like transport networks and population density, are unstable factors for large-scale and long-term urban growth, which will have significant effects on urban development under the management of China’s government. For example, a newly built railway or super highway will change the driving mechanisms, and attenuating the constraints of the household registration system will attract a large-scale population inflow. Therefore, some relatively stable factors are selected to be spatial driving variables in the development of urban agglomeration, including (a) elevation, (b) slope, (c) distance to the main canal (the Pearl River), (d) distance to the core nodes (Guangzhou, Shenzhen, Hong Kong, and Macao) of the urban agglomeration network, (e) distance to several important nodes, and (f) distance to other nodes. These spatial driving variables are used as independent variables in the derivation of transition potential distribution, and the urban growth of the GHMGBA during 2000–2010 is used as the dependent variable.

The relationship between the dependent variable and independent variables is always interpreted by data mining approaches, in which logistic regression (LR) is the most widely used [27,28]. LR modeling defaults the driving mechanisms to being spatially homogeneous, which makes each spatial driving variable have only one spatially stationary regression coefficient [57]. However, the spatial heterogeneity of urban growth would become more significant on a large scale than on a small scale, as each city may have its characteristics. Therefore, GWR is adopted to estimate the transition potential of each cell to incorporate the spatially heterogeneous features of urban growth. This study resamples the binary change (0 means unchanged and 1 means urbanized) of 30 m-resolution urban growth from 2000 to 2010 into a continuous density distribution with a resolution of 100 m. The spatial driving variables are also linearly normalized to a range between 0 and 1 to participate in GWR modeling. Furthermore, the probability distribution obtained from GWR modeling needs some modification to be the transition potential distribution, including (1) assigning 0 to the cells in ecological conserved areas, (2) assigning 1 to existing urban cells, and (3) inserting key development areas and integrating them into $\Delta x_W$ to steer the urban growth trend. In the historical development of the GHMGBA, areas near line A (in Figure 6) are the focus areas of urban growth, while the new development strategy proposed by the government plans to promote the development of areas near line B and coordinate with areas near line C. Since the urban land distribution had few changes near line B from 2000–2010, it is difficult for GWR modeling to acquire a probability distribution with high-probability cells near line B. Therefore, some subjective experiences and inferences need to be integrated into the transition potential distribution as no model can perfectly predict future dynamics. This study labels the nodes near lines B and C and makes their $\Delta X=1$. With the above derivations and modifications, the final transition potential distribution with key development areas has been obtained and is shown in Figure 6.

3.3. Carbon Emissions Caused by UAG

Although the process of carbon emissions caused by land-use change is complex, this study focuses on the change in carbon emissions caused by urban growth and its encroachment on other land uses. Therefore, carbon emissions caused by land-use change are estimated using the coefficients proposed by the IPCC [13,68,72].

$$E_n(k)=A_n(k)·B(k)$$

(7)

where $E_n(k)$ is the carbon emission of the k-th land use in region n, $A_n(k)$ is the area of the k-th land use in region n, and $B(k)$ is the carbon emission coefficient of the k-th land use. A recent study has estimated the carbon emissions related to land use in the GHMGBA with a combination of remote sensing data and socioeconomic data. It has estimated the carbon emissions from different land-use types in the GHMGBA based on the logarithmic mean Divisia index (LMDI), which is one of the most effective and commonly used methods for analyzing the factors affecting carbon emissions. According to its results for 2010, if the cultivated land, forest land, grassland, and water bodies are developed into urban land, the reduced carbon sink should be 397.09 t/ha, 419.71 t/ha, 393.31 t/ha, and 459.21 t/ha, respectively. Although the implementation of carbon-neutral planning may decrease these values from 2010–2050, they can still be used as references for early warning.

4. Results and Analysis

4.1. Predicted Urban Land Quantity

The prediction accuracy of the network automata model has been tested before applying it to predict urban land quantity. Four prediction models are established to investigate the performance of the proposed network automata model, including (1) ErS1: urban growth is driven by the historical inertia of 2000–2010; (2) ErS2: urban growth is driven by the historical inertia of 2000–2010 and the S-curve of urbanization; (3) ErS3: urban growth is driven by the historical inertia of 2000–2010 and the effects of urban agglomeration network; and (4) ErS4: urban growth is driven by the three factors, i.e., the proposed network automata. With the four models, the urban land quantity of 36 sub-regions in 2020 has been simulated and is compared to the observation in 2020 to show their errors (Figure 7). Ideally, if the model can perfectly predict the urban land quantity of each sub-region, the error line should be a straight line with a value of 0. As can be seen from Figure 7, the error lines of the four models are fluctuant, and their mean error rates are 17.39%, 15.82%, 12.86%, and 9.83%, respectively. The results show that the proposed network automata model has the best performance in predicting urban land quantity, and its error rate is less than 10%, which means that it is more reliable than the other models in predicting urban land quantity and can be used to predict the urban land quantity of each sub-region in 2050.

The tested network automata model is used to predict the urban land quantity in each sub-region for the period 2020–2050, which is shown in

Table 1. Predicted urban land quantity of sub-regions in 2030, 2040, and 2050 (km²).

ID	Sub-Region	Q2030	Q2040	Q2050	ID	Sub-Region	Q2030	Q2040	Q2050
1	GZ-1	512.13	545.92	563.01	19	HZ-3	339.55	395.58	428.36
2	GZ-2	390.54	414.92	419.16	20	HZ-4	365.26	424.27	462.13
3	GZ-3	445.5	474.98	489.91	21	HZ-5	96.57	126.73	143.86
4	GZ-4	227.84	261.28	277.29	22	JM-1	209.98	228.38	238.05
5	GZ-5	407.01	447.05	466.74	23	JM-2	313.36	357.79	381.08
6	GZ-6	437.93	492.96	521.11	24	JM-3	198.74	227.4	242.85
7	GZ-7	193.11	229.04	248.81	25	JM-4	378.17	438.71	478.2
8	SZ	1278.72	1278.72	1278.72	26	JM-5	227.32	267.78	291.12
9	ZH-1	226.46	242.04	245.71	27	JM-6	155.52	180.48	197.85
10	ZH-2	391.22	436.13	453.29	28	ZQ-1	116.93	129.24	136.54
11	FS-1	929.31	988.28	1018.69	29	ZQ-2	200.68	229.45	245.55
12	FS-2	576.59	617.08	637.86	30	ZQ-3	218.59	254.41	276.76
13	FS-3	301.17	337.79	356.48	31	ZQ-4	59.16	69.04	75.17
14	FS-4	173.73	197.56	210.64	32	ZQ-5	64.18	76.15	83.39
15	DG	1821.86	1920.35	1920.35	33	ZQ-6	101.41	119.98	131.77
16	ZS	858.94	941.9	980.3	34	ZQ-7	50.87	58.88	65.2
17	HZ-1	368.02	410.36	434.26	35	HK	380.82	404.98	410.72
18	HZ-2	389.66	426.21	440.19	36	MO	25.83	26.02	26.02

. As shown in the prediction results, the urban land in Shenzhen will stop growing in 2030, and the urban land in Dongguan and Macao will stop growing in 2040. This is because their area is relatively small but they have experienced a rapid urbanization process before 2020, which leads to a high radix of urban land development ratio for them. The predicted urban growth during 2040–2050 is far less than that during 2030–2040 in several sub-regions such as GZ-2, GZ-4, GZ-5, ZH-1, ZH-2, ZS, HZ-2, and HK, where their decrease rates are up to 50%. Furthermore, the urban land development ratios of sub-regions in different periods are diagrammed to investigate the urban growth characteristics of the GHMGBA (Figure 8). As a whole, although the urbanization of China has entered a later stage with no more large-scale urban growth, the GHMGBA still has its own development needs. The quantity growth of urban land in sub-regions of the GHMGBA has an obvious trend of hierarchical convergence. After quartering the urban land development ratios, it can be found that in 2050, there are 10 sub-regions ranked in the first level (0.75~1), 8 sub-regions ranked in the second level (0.5~0.75), 11 sub-regions ranked in the third level (0.25~0.5), and 7 sub-regions ranked in the fourth level (0~0.25). These characteristics indicate that hierarchically converging is an inevitable trend for the prediction of urban land quantity in urban agglomerations, especially in the case of implementing the intensive land-use system in China. Combined with the above results, it can be speculated that the urban growth in sub-regions is different, and there are several sub-regions in the GHMGBA that will reach the later stage of urbanization in 2050. Nevertheless, the total urbanization of the GHMGBA is still in an acceleration phase and will enter the stable stage in 2035~2050, which is partly because the implementation of the “Development Plan for the GHMGBA” will provide new impetus for its development. Finally, according to the prediction results in this study, all the urban land area in the GHMGBA will be about 15,277.14 km² in 2050, which accounts for 27.28% of the total land area and is lower than the red line (30%).

4.2. Predicted Land-Use Patterns in 2050

Based on predicted urban land quantity growth, urbanized cells can be spatially allocated through the spatial module, i.e., the CA model, to complete the urban agglomeration growth simulation. The simulation performance of the CA model is detected through the kappa index (Kappa) and figure of merit (FoM); the former measures the quantity consistency between two maps, while the latter measures the change consistency. Based on simulated and observed land-use maps in 2020 (Figure 9), the Kappa of simulated urban land is calculated to be 0.4416, which seems to be not high enough. This is because Kappa is easily affected by the proportion of urban land to the total area, and urban land accounts for a small percentage of the GHMGBA [61,73]. FoM is relatively better than Kappa in terms of reflecting the simulation performance as it focuses on the urban land change in the simulated and actual urban growth [74]. It is obtained by calculating the ratio of intersection and union of simulated and observed urban growth and is 0.2926 in this study, which is satisfactory and indicates the good performance of urban growth modeling.

The CA model is used to predict the future land-use patterns in the GHMGBA in 2050 (Figure 10) after testing the capability of the CA model to simulate urban growth. As can be seen from Figure 10, depending on the key development areas prescribed in the “Development Plan for the GHMGBA”, the urban land patches of the GHMGBA will be evenly distributed in 2050 and gradually form an obvious “Ω” structure, which can be confirmed by the significant development of western sub-regions. Furthermore, five location cases are selected to verify the effectiveness of predicted land-use patterns, including (a) the development zone in Zhaoqing, (b) Sino-Singapore Guangzhou Knowledge City, (c) the central urban area of Huizhou, (d) the Nansha Development Zone in Guangzhou, and (e) the Doumen district of Zhuhai. The five locations in the predicted map (2050) are enlarged in Figure 11 and compared with the base map, which is the latest high-resolution satellite image of 2020. The comparison shows that road networks or bulldozed ground happen almost in all of these locations, indicating they are experiencing the urbanization process. These findings indicate that the proposed framework has a good ability to predict the urban development of the GHMGBA. This is because the network automata can reconstruct the development mode of the urban agglomeration. Although the effect of this reconstruction cannot be observed in the short term (e.g., urban land-use change from 2010–2020), it would make a significant difference in long-term prediction when combined with government planning.

4.3. Early Warning of the Carbon-Neutral Pressure

With the predicted land use map, the loss of carbon-sink land can be detected (Figure 12a) and the carbon emissions caused by urban agglomeration growth can be further estimated. As can be seen from the figure, most of the lost carbon-sink land is distributed in the middle area of the GHMGBA, such as Guangzhou and Foshan. Combining the estimated carbon emission intensity of each land-use type in 2010 [72] and the simulated land-use change between 2020 and 2050, the total carbon emissions caused by urban growth in the GHMGBA until 2050 is about 2.17 × 10⁸ tons. These carbon emissions are deeply related to land-use distribution, which makes the carbon-neutral pressure differ in the spatial domain. A statistic on the land-use-related carbon emissions has been carried out at the city level, and the results are graded by the natural breaks method (Jenks) and shown in Figure 12b. For a city, the more carbon emissions there are, the more pressure there is to be carbon neutral. The distribution of carbon-neutral pressure at the city level has an obvious feature. That is, the carbon-neutral pressure in cities in the central region is relatively large, such as Guangzhou, Foshan, Huizhou, and Jiangmen, and they form a “Ω” structure. Although a large amount of carbon emissions lead to high carbon-neutral pressure, it also shows that these cities have great development potential. In addition, the carbon-neutral pressure of Dongguan, Zhongshan, and Zhaoqing is at a medium level, and the carbon-neutral pressure of Hong Kong, Macao, Shenzhen, and Zhuhai is at a low level, which indicates that these four cities are likely to take the lead in achieving the carbon-neutral goal. By estimating carbon emissions from projected urban growth, the early warning of carbon-neutral pressures can be provided to support carbon-neutral targets. This early warning approach allows the GHMGBA sufficient time to formulate development plans to achieve its carbon-neutral targets, such as upgrading its industrial structure, limiting uncontrolled urban growth, and increasing carbon-sink land.

5. Discussion

With the deterioration of the climate and environment, governments and scholars have placed an increasing emphasis on the negative effects of carbon emissions and expect to accurately quantify them. The sources of carbon emissions are diverse, but most of them are related to anthropogenic activities, in which land-use change is a typical representation of human–environment interactions [15]. Land-use change has a great impact on carbon emissions, not only because of urbanization and consequent industrialization but also the quantity and quality change of carbon-sink land [75]. In particular, the increase in urban land and the decrease in forestland, grassland, and cultivated land together aggravate carbon emissions [76]. Therefore, reasonable land-use planning is an effective way to achieve carbon neutrality. Even though numerous models have been developed to support land-use planning, there are still some issues to be addressed before applying them to estimate the carbon emissions caused by land-use change.

Firstly, achieving carbon neutrality is usually a long-term and large-scale problem of regional planning, in which many counties, cities, and even provinces are covered [44,45]. Although there are more and more studies modeling land-use change in a large-scale region, the methods they use are always similar to those in small regions. With these methods, it is usually difficult to depict large-scale land-use change. This is because large-scale urban growth modeling needs to integrate regional interactions into the simulation process, which is often omitted in urban growth modeling for small regions [47,53]. Some scholars have recognized this issue and proposed several solutions to model large-scale land-use change, such as incorporating urban flows and data fields into land-use modeling [54,55,56]. However, these frameworks focus on the micro-allocation of cells. That is, they tend to quantify or integrate the regional interactions into a spatial distribution to direct the state change of cells. Although regional interactions do have an impact on the state change of cells, they have a more direct and explicit effect on the scale development of cities within the large-scale region. It is necessary to associate regional interactions with urban land quantity to depict the complex relationships between sub-regions in large-scale urban growth modeling. Furthermore, the interaction among sub-regions is not a simple bidirectional effect or a hierarchical structure, it is a complex network coupling the above two [77]. Therefore, this study proposes a network automata model to characterize the complex long-term coevolution of urban scale in urban agglomeration. This model framework integrates the historical inertia of urban development, urbanization law (the S-curve), and the effects of the urban network. It is a novel solution to enhance the ability to simulate and predict long-term quantity changes in large-scale regions, thus improving the suitability of land-use models to support carbon-neutral planning.

Secondly, another problem will emerge as the region of urban growth modeling becomes larger. That is, the difference between the development of sub-regions will be more distinct in large-scale regions. This difference is mainly manifested in two aspects, including quantity and driving mechanisms. As for the former, each sub-region has its own urbanization stage, so different amounts of urban land are needed in different sub-regions [31]. A partitioned modeling framework with accurate quantity control should be coupled with the network automata model to describe the spatial heterogeneity of urban growth quantity. Consequently, an appropriate zoning method is important for the rationality and practicability of large-scale urban growth modeling. Recent studies on urban growth modeling have introduced several zoning methods based on different features, including topographical features, land-use attributes, and administrative division [78]. However, zoning based on topographical features or land-use attributes is defective in practicability for land-use planning, and zoning based on the administrative division is greatly affected by the size and administrative level of different sub-regions. Therefore, this study adopts the ideas of NUTS (EU) and TPU (HK) to divide the GHMGBA into some sub-regions with similar sizes. In terms of the latter, many studies expect to characterize the spatial heterogeneity of driving mechanisms in a similar zoning way. However, zoning can only depict the macro differences of sub-regions, and it is difficult to explore the spatial heterogeneity of driving mechanisms at the micro-scale. In addition, a recent study has proved that coupling the zoning strategy, accurate quantity control, and GWR modeling will result in a CA model with better simulation performance [61]. Therefore, this study uses the GWR model to interpret the spatially heterogeneous driving mechanisms at the micro-scale. By coupling the partitioned quantity control based on the network automata model and the driving mechanisms derived by the GWR model, the spatial heterogeneity of large-scale urban growth can be well interpreted in simulation and prediction, thus improving the effectiveness of supporting carbon emission estimation.

Thirdly, no matter how good the performance of urban growth modeling is, it is always difficult to completely model the real process, especially in the prediction of future land-use patterns. This is because urban growth is a process greatly affected by human decisions [79]. Many scholars expect to model the decision-making process of humans to enhance the authenticity of the simulation. ABM modeling is a good example [46]. However, quantifying the decision-making process of humans to be model rules is very challenging. That is also one of the reasons why the ABM model is not the most popular land-use model. In addition, CA modelers have tried to incorporate some human decision-making into urban growth modeling, such as controlling urban land quantity according to land-use policies, adopting boundaries to restrict disorderly urban growth, and generating random seeds in the planning area to guide the change of land-use patterns [80]. However, these approaches cannot be integrated into the self-organization of urban growth, which is the key process for producing and directly affecting land-use patterns. Therefore, this study proposes a novel method to improve the neighborhood component that characterizes self-organized urban growth by incorporating the distribution of key development areas into its definition. This optimization will strengthen the effect of land-use planning in the simulation and prediction of self-organized urban growth, thus improving the ability of CA models to produce desired land-use patterns.

6. Conclusions

This paper proposes a framework based on urban growth modeling to provide early warnings about the carbon-neutral pressure of urban agglomerations. The framework has two core modules, the quantity module and the spatial module, to accomplish urban growth simulation and prediction for large-scale regions over a long time period. The former integrates the inertia of historical land-use change, the driving effects of the urbanization law (S-curve), and the traction of the urban agglomeration network to model the long-term quantity change of urban land in a large-scale region. The latter couples a partitioned modeling framework, spatially heterogeneous rules, and quantified land-use planning orientation to build a CA model to accurately allocate the urbanized cells in a large-scale spatial domain. These two modules are both designed to improve the ability of the land-use model to simulate urban growth on a large scale and over the long term, which are both basic characteristics of carbon-neutral planning.

Taking the GHMGBA as an example, the historical urban growth from 2000 to 2010 is used to calibrate the proposed framework, and its performance is validated by simulating the urban growth from 2010 to 2020. The calibration and validation results indicate that the proposed framework is reliable and well-performed both in the modeling of urban land quantity and the spatial allocation of urbanized cells. With the proposed framework and IPCC coefficients, the carbon emissions caused by the urban growth of the GHMGBA from 2020 to 2050 are estimated. The results show that about 2.17 × 10⁸ tons of carbon will be released as a result of the urban growth of the GHMGBA due to two aspects: the increase in carbon-source land and the decrease in carbon-sink land. Among the cities in the GHMGBA, Guangzhou, Foshan, Huizhou, and Jiangmen are under great pressure to achieve carbon-neutral targets in the future, while Hong Kong, Macao, Shenzhen, and Zhuhai are relatively easy to bring up to the standard.

However, there are still some limitations in this study. First, urban growth is not only a typical geographical process with significant spatial heterogeneity but also the evolution of complex geographical systems, which have obvious nonlinearity. The proposed framework is unable to capture the nonlinear driving mechanisms of urban growth. Second, the sources of carbon emissions are very diverse. This paper only models the carbon emissions caused by land-use change with remote sensing data and socioeconomic data. In the future, a comprehensive framework should be developed that incorporates additional sources into the modeling to accurately predict carbon emissions. Finally, although the performance of the proposed framework has been validated and assessed by the evidence, its generalizability should be further tested in other study areas, and its performance needs a comprehensive comparison with other techniques.