Multi ‐ Scenario Simulation of Urban Growth Under Integrated Urban Spatial Planning: A Case Study of Wuhan, China

: Although many publications have noted the impact of urban planning on urban develop ‐ ment and land ‐ use change, the incorporation of planning constraints into urban growth simulation has not been adequately addressed so far. This study aims to develop a planning ‐ constrained cellu ‐ lar automata (CA) model by combining cell ‐ based trade ‐ off between urban growth and natural con ‐ servation with a zoning ‐ based planning implementation mechanism. By adjusting the preference parameters of different planning zones, multiple planning ‐ constrained scenarios can be generated. Taking the Wuhan Urban Development Area (WUDA), China as a case study, the planning ‐ con ‐ strained CA model was applied to simulate current and future urban scenarios. The results show a higher simulation accuracy compared to the model without planning constraints. With the weak ‐ ening of planning constraints, urban growth tends to occupy more ecological and agricultural land with high conservation priority. With the increase in preference on urban growth or natural conser ‐ vation, the future urban land pattern will become more fragmented. Furthermore, new urban land beyond the planned urban development area can be captured in future urban scenarios, which will provide certain early warning. The simulation of the current urban spatial pattern should help plan ‐ ners and decisionmakers to evaluate the past implementation of urban planning, and scenarios sim ‐ ulation can provide effective support for future urban planning by evaluating the ( j ) cultivated


Introduction
Urbanization is one of the greatest social transformations in human history [1,2]. In order to meet the unprecedented demand of urban population growth and economic development, urban land presents a violent and disorderly expansion trend [3]. Nowadays, the average growth rate of urban land is twice that of their populations [4]. Although urban land is only a small part of the Earth's surface, its rapid growth has been identified as a major contributor to natural habitat loss, cropland reduction and ecosystem services degradation [5][6][7], posing great threat to the world's sustainable development [8,9]. Thus, it raises continuous competition for land between urban development and natural conservation [10], which is especially relevant to China due to its limited land resources and huge population living in urban areas [11]. Currently, Chinese local governments have put forward a series of administrative regulations and spatial planning policies, including major function-oriented zoning, land-use planning, urban spatial regulations, ecological redlines and so on. These policies and regulations have played a vital role in restraining the decrease in ecological and agricultural land and guiding urban compact development [12][13][14], although illegal urban encroachment on protected areas frequently occurs [15,16]. Therefore, simulation of future urban expansion and encroachment incorporating different planning policies is of great significance for authorities and planners to make effective decisions in optimizing urban spatial patterns and achieving urban sustainable development.
Since the first application of Tobler in 1979 [17], cellular automata (CA) and its derived models have been widely used in urban growth modeling and land-use simulation due to their simple and flexible structure, among which the transition rule is a core part [7,18,19]. The process of urban growth and urban land change is influenced by multiple drives across different scales, including socio-economic factors, geophysical factors, proximity factors, environmental factors and planning factors [18,[20][21][22]. In order to link these multiple driving factors with the transition rules of the CA model, scholars have explored various methods including multivariate statistics methods, intelligent algorithms, agentbased models and multi-criterion decision analysis [3,[23][24][25]. Although the combination of these models and cellular automata has made many important achievements in urban land-use simulation, the applicability of these models in reproducing urban patterns resembling reality and solving actual planning problems in complex urban systems is still challenging [18,20].
In recent decades, many researchers have increasingly recognized the important influence of urban planning and other policies on urban growth and land-use change [15,26,27]. Thus, scenario simulation of urban land use incorporated with planning information and conventional urban CA models has attracted considerable attention from scholars and planners. Through the constraint layers or exclusion layers in urban CA models, policy information including supportive planning and restrictive planning can be implemented to simulate and predict urban growth by either a Boolean-constraint mechanism or a gradual mechanism [27][28][29]. With regard to Boolean-constraint, for example, He et al. consider the primary farmland protection area as mandatory constraints to prohibit urban encroachment [12]; Zhou et al. set water areas, historical protection areas, restricted development areas, prohibited development areas and large-scale parks in planning as constraint layers to simulate urban spatial patterns under the planning intervention scenario [3]. While the binary planning constraint design is too absolute and simplistic [30], some studies have attempted to simultaneously consider the gradual constraints into their urban CA models. Onsteda et al. incorporated zoning policy information into the CA-based SLEUTH model by assigning gradual values to the different zoning categories, which obviously improves model performance with respect to a model without zoning [26]. Liang et al. proposed a random seeding mechanism to simulate the potential effect of planning development zones [20]. Wang et al. integrated a series of urban spatial planning to construct synthetical planning constraint layers in urban CA model by a multi-criteria evaluation (MCE)-based weighting method [31]. These studies have achieved remarkable results and provide good references for building a realistic planning constraint mechanism.
However, in previous studies, the planning constraints in a planning zone are often set to a unified value, ignoring the spatial heterogeneity of different cells in the same planning zone [32,33]. In the real complex urban systems, the planning influence on the cells in a specific planning zone is often not homogeneous. The application of planning constraints in urban CA models should not be restricted to considering the zone effects; it should take the natural environment and socio-economic factors based on the land-use cell into account [34][35][36]. For example, in the planned ecological protection area, cells with high ecological importance tend to be retained, while cells with low ecological importance located in urban fringe areas are more likely to be occupied. Currently, the influence effects of planning policies based on the trade-off between urban growth and conservation has received little attention. In addition, existing studies often set the restriction or impact coefficients of different planning zones based on subjective or historical experience [27,37]. Once these parameters are determined, they will not change with time in the simulation process, which is not conducive to the dynamic urban growth scenario simulation.
This study aims to explore the potential influence of integrated urban spatial planning on urban growth and realize the planning-constrained urban multi-scenario simulation. Thus, we propose an innovative synthetical planning-constrained strategy to quantitatively characterize the influence of integrated urban spatial planning on realistic urban growth process by combining the cell-based trade-off between urban growth and conservation with a zoning-based planning implementation mechanism. By incorporating the strategy into urban CA models, a planning-constrained CA model is developed. Taking Wuhan city, one of the mega cities in China, as a case study, we applied the planningconstrained CA model to assess the impact of urban planning implementation on the urban growth process and simulate the urban spatial patterns under different future scenarios. The results should help planners and authorities evaluate the guiding effects of urban planning implementation.
The rest of the sections of this article are organized as follows. Section 2 describes the study area and the datasets used in this study. Section 3 introduces in detail the methodology of the planning-constrained CA model. Section 4 presents the multi-scenario simulation results under planning constraints. The discussion and conclusion are given in Sections 5 and 6.

Study Area
Wuhan, the capital city of Hubei Province, is the key city in Central China and one of the core cities in the Yangtze River economic belt. As shown in Figure 1, Wuhan administers seven central districts (Hanyang, Hongshan, Jiangan, Jianghan, Qiaokou, Qingshan, Wuchang) and six suburban districts (Caidian, Dongxihu, Hannan, Huangpi, Jiangxia, Xinzhou), covering a total area of 8569.15 km 2 . In past decades, Wuhan has experienced an unprecedented urbanization process. By the end of 2020, the total population of Wuhan was 12.32 million, of whom 10.39 million were urban residents. Rapid urbanization has led to a large amount of urban land expansion and encroachment on other lands, resulting in strong contradictions between urban growth and natural conservation [38,39]. From 2002 to 2019, the area of urban construction land increased by 98.7%, from 1338.75 km 2 to 2646.26 km 2 . Therefore, Wuhan has formulated a series of spatial planning policies, which have greatly affected the process of urban growth [11,13,28]. In this study, we will specifically target the Wuhan Urban Development Area (WUDA), which was designated as the centralized area of urban development in the 2010-2020 Wuhan Master Plan. The WUDA spans a total area of 3261 km 2 , consists of all central districts and parts of the six suburban districts.

Datasets
The datasets used in this study consist of land-use, socioeconomic, natural eco-environment and planning datasets, including the following: . These planning datasets were integrated to produce urban growth scenarios based on dynamic planning constraints.
All these datasets were preprocessed and rasterized with a spatial resolution of 30 m in ArcGIS 10.2 software as input of the urban CA model.

Methodology
A planning-constrained CA model is developed in this study to simulate urban growth and encroachment under different scenarios. Figure 2 shows the model procedure and flowchart. Firstly, cell-based conservation priority and urban growth potential can be calculated by multi criteria analysis and a logistic regression model, respectively. In terms of conservation priority, it consists of ecological land conservation priority and cultivated land conservation priority. The ecological priority is mainly based on relevant indicators such as ecosystem functional importance and ecological vulnerability, while the cultivated land conservation priority mainly depends on the agricultural production suitability of cultivated land. The urban development potential is determined by natural, socioeconomic and location-based factors. Secondly, we defined the synthetical planning-constrained mechanism by integrating the cell-based tradeoff between urban growth potential and conservation priority with zoning-based planning preference. By adjusting the preference parameters of different planning zones, multiple planning-constrained urban growth scenarios can be generated. Thirdly, we combine the synthetical planning-constrained mechanism with urban development potential, neighborhood effect and random disturbance to construct the transition rule of urban CA model. Markov chain is used to predict the amount of urban land-use change. Kappa and FoM coefficients are used to calibrate and verify the model. Finally, future urban growth pattern under different scenarios are simulated. These phases are presented separately as following.

Conservation Priority Evaluation
Urban growth is often at the expense of the loss of ecological land and cultivated land [12,40,41]. In this study, conservation priority includes ecological land conservation priority (ELCP) and cultivated land conservation priority (CLCP). For an ecological land unit such as forest land, grassland, water area and garden land, its conservation priority depends on relevant factors of ecosystem functional importance and ecological vulnerability [42,43], while for a cultivated land unit, its conservation priority should also consider agricultural production suitability [44].
The ecosystem functional importance refers to the ecological conditions and utility formed in the ecosystem and its ecological process, which are conducive to human survival and development [45,46]. According to the ecological status of the study area, biodiversity conservation (BC), carbon storage (CS), water conservation (WR), soil conservation (SC) and flood regulation (FR) are selected as the factors for the evaluation of ecosystem functional importance. Ecological vulnerability refers to the sensitive response and self-recovery ability of ecosystem relative to external interference at a specific time and space scale, wherein a higher vulnerability indicates worse self-recovery ability of the ecosystem [47]. Here, ecological vulnerability includes soil erosion sensitivity (SES), land desertification sensitivity (LDS) and stony desertification sensitivity (SDS). The acquisition and calculation of ecological related factors are mainly based on the Technical Guide for Delimitation of Ecological Protection Red Line issued by the Ministry of Environment Protection of China in 2015 and partly adjusted according to the study area. Table 1 presents the methods for ecosystem functional importance evaluation. Cultivated land quality grade data were used to describe the agricultural production suitability. Then, through ArcGIS 10.2 software, each factor was reclassified to values ranging from 0 to 1, by using the standardized function. The higher the value, the higher the conservation priority. The weight of each factor is obtained by the analytic hierarchy process (AHP) model. Finally, the conservation priority equation can be expressed as: where Ci is the conservation priority for cell i, ranging from 0 to 1; we and wc are the weights of ELCP and CLCP, respectively; we = 1.0 and wc = 0.0 for an ecological land unit, while we = 0.5 and wc = 0.5 for a cultivated land unit; xk represents the value of kth factor, and bk is the corresponding weight. M is the total number of variables. Figure 3 illustrates the results for conservation priority evaluation

Models Parameter explanation
Biodiversity conservation (BC) The habitat quality model of InVEST: 1 is the habitat quality index of patch group x in (dimensionless). and are the habitat suitability score and the total stress level of grid x in , respectively. k is the scale factor (constant).

Carbon storage (CS)
Carnegie-Ames-Stanford approach model (CASA): APAR is the ecosystem net primary productivity of pixel x in tst month (gC•m-2•month-1), and are its corresponding absorbed photosynthetic effective radiation and actual light energy utilization, respectively. , , and are the total solar radiation, the absorption ratio of vegetation layer to incident photosynthetic effective radiation, the influence coefficient of water stress of pixel x in tst month, respectively. T1 and T2 indicate the stress effect of low and high temperature on light energy utilization , is the maximum light energy utilization under ideal conditions. ; SC is the quantity of soil conservation (t/hm2•a), and are the quantity of potential and actual soil erosion, respectively. R, K, LS, C and P are the rainfall erosivity, soil erosion factor, slope factor, vegetation cover factor and the soil conservation practices factor, respectively. Flood regulation (FR) The assessment is based on the water level and area of the water body.
The main rivers of the Yangtze River and Han River, large lakes (>1 km2), large-and medium-sized reservoirs are extremely important flood control ecological function areas; others belong to important areas. Soil erosion sensitivity (SES) √ C SES is the sensitivity index of soil erosion (dimensionless). R, K, LS and C are the sensitivity grades of rainfall erosivity, soil erosion, slope factor and vegetation cover, respectively.
Land desertification sensitivity (LDS) √ LDS is the sensitivity index of land desertification (dimensionless). I, W, K and C are the sensitivity grades of dryness index, sand-driving wind days, soil texture and vegetation cover, respectively.
√ SDS is the sensitivity index of stony desertification (dimensionless). D, S and C are the sensitivity grades of exposed area percentage of carbonate stony, terrain slope and vegetation coverage, respectively.

Urban Growth Potential Evaluation
While urban growth is affected by various factors, the literature shows that the following factors are crucial: natural environment, socio-economic status and location conditions [22,34,48]. A variety of methods, such as multivariate statistics, intelligent algorithms and multi-criterion analysis, were used by researchers to define the contribution of diverse driving factors to historical urban growth [18,[49][50][51]. In this study, a widely used binary logistic regression (BLR) method [6] was used to produce the urban growth potential map. The BLR equation can be expressed as: where Ui is the urban growth potential for cell i, ranging from 0 to 1. B0 is the constant, xk represents the kth independent variable, and bk is the corresponding weight. N is the total number of independent variables. Based on previous researches [15,51], we selected 14 driving factors, such as dem, slope, distance from water area, distance from highway, distance from main road, distance from secondary road, distance from branch road, distance from subway station, distance from city center, distance from parks, distance from commercial center, land value, POI density and population density, to evaluate their impacts on Wuhan urban growth and produce the urban growth potential map. In addition to the population density data from 2015, other factor layers were acquired in 2017. We mapped all factors using ArcGIS 10.2 software and processed them by standardized methods. Thirty thousand sample points were extracted from the actual urban land-use change map (2009-2017) and the variable maps by a stratified random sampling method to construct the logistic regression equation. Table 2 presents the variables and their coefficients of BLR model. Figure 4 illustrates their spatial distributions and the final urban growth potential map.

Planning-Constrained Mechanism
Urban spatial planning, widely used in urban growth management, has an obvious impact on guiding and controlling the urban development [3,20]. Planning constraints have become an indispensable part of urban simulation. Different planning zones represent the different preferences of the government and planners for urban development and natural conservation. Previous studies established the planning constraint mechanism through binary planning constraints or preset planning zoning preferences [32,33], ignoring the cell-based trade-off between urban growth and conservation in the same planning zone [35]. In this study, by integrating the cell-based tradeoff between urban growth potential and conservation priority with zoning-based planning preference, a synthetical planning-constrained mechanism was defined as follows.
, , 1 where PCi is the planning constraints for cell i, ranging from 0 to 1. Ui and Ci are the value of urban growth potential and conservation priority, respectively. Zj represents the planning preference on urban growth in planning zone j, ranging from 0 to 1. Zj = 0 means planning zone j is completely unbuildable, Zj = 1 means planning zone j has no planning restrictions on urban development.
In past decades, Chinese local governments have issued a series of spatial planning policies to restrain the decrease in ecological and agricultural land and guide urban compact development, among which Overall Land Use Planning, Urban Master Planning, Primary Farmland Protection Planning and Ecological Protection Planning played extremely important roles [13,28,52]. By integrating these plans of the study area, we classified the study area into three categories of planning zones: the urban development zone (UDZ), ecological and agricultural conservation zone (EAZ) and completely prohibited undevelopable zone (PUZ). The UDZ refers to areas where cities and towns are intensively developed and constructed and can meet the needs of urban production and life. The value of the ZUDZ is set to 1, considering a higher preference on urban growth in the UDZ. The PUZ includes primary farmland protection area, mountains, water bodies and natural reserves wherein urban development is strictly prohibited. Thus, the value of ZPUZ is set to 0 in this study. The EAZ refers to agricultural and ecological areas besides the PUZ that need to retain their original appearance, strengthen natural conservation and ecological construction and restrict urban development. The value range of the ZEAZ is [0,1]. By adjusting the planning preference parameter in the EAZ, various planning constraints can be generated. Furthermore, considering the timeliness of planning, planning in different time stages should be applied to build the corresponding planning zones. In this study, the planning zones of 2009-2017 and 2017-2025 were established, respectively, and applied to simulate the urban growth process of corresponding phases.

Planning-Constrained CA Model
CA is a dynamic system with a finite set of state elements, which evolves in discrete time intervals based on transition rules [53,54]. CA-based models have been widely used in urban growth modeling and land-use simulation due to its simple and flexible structure. In general, the conceptual formula of the urban CA model can be described as follows [48]: where and denote the state of the cell i at the moment of t+1 and t, respectively; Ui is the urban growth potential of cell i, Ω is the neighborhood effect, Con is a constraint function, R is the random perturbation, and f represents the function of transition rules.
Specifically, the urban growth potential has been described in subsection 3.1.2. The neighborhood effect is defined and calculated based on the ratio of adjacent developed cells by using a regular n × n Moore neighborhood [25], which is expressed in Equation 6. Con can be replaced by PCi, which has been introduce in subsection 3.2. R represents the uncertainties and random perturbations in the urban land evolution process, which is calculated by Equation 7.
where c is the number of adjacent developed cells, γ is a random number ranging from 0 to 1, and α is a parameter for adjusting the stochastic magnitude, which was set to 1 to represent a small probability of randomness in this study. By combining the aforementioned planning-constrained mechanism with urban development potential, neighborhood effect and random disturbance, the final transition probability of urban CA model can be expressed as follows: Furthermore, the amount of urban land development was calculated by the Markov chain model, which has been widely and successfully used in the quantitative prediction of land-use change [3,55]. The Kappa index [56] and FoM coefficient [57] were used to calibrate the urban CA model and assess the simulation outcomes. By a pixel-by-pixel comparison between the simulation result and reference map, the Kappa index and FoM coefficient can be calculated as follows: where PO is the observation consistency, which refers to the ratio of the number of correctly simulated cells to the total number of cells; PC is the expected consistency; N is the total number of cells; a0 and a1 are the number of non-urban and urban cells, respectively; b0 and b1 are the number of cells simulated as non-urban land and urban land, respectively; A is the number of error cells observed as urban but simulated as persisted non-urban; B is the number of correctness cells observed and predicted as change; C is the number of error cells observed persistence but simulated to be urban land.

Multi-Scenario Design of Planning Constraints on Urban Growth
Multi-scenario analysis of planning constraints on the urban growth process has been conducted in this study. By adjusting the planning preference parameters of different planning zones, multiple planning-constrained urban growth scenarios can be derived. Considering the actual situation, the planning preference parameters of the UDZ and PUZ were preset to 1 and 0, respectively. To model the actual planning preference in the EAZ, 11 values, ZEAZ∈ {0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1}, were selected to generate different planning constraints. Among these values, 0 means corresponding planning zone is completely unbuildable, and 1 means corresponding planning zone has no planning restrictions on urban development. Then, multiple scenarios of urban land patterns were simulated for Wuhan in 2025 through integrating the multiple planning constraints into the urban CA model.

Urban Growth and Encroachment in 2009-2017
Based on the land-use maps in 2009 and 2017, we obtained the spatial pattern of urban growth in the WUDA from 2009 to 2017. Then, by superimposing it with the integrated urban spatial planning map, the spatial distribution of urban growth in different planning zones was produced, as shown in Figure 5. In addition, through further statistical analysis, the area and proportion of urban growth and encroachment on other lands in different planning zones are calculated, as shown in Table 3. During 2009-2017, the urban growth of the WUDA presented obvious infilling and edge sprawling characteristics, with a total area of .21 ha. At the same time, the expansion of urban land has led to the occupation of cultivated land, forest land, garden land, grassland and water areas. Among them, the occupied area of cultivated land accounts for 57.96% of the increased area of urban land, the occupied area of water area accounts for 23.69%, and the occupied areas of forest land, garden land and grassland account for 4.25%, 1.96% and 1.96%, respectively. These show that the urban growth in the WUDA during 2009-2017 is mainly at the cost of the reduction in cultivated land and water area. In addition, we found that the integrated urban spatial planning of Wuhan has an obvious guiding and restricting effect on the expansion of urban land at this stage, with 73.26% of the new urban land occurring in the UDZ. However, the urban land expansion beyond the preset UDZ should not be ignored, with 23.53% of the new urban land appearing in the EAZ, and even 3.21% in the PUZ.

Calibration and Validation of the Simulation Model
This paper analyzes the impact of urban spatial planning policies on urban growth during 2009-2017. By adjusting the planning preference coefficients of different planning zones in this time stage, we obtain the corresponding planning-constrained mechanism. Considering the practical significance of planning zones, the planning preference parameters of the UDZ and PUZ were preset to 1 and 0, respectively, while 11 values ranging from 0.0 to 1.0 on 0.1 intervals are used to generate different planning-constrained scenarios in the EAZ. By using the planning-constrained CA model described in Section 3.3, the urban growth process and spatial pattern from 2009 to 2017 under different planning constraints were simulated. Finally, the actual land-use map of 2017 is used to verify the simulation results. Figure 6 shows the changes in the Kappa index and FoM coefficient with the adjustment of planning constraints during the verification process. The results show that with the increase in planning preference coefficient of AEZ planning zone, Kappa index and FoM coefficient increase first and, then, decrease, which means that the simulation accuracy has the same change trend. The planning preference coefficient indicates the preference degree of different planning zones for urban growth relative to ecological or agricultural protection. The results show that the actual urban growth violates the urban spatial planning to a certain extent and pays more attention to urban development than ecological and agricultural protection. When the ZEAZ value is 0.7, the accuracy of the model is the highest. This shows that from 2009 to 2017, the urban growth of the WUDA is based on the development strategy that urban growth takes precedence over ecological and agriculture protection. Although this development strategy considers the trade-off relationship between urban development and ecological and agricultural protection to a certain extent, it obviously pays more attention to the economic benefits brought by urban development and largely ignores the ecosystem functional importance, ecological sensitivity and agricultural production suitability in the EAZ. This will cause higher risk of ecological and agricultural loss and further lead to the deterioration of ecological environment and the reduction in grain production.  Figure 7 shows the spatial distribution of urban growth probability under different planning constraints, which is calculated by the method described in Section 3.2. It can be found that with the adjustment of planning zoning preference parameters, the transition probability map of urban growth shows obvious differences. Overall, with the increase in ZAEZ, the area of land with high transition probability increases significantly, showing obvious influence of planning constraints. From the perspective of the planning zone, the land located in the UDZ shows significantly higher transition probability and remains unchanged with the increase in planning zoning preference parameters, mainly because the planned UDZ is often located around the existing urban land, and it is not affected by the spatial containment function of the planning zone. The land in the EAZ shows obvious spatial heterogeneity, which is mainly located at the urban fringe and suburban area. The farther away from the urban center, the lower the transition probability. At the same time, with the reduction in the planning zone preference parameters, the effect of planning constraints enhanced, resulting in lower transition probability. When the ZEAZ = 0, it means that the EAZ is subject to the strictest planning constraint, and the urban growth probability of all land located in this area is 0. When ZEZA = 1, it means that the EAZ is not affected by the planning constraint; thus, its transition probability only depends on its urban growth potential. When the ZEAZ < 1, it means that the transition probability of the EAZ area depends on the joint action of conservation priority and urban growth potential, and the land with higher urban growth potential and lower conservation priority is easier to convert to urban land. The transition probability of the land located in the PUZ is always 0, because it is subject to the strictest planning constraint.

Multi-Scenario Simulation of Urban Growth Based on Planning Constraints
Based on the urban growth potential map and the planning-constrained CA model introduced in Section 3.3, we simulate the spatial distribution of the WUDA urban growth in 2025 under different planning-constrained scenarios. We used the Markov chain model to predict the urban land demand in 2025, showing that the total area of urban land would be 1342.21 km 2 . Figure 8 illustrates the simulation results under different planning constraints. During 2017-2025, urban growth in the WUDA continued to occur around the existing urban land, showing obvious characteristics of edge and infilling expansion. At the same time, the urban growth in this stage is more distributed in the new town groups far from the urban center. The central urban area shows weak urban growth due to strict urban spatial planning policies and few land areas for urban growth. With continuous changes in the ZEAZ value, the urban spatial pattern evolved in a steady and gradient mode; however, a sharp transformation was observed with significant changes to the parameter. In comparison to the weak planning-constrained scenarios (when Z > 0.5), the strong planning-constrained scenarios (when Z ＜ 0.5) showed an obvious reduction in urban encroachment on farmland with high agricultural production suitability and ecological land with high ecological importance. This shows that the strong planning constraint scheme proposed in this paper can effectively avoid the occupation of cultivated land and ecological land with high conservation priority, which is conducive to the urban sustainable development. Table 4 presents the landscape indexes of urban growth simulation results under different planning constraints, including number of patch (NP), mean patch size (MPS), mean Euclidean nearest neighbor distance (ENN_MN), mean perimeter area ratio (PARA_MN) and aggregation index (AI). Overall, NP and PARA_MN decrease first and, then, increase, MPS, AI and ENN_MN show the opposite trend. When the ZEAZ ≤ 0.5, NP and PARA_MN values decrease gradually, MPS, AI and ENN_MN values increased gradually with the increase in the ZEAZ. When the ZEAZ > 0.5, with the increase in the ZEAZ, NP and PARA_MN values increase slowly, while MPS, AI and ENN_MN values decrease gradually. This indicates that when there is no obvious preference on urban growth or ecological land and cultivated land protection in the EAZ, the patch number of urban lands is the least, the patch shape is more regular, and the spatial distribution is the most concentrated. With the increase in preference on urban growth or ecological land and cultivated land protection in the EAZ, the patch number of urban lands increases, the patch shape is more irregular, and its spatial distribution is more fragmented. Further, the urban spatial distribution under weak planning constraints is more concentrated than that under strong planning constraints.   Based on the spatial structure of three main urban areas (Wuchang, Hankou and Hanyang) and six new city clusters (North, East, Southeast, South, Southwest and West) in the WUDA determined by the master plan of Wuhan, we have made zoning statistics on the simulation results of urban growth in the study area in 2025. Table 5 shows the area of urban growth and its proportion in the UDZ across different regions under multi-scenarios. With the increase in Z value in the EAZ, the proportion of urban growth in the UDZ will decrease in 2025, from 100% to 68.24%, which clearly shows the impact of planning constraints proposed in this paper on urban growth. In addition, we note that the South and Southwest new city clusters show more obvious changes, while the Hanyang main urban area has the least change. This shows that with the weakening of the impact of planning constraints, the urban growth of the South and Southwest new city cluster is more likely to exceed the UDZ preset in the urban spatial plans, while the urban growth of Hanyang tends to occur in the UDZ.

Discussion
CA-based models have been widely used to simulate urban growth and land-use change in past decades. Although there are many practices of using technologically sophisticated methods to construct the urban CA model, the incorporation of planning constraints into modeling has not been adequately addressed so far. By combining the cellbased trade-off between urban growth and conservation with a zoning-based planning implementation mechanism, we proposed a new planning-constrained CA model to simulate urban growth and encroachment under different planning scenarios and explore the potential influence of integrated urban spatial planning on urban growth process. In comparison to the earlier publications [15,[26][27][28], our research incorporated the synthetical urban planning constraints into the CA model, so as to make the predicted future urban scenario closer to reality and provide a useful evaluation tool for the urban planning implementation.
Taking the Wuhan Urban Development Area as a case study, the planning-constrained CA model was applied to simulate the spatial pattern of urban growth from 2009 to 2017 and predict the urban scenarios under different planning constraints in 2025. We find that a higher simulation accuracy can be achieved by considering planning constraints in the simulation. When the ZEAZ value is 0.7, the accuracy of the model is the highest, which indicates that the actual urban growth of the WUDA from 2009 to 2017 violates the urban spatial planning to a certain extent and pays more attention to urban development than ecological and agricultural protection. In future scenarios, with the weakening of planning constraints, urban growth tends to occupy more ecological and agricultural land with high conservation priority. With the increase in preference on urban growth or ecological land and cultivated land protection in the EAZ, the spatial distribution of future urban land becomes more fragmented. This shows that the planning constraint scheme proposed in this paper can effectively guide urban compact development and protect the ecological environment and food security. Furthermore, in the simulation process, the location and quantity of new urban land beyond the planned urban development area can be captured, which will provide certain early warning.
At present, the Chinese government and local governments are formulating a series of new spatial planning policies, one of which is to optimize the regional spatial pattern. It mainly includes the delineation of the three lines (ecological protection red line, permanent basic farmland protection red line and urban expansion boundary line), as well as the determination of planning zoning such as urban development area, rural development area, ecological protection area and agricultural protection area. Policymakers are faced with the problem of how to evaluate spatial planning policies and how to integrate ecological priority into urban spatial planning [32,33]. The methods and models proposed in this study can provide effective support.
However, there are several limitations in our work. First, the Markov model was used to predict the quantity of future urban land in this study, neglecting the impact of different policies and planning schemes, which should be addressed in further research. Second, this study mainly considers the natural and socio-economic driving factors to construct the transformation rules of the urban CA model, ignoring the impact of people's behaviors and decisions. Future research needs to focus on multi-objective human decision-making behaviors and their impact on the urban growth process. Finally, this study predicts the spatial pattern of urban growth based on the 30 m grid, ignoring that the actual land development process is based on the plot. Future research can consider building the conversion rules of urban CA model based on the land-use plot, so as to produce more accurate simulation results.

Conclusions
Spatial planning policies are important factors affecting urban development. How urban growth is impacted by urban spatial planning policies needs to be understood in order to better simulate the future urban patterns [28]. Many other publications have made significant contributions to the urban CA model and application [18]. However, so far, the scenario simulation considering the actual urban spatial planning policies has not been adequately addressed, which weakens the ability of CA to provide insight into the future urban spatial patterns and inform the sustainable development strategy [20]. In this paper, we developed a planning-constrained CA model to simulate urban growth and encroachment under different planning scenarios by integrating cellular automata and actual urban spatial planning. In the planning-constrained CA model, a synthetical planning-constrained mechanism was proposed by combing the cell-based tradeoff between urban growth potential and conservation priority with zoning-based planning preference. By adjusting the preference parameters of different planning zones, multiple planningconstrained scenarios can be generated. Taking the Wuhan Urban Development Area as a case study, the planning-constrained CA model was applied to simulate current (2017) and future (2025) urban scenarios. The main conclusions were: 1. The planning-constrained CA model demonstrated a higher simulation accuracy compared to the model without planning constraints.
preference on urban growth or ecological land and cultivated land protection in the EAZ, the future urban land pattern becomes more fragmented. 4. Location and quantity of new urban land beyond the planned urban development area can be captured in future urban scenarios, which will provide certain early warning.
The method proposed in this study can provide effective support for planners and decisionmakers to evaluate the impact of actual urban spatial planning policies on urban growth and generate urban future scenarios in which urban development and natural protection are coordinated, so as to realize the urban sustainable development.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.