Effect of zoning plans on urban land-use change: A multi-scenario simulation for supporting sustainable urban growth study the and patterns of Sustainable Cities and Society

Even though urban land-use change simulations provide useful information for decision makers, planning is generally weakly integrated into land-change modelling. However, the increasingly digitally available zoning data from statutory planning offers new opportunities. This study aims to reveal the potential effectiveness of statutory planning in terms of sustainable urban development by integrating zoning regulations in a multi-scenario simulation. Specifically, the gross floor area that can be built per parcel, as defined in the zoning plan, supports the allocation of varying degrees of urban densities. Using the CLUmondo modelling framework that couples cellular automata and multivariate logistic regression, we simulated urban growth in Madrid, Barcelona, Valencia, and Zaragoza Spanish Functional Urban Areas from 2012 to 2030 in four scenarios. The scenarios reflect the degree of planning intervention, ranging from high intervention to unrestricted develop- ment, while consider Spanish legislation and urban agenda 2030 sustainability goals. Simulations shows that by shifting growth to zones with urbanization projects almost 4200 ha of grassland and cropland could be saved from overbuilding, and 3800 ha by shifting it to zones without urbanization project. The simulation results provide critical information to support decision-makers and planners in revising plans and designing new plans to achieve sustainable urban development.


Introduction
Urbanization has increased all over the world generating intense changes in the spatial distribution, quantity and type of land-uses (Xu, Yang, Dong, Liu, & Qiu, 2016). Urban sprawl is a characteristic of this global urbanization (Huang, Huang, & Liu, 2019), constituting a challenge for sustainable urban development. Urban areas are not only hotspots for economic opportunities that could improve quality of life (Zhang, Ban, Liu, & Hu, 2011;Zhou, Dang, Sun, & Wang, 2020), but also associated with land degradation, reduction of natural areas, human health problems and socio-economic issues (Oliveira, Tobias, & Hersperger, 2018). The importance of urban areas for our future has been recognized in the UN's Sustainable Development Goals (SDGs), where two goals are directly associated to future urban development. In Europe, these goals have been subsequently transposed at the continental and national levels, including, among others, in Spain.
The 2030 Spanish Urban Agenda defines that spatial planning should make a rational use of soil, while urban sprawl should be avoided, revitalizing existing urban areas. Spain suffered one of the largest real estate bubbles in all of Europe before the last economic crisis in 2008. Legislative changes that allowed Spanish municipalities to obtain a significant increase in their income by creating new buildable land (Varela-Candamio, Rubiera Morollón, & Sedrakyan, 2019), unsustainable investment capital, and sometimes the gamble of the construction sector, generated urbanization pressures. These factors lead to an unprecedented expansion of urban land and sprawl (Díaz-Pacheco & García-Palomares, 2014). The definition of strict planning policies and their implementation are required to prevent future land degradation, keeping urban compactness in Spanish urban areas.
The simulation of urban land-use dynamics provides relevant information for planners (Karakus, Cerit, & Kavak, 2015) given that they are at the core of sustainable development goals. The scenario approach is appropriate for urban development as future urban key issues are unknown (Goodspeed, 2020). The scenario modelling is based on different assumptions, not only forecasting the most likely future, but describing alternative urban land-use changes that can assist on urban expansion management and plan development. Explorative scenarios offer the possibility to study the magnitude and spatial patterns of likely land-use changes (Verburg, Eickhout, & Meijl, 2008) while providing the possibility to evaluate policy proposals (Schulp, Nabuurs, & Verburg, 2008). The integration of planning policies in scenario land-use modelling plays a relevant role in urban expansion, which can help improving simulations accuracy (Onsted & Chowdhury, 2014;Verburg & Overmars, 2009). Even so, this is often disregarded Zhou et al., 2020).
The definition of likely scenarios, and subsequent implementation in a modelling framework, should account for driving factors and actors involved in land-use change at local and regional scales (Huang et al., 2019;Liang et al., 2018). In this sense, CLUMondo modelling framework was chosen to integrate the local "bottom-up" effects and the regional "top-down" perspective. Local effects refer to the relationship between land-use patterns and spatial driving factors, while regional effects are related to future socio-economic demands that define urban growth rates (Van Asselen & Verburg, 2013;Verburg & Overmars, 2009). The integration of both local dynamics and future demands better characterize future urban development patterns (Liang et al., 2018;Zhou et al., 2020). CLUMondo include cellular automata (CA) and suitability models that have been selected in recent years to model urban expansion at micro-scale (Asgarian, Soffianian, Pourmanafi, & Bagheri, 2018;Dadashpoor, Azizi, & Moghadasi, 2019;Saganeiti, Mustafà, Teller, and Murgante (2020); Zhou et al., 2020) urban dynamics or delimit urban growth boundaries (Huang et al., 2019;Liang et al., 2018). The expansion of urban spatial patterns are frequently restricted by protected land, cultivated land, or areas with cultural or planning constraints (Zhou et al., 2020), which should be integrated in urban growth simulations. Planning can be implemented in modelling using either hard or gradual restrictions . Different studies have addressed the inclusion of spatial planning policies in land-use change modelling (Geneletti, 2013;Liang et al., 2018;Lin & Li, 2019;Poelmans & Van Rompaey, 2010;Zhou et al., 2020). However, the integration of zoning constraints is generally established as hard or Boolean-based restrictions (e.g., whether urban development is allowed or not), while not accounting for the spatial heterogeneity or gradual characteristics within planning zones (e.g., whether planning regulations allow low, medium or high urban density). Little research has focused on the use of softer constrains based on weighs or the definition of land-use intensities within defined zoning areas, though these could improve real patterns simulations in urban areas (Li & Yeh, 2000;Zhou et al., 2020). In this sense, differentiated zoning regulations that are related with land-use intensity, such as population and floor area (van Vliet, Verburg, Grȃdinaru, & Hersperger, 2019) are analysed. Spanish General Zoning plans are suitable to explore these aspects within land-change modelling because they define, at the municipal level, three main land classes: urban, buildable and rural land (García-Ayllón, 2018;SIU, 2021).
In light of the above considerations, the overall aim of the present study is to model urban land-use changes using a multi-scenario approach that integrates digitized zoning plans for the Functional Urban Areas (FUAs) of Madrid, Barcelona, Valencia, and Zaragoza. The following specific objectives are addressed: i) to analyse the role of planning by defining three future scenarios that integrate digitized zoning plans and one scenario that assumes almost no planning intervention; ii) to introduce zoning constraints that reflect different degrees of urban densities; iii) to generate a transferable spatially-explicit modelling framework to integrate planning into land-use change simulations.

Study area
The study area encompasses four Functional Urban Areas (FUAs) in Spain (Fig. 1), each defined by the city and its commuting zone according to the Digital Spanish Urban Atlas and Eurostat from the  The Madrid FUA is located in the centre of Spain, has a population of 6,120,254 inhabitants (2018, Spanish Statistical Office, INE), and is divided into 52 municipalities (2,887 km 2 ). Since the 1990s, urban growth has increased significantly, being one of the European Union's urban growth hotspots. This growth has occurred predominately in six rings of suburban municipalities (Díaz-Pacheco & García-Palomares, 2014), mostly in a sprawling pattern due to the increase in semi-detached and terraced houses in small municipalities (Moliní & Salgado, 2012).
The Barcelona FUA is located on the north-east coast of Spain, has a population of 5,108,383 inhabitants (2018, INE), and is divided into 124 municipalities (2,341 km 2 ). It is characterized by a polycentric metropolitan system with a compact centre, two metropolitan rings, and 7 sub-centres (Catalán, Saurí, &Serra, 2008, andTombolini et al., 2015). The urban growth trend has been associated with urban spillover in the outer part of the metropolitan region, an increase in retail activities along the urban fringe, and the presence of detached family houses (Tombolini et al., 2015). According to INE, an annual population growth rate of 0.7 % is projected until 2033, which would potentially lead to urbanization in fringe areas (Serra, Vera, & Tulla, 2014).
Valencia FUA is located on the centre-east coast of Spain, and has a population of 1,552,783 million inhabitants (2018, INE), half of whom live in Valencia central city. The Valencia FUA is divided into 45 municipalities with an area of 634 km 2 . Although the population has grown slowly at an annual rate of 0.7 % since 1990 (INE), the GFA (gross floor area) has increased with an annual growth rate of 3.6 %. The urbanization process is characterized by the progressive growth of Valencia city and growth along two main axes (north-south, and coast-interior) as the most populated municipalities outgrow the first urban areas.
The Zaragoza FUA is located in the north-east of Spain, between the two main Spanish FUAs, and has a population of 744,862 inhabitants (2018, INE). Although the Zaragoza FUA is divided into 15 municipalities with an area of 2,205 km 2 , almost 90 % of the population live within the Zaragoza municipality. The metropolitan area has grown mainly in the last two decades due to lower land prices in municipalities close to Zaragoza city (De Miguel, 2014). While the GFA has increased since 1990 at an annual growth rate of 3.87 %, the population has grown only at a rate of 0.8 % a year (INE). Thus GFA growth has strongly exceeded the demand for living space (De Miguel, 2014).

Multi-scenario approach for land-use simulation
The modelling of future land-use transitions was carried out using a scenario-based approach in CLUMondo. The approach considers key local-scale driving factors, future demands to account for regional effects and zoning plans constraints to simulate the magnitude and location of future land-uses. Simulations were performed for the year 2030, a time frame associated with the 2030 Spanish urban agenda and the Sustainable Development Goals (SDGs). CLUMondo requires four main inputs: the total demand per year associated with each specific scenario (Sections 2.3 and 2.4), the conversion settings, defined as a conversion matrix, to promote or restrict specific transitions between land-use types (Section 2.4), the initial state of the study area, which refers to the EUA 2012 land-use base map (Section 2.5.1), and location probability maps for each land-use type (Section 2.6) (Ornetsmüller, Verburg, & Heinimann, 2016) (Fig. 2). GFA and population demands were defined for each scenario based on urban plans and population growth scenarios for Spain (at the provincial level). Conversions were specified in a conversion matrix (Appendix A Tables A8 and A9) while spatial restrictions were based on Spanish General zoning plans and integrated in the conversion matrix (Appendix A, Tables A8 and A9, and Fig. A5). The zoning plans selected for the analysis constitute a harmonized and updated version of the current municipal zoning plans of Spain developed by the SIU "Sistema de Información Urbana" after the Land and Urban Rehabilitation law of 2015. The zoning plans reflect the medium to long-term future of urban areas. With minor updates within the plans, it is possible to adapt the plans to new needs, sometimes exploited by the increasing power of interest groups (Janin Rivolin, 2017). Though, these updates can only make small changes at a very high level of detail that are not expected to affect the boundaries and land-uses used in the simulations. These plans define urban, development, and rural land classes at municipal level. Urban land refers to urbanized areas with access to essential utilities. Development land includes two sub-classes: buildable areas with an approved partial plan and an urbanization project (hereinafter: buildable with urbanization project), and buildable areas for which no urbanization project has yet been defined (hereinafter: buildable without urbanization project). Rural land includes croplands and forests with building only allowed for agriculture and livestock. The initial state of the study area was defined by the European Urban Atlas (EUA) in 2012 and associated statistics for each FUA. Local-scale drivers of land-use change were selected to compute probability maps for each land-use type (2.5.2). CLUMondo modelling framework was selected to allocate future land-uses. Simulations were carried out from 2012 to 2018 to validate simulation performance, and subsequently simulations were performed for the year 2030 (Section 2.7 and Fig. 2).

Scenario storylines
Four future land-use demands scenarios were defined for the FUAs to characterize uncertainty in the development of these drivers. Storylines were created considering probable development scenarios related to zoning plans, current Spanish legislation and sustainability goals defined along two axes: a high market-oriented vs. high planningintervention axis, and an axis of short-term economic growth vs. longterm sustainable growth (Fig. 3).
The sustainable development scenario (S1) is characterized by low GFA growth that is limited to areas that are currently under development according to zoning plans. Capital investment is associated with local companies while the growth in population is slow. There is high ecological awareness which leads to avoiding new built-up areas being located in protected areas or potentially floodable areas, even if these are considered buildable by zoning plans. There is an emphasis on proximity food production to urban areas and self-sufficiency, which fosters cropland as protected areas.
The business-as-usual scenario (S2) is characterized by medium GFA growth in the range of on-going trends. Buildable land is limited to buildable areas with urbanization project included in zoning plans, with awareness of protected reserves and floodable areas. Zoning defines specific or sectored plans with the maximum buildable GFA. Capital investment is associated with local and national companies, and population growth is in line with by provincial trends (INE).
The strong development scenario (S3) is characterized by high GFA growth rates. Growth is restricted to buildable areas without urbanization project designated in zoning plans. Modifications to current zoning plans are assumed. International capital investment is expected, while population growth is slightly higher than provincial trends as defined by INE. Low concern for sustainability is expected, but protected areas are respected.
The unrestricted development scenario (S4) prioritizes a high degree of market liberalization characterized by high GFA growth that surpasses population demands. S4 follows a rapid economic growth pattern with almost no planning intervention. S4 thus foresees a "construction bubble" that emulates previous liberalizing legislations that provide capital gains to municipalities by fostering excessive urban growth (Diaz Pacheco & Hewitt, 2010). Strong international capital investment is taken on while population growth is high. Although there is no environmental awareness, natural reserves defined at a European level are respected. The storylines do not consider a scenario with high market orientation and long term sustainability, as generally companies try to have higher profits at short term. A scenario related to high planning intervention a fast growth was neither considered, as planning should help prevent land degradation at short term.

Scenario quantification
Scenario quantification was associated with two model settings: the demand for goods and services, and spatial restrictions to promote or restrict specific transitions from one land-use type to another. The four defined scenarios consider two demands that influence urban growth and urban land-use density, namely, GFA and population. Scenario quantification was addressed in three steps: i) quantification of the initial GFA and population supply, ii) definition of spatial restrictions with three urban densities, iii) quantification of future GFA and population demands.
The quantification of the initial GFA and population supply was determined for each land-use type and FUA (Appendix A, Table A6). The GFA per building derived from the Spanish cadastre, and population density from the 2011 census (INE) were used to compute per land-use initial demands.
The definition of spatial restrictions was based on zoning plans for S1-S3, while S4 assumes almost no planning intervention in an urban sprawl context characterized by medium housing densities at the urban fringe. Zoning plans were used to delimit initial land-use areas, enabling transformations to urban or industrial in S1-S3 (Appendix A, Tables A8 and A9). Zoning plans include the total GFA in planned development areas, which enabled the definition of three different urban densities: low, medium and high. The thresholds between urban densities were computed based on the initial GFA supply per urban land-use from the EUA in 2012. The average GFA within "continuous urban" and "dense urban" land-uses defined high-density urban areas, average GFA within "medium dense urban" land-use determined medium-density urban areas, and average GFA within "low dense urban" and "very low dense urban" land-uses defined low-density urban areas.
These thresholds were applied to define low, medium and high densities for developing areas in S1 and S2. The development areas for S3 defined by zoning plans required the creation of urbanization projects that define GFA (buildable without urbanization project). The densities from the closest development areas were assigned using GIS techniques and verified by visual interpretation using high-resolution ortophotos from the PNOA, using positive spatial autocorrelation principles to define S3 densities (Dungan et al., 2002). Although S4 assumes almost no planning intervention in an urban sprawl context, medium  and low urban density areas were defined based on previous locations of urbanized areas for each FUA (Zhao, Weng, & Hersperger, 2020). The definition of medium and low urban densities for S4 was carried out using adaptive buffers for each FUA, considering the area that new buildable land represents with respect to the total area of the FUA. A linear interpolation of the percentage of medium and high density buildable land in S1, S2 and S3 was carried out to determine the maximum area with medium density for S4. Low densities in S4 were allowed for these non-medium density areas (Appendix A, Fig. A5). The quantification of future GFA demands was based on historical trends and literature findings, while future population demand was based on INE population projections and historical trends ( Table 1). The quantification of GFA demand for S1 includes the 2009− 2019 trend after the largescale construction estate bubble in 2008. S2 accounts for the 1990− 2019 historical trend characterized by legislative and social changes that increased urban development (Diaz Pacheco & Hewitt, 2010) after Spanish EU membership. The growth of urbanization in Spain is linked to the economic growth that began in 1975 (Moliní & Salgado, 2012), which generated a relevant increase in population in less than twenty years (Marroquin, Morollon, & Rivero, 2013). Consequently, the historical trend between 1975 and 2019 was taken into consideration for S3. S4, characterized by the highest growth rates in terms of GFA and population, accounts for the 2000− 2008 historical trend, a period with unprecedented urban expansion and city sprawl previous to the economic crisis (Díaz-Pacheco & García-Palomares, 2014). The quantification of population demands for S1 and S2 was determined by INE population projections for 2033. The population demand for S3 and S4 was computed using historical trends, taking into consideration the same historical periods used for GFA demand in S3 (1975− 2019) and S4 (2000− 2008).

Data sources 2.5.1. Land use data
The European Urban Atlas (EUA), developed by the European Environmental Agency, provides high-resolution, reliable and intercomparable LUC data. Data from the 2012 and 2018 EUA were selected. The EUA identifies 27 LUC classes distributed within 5 thematic groups (artificial, agricultural, natural and semi-natural, wetlands and water). The EUA LUC data was transformed to raster format at 30-m resolution, considering the central pixel value. The EUA classes were reclassified into 15 classes (Appendix A, Table A1). Ten classes describe artificial surfaces, one class represents agricultural areas, three classes refer to natural and semi-natural areas, and one class integrates wetlands and water.

Driving factors
Three groups of spatial drivers commonly used to model future landuse allocation at regional scales were determined: environmental, socioeconomic and proximity (Pazúr & Bolliger, 2017;Zhou et al., 2020) ( Table 2). The drivers were computed in raster data format at 30-m resolution to match the EUA data. The transformation of shape data format to raster format and subsequent resampling of the data were carried out to generate model input drivers.
Eight environmental drivers were considered: precipitation, solar radiation, geology, distance to flood plain area, protected areas, elevation, slope and aspect. Average precipitation from climate stations across Spain was provided by the Ministry of Agriculture and Fisheries, Food and Environment. The global horizontal irradiation product, provided by the Global solar atlas was used to describe average solar radiation. The Geology map describes different lithology classes within each FUA (Spanish Geological and Mining Institute). The Euclidean distance to floodable areas, defined as areas with a 500-year return period, was computed to determine the distance to flood plain areas. Protected areas with high protection grades, defined at local, regional, national and European levels, were identified, generating a dichotomous variable. Slope and aspect were based on the digital elevation model provided by the Spanish National Plan for Aerial Orthophotography (PNOA).
Five socioeconomic drivers were included to characterize the demographic and labour market profile. Population density was derived from the Global Human Settlement Population layer of 2015, while for future scenarios, we considered existing projections of population change (INE). Household density and the time required to reach the workplace was derived from 2011 Spanish census. The change in the number of companies and employed persons from 2012 to the present was provided by INE and computed using monthly reports from the Spanish General Treasury of Social Security.
Proximity drivers were created using Euclidean distance, which constitutes a commonly used approach (Liang et al., 2018;Zhang et al., 2011). Distances to major and urban roads were selected to determine proximity to transit lines. Distance to public transport stops, to bus lines and to train lines were computed using transport city datasets (Open Street Maps). Distance to green areas was included due to the association of parks and green areas with a higher quality of life. The distance to buildable land with urbanization project, and distance to buildable land without urbanization project as well as distance to the centre of the main city in an urban region, and distance to secondary settlement centres were computed as Euclidean distances.

Probability maps
The location probability of each land-use type was predicted by relating the current distribution of land-uses with a set of driving factors (Ornetsmüller et al., 2016). A multivariate logistic regression model (MLR) was used to predict, at pixel level and for each land-use, the probability of occurrence of a land use (Price et al., 2015). Prior to modelling, random samples were taken considering a minimum distance of two pixels between each sample to reduce spatial autocorrelation (Dungan et al., 2002). A balanced number of samples, in terms of presences and absences, was determined for each land-use type considering the area that represents the total FUA area. For frequent land-uses, a maximum of 3,000 land-use presences and 3,000 absences were considered. MLR was developed for the 15 sampling sets associated to 15 different land-uses in each FUA. The model for each land-use type included a set of best fitted drivers according to AIC value and variable significance. Model performance was validated using the Area under the Receiver Operating Curve (AUC) (Pontius & Schneider, 2001). AUC values range between 0.5 and 1, indicating random allocation and perfect fit, respectively. The probability models were the same across the scenarios, but different for each FUA.

Allocation and validation
The CLUMondo land system model was applied at 30-m resolution to simulate land changes based on different demands at the same time (Van Asselen & Verburg, 2013). Land-use for a certain location simulated by CLUMondo takes the following into account: local probability of each land use based on the MLR, GFA and population demands associated to each scenario and FUA, zoning restrictions, and specific LUC conversion settings. Zoning is integrated as spatial layers that restrict urban and industrial development for each scenario and FUA. Concretely, up to four raster layers where defined to account for the different intensities of new buildable land (see 2.4. Scenario quantification). These layers were included in the conversion matrix (see Appendix A, Tables A8 and A9) to restrict spatially the conversion of specific land-uses (see Appendix A, Fig. A5). CLUMondo only requires one initial land-use map for simulating future land-uses (Ornetsmüller et al., 2016). Consequently, historical land-use maps were not considered for model calibration. Firstly, we run the simulations from 2012, using the EUA 2012 as base map, to predict the land-uses in 2018. Secondly, an independent validation of CLU-Mondo simulations was performed considering the assessment of Fig. 4. Area of land-use change (km 2 ) by land-use type for scenarios S1 -S4 in FUAs, for 2012 to 2030. For explanation of land use classes, see legend in Fig. 1.

Table 3
Summary of location and pattern accuracy from simulated land-uses. Values are expressed in percentage. AF stands for aggregation factor, WPP for well predicted performance, OP for overall model performance, PA for Producer's accuracy, UA for User's accuracy, and FIS for fuzzy inference system.  location accuracy (Diogo & Koomen, 2011;Pontius et al., 2018) and pattern accuracy (Power, Simms, & White, 2001 (Pontius et al., 2008). Model performance was summarized using hits (correct pixels based on observed change predicted as a change), false alarms (erroneous pixels based on observed persistence predicted as change), misses (erroneous pixels based on observed change predicted as persistence) and correct rejections (correct pixels based on observed persistence predicted as persistence) (Pontius et al., 2018). The producer's accuracy, user's accuracy (Pontius et al., 2008), well-predicted performance and overall performance were computed (Diogo & Koomen, 2011). Model performance statistics were computed at three aggregation levels (1-3). Pattern accuracy validation does not consider small displacements in pixels to be disagreement even though the land use simulated patterns are essentially the same (Power et al., 2001). The pattern validation compares the observed map of 2018 'EUA 2018 LUC map' and simulated map 'Simulated 2018 LUC map' using the fuzzy inference system (FIS) statistic, which was computed with Map comparison kit.

Probability model performance and allocation validation
The number of drivers selected to assess land-use change using the MLR regression for the four FUAs ranged from 2 to 9 (Appendix A, Tables A2-A5), with population density, distance to roads, distance to centres, elevation and slope drivers those selected most often. Population density was the main driver for urban land-uses followed by distance to roads, while slope was the main driver for natural and agricultural land-uses. MLR performance showed generally good power for all land-use types and the four FUAs, with an average AUC value of 0.8 (Appendix A, Tables A2-A5, and Figs. A1-A4). Models associated with urban land-use types showed an average value of 0.9 AUC. Similar excellent MLR performance was found when modelling airports, green urban areas and water. MLR probability of industrial and forest landuses had good performance on average, while lower AUC values were found when modelling roads, construction areas, croplands, pastures, and natural grassland, indicating fair model performance.
The independent validation of the location accuracy at pixel level showed good overall performance, with 93 % correctly allocated landuse (Table 3, and Appendix A Table A10). Although the model correctly simulated the persistence of land use, simulations showed low levels of agreement between the observed and simulated land changes at pixel scale. These low levels of agreement may be related to the small net change between 2012 and 2018 (Pontius et al., 2008). The pattern accuracy validation determined that simulations properly predicted land-use and change patterns similarly, with an average value of 79.5 % for the FUAs.

Future land-uses under defined scenarios
The S1 sustainable scenario shows very little increase in urban land uses, and these changes take place mainly at the expense of construction areas, generally characterized by the presence of bare soil or scarce vegetation (Fig. 4). The scenario also simulates an increase in cropland and forested areas at the expense of pastures and natural grassland. S2 shows greater transition to urban land uses through the conversion of construction areas, natural grassland and pastures with values ranging from 16 km 2 in Zaragoza up to 165 km 2 in Madrid. We should stress that specific change values vary for different FUAs according to total FUA size (see Section 2.1 and Fig. 4). These trends are emphasized in S3, which is associated with a strong development that reaches up to a 190 km 2 increase in the case of Madrid. The growth of urban land uses shows a predominance of medium urban densities at the fringes or in new development areas (Figs. 5 and 6). S4 is characterized by low planning intervention and a high degree of market liberalization. The simulated maps show, as expected, the most extreme land-use changes, in particular a strong trend toward sprawl with predominantly medium-and low-density urbanized areas creating new areas of development far from existing centres.
Generally, the urbanization process mostly occurs at the expense of natural grassland, construction areas and pastures in Madrid and Zaragoza, while croplands are also affected in Barcelona and Valencia in scenarios S2 to S4. Simulations predict an increase in forested areas in Valencia and Zaragoza at specific hotspots in the south of Valencia and northwest of the Zaragoza FUA (Fig. 6). Industrial areas show a small decrease in S2 and S3, with up to 15 km 2 (in Madrid) being converted from industrial areas to other uses. This may be due to the regeneration of brownfields. Scenarios S1 to S3 ensure dense and medium-dense urban densities with low presence of low densities in Valencia and Zaragoza. A similar presence of low to high densities is observed in Barcelona, which could be associated with urban saturation (Serra et al., 2014), while a higher prevalence of low densities is observed in Madrid in all scenarios.
Urbanization is mainly expected to take place in the external rings in Madrid and within coastal municipalities, rings and sub-centres such as Terrasa, Rubí, Sabadell (locations of towns and labels are included in S2 in Figs. 5 and 6) or those close to the coast such as Mataró and Vilanova i la Geltrú in the case of Barcelona (Fig. 5). In contrast, in Valencia and especially Zaragoza, urban growth is expected at the fringes of the main city (Fig. 6). Urban growth in Madrid is mostly simulated in the first rings located south and north east of Madrid city, as the northwest area is protected. The greatest increase is still modelled at the border of Madrid municipality close to Alcorcón, between Getafe and Coslada, and south of Alcobendas. In Valencia, new built-up areas are found in the main city, the municipalities located in the north and south of Valencia city, and municipalities in the interior, including, Torrent, San Antonio de Benágeber, Paterna, or el Puig de Santa María. Especially the surroundings of Valencia tend to the creation of a conurbation. Zaragoza shows clear development of new urban areas along the previous development axes associated with the Ebro (west-east), Gallego (north-south) and Huerva (south-north) rivers. Zaragoza remains the most relevant urban area, while municipalities in the adjacent south (e.g., Cuarte and Cadrete) also show a significant increase in built-up area. In summary, the simulated patterns show a growth of urban land-uses at the external rings and sub-centres in Madrid and Barcelona, while urban growth is expected at the fringes of urban centres in Valencia and Zaragoza.

Multi-scenario comparison
The comparison of scenarios shows that the integration of spatial planning constraints in land-use models plays a relevant role in guiding sustainable urban development in Spanish FUAs. The magnitude of urban area increase is very different between the simulated scenarios. Under S1, the increase of medium-dense and dense urban land-uses dominates. Some neighbourhoods show densification of existing urban areas as for example in the urban core of Valencia and Zaragoza (De Miguel, 2014), over Terrasa and Sabadell sub-centers in Barcelona or Alcobendas and Coslada in Madrid. Comparing S2 and S3 we found that several neighbourhoods have the same development trend because in both scenarios the same zoning constraints and maximum gross floor area settings restrict the land-use simulated changes. However, the urban expansion rate is higher in S3 as can be observed in Zaragoza where simulated land uses tend to create a ring around the main city. A similar pattern is visible in Madrid, that growths within external rings, with higher presence of planned low densities that increase urban land patches (Moliní & Salgado, 2012). These spatial patterns match the sprawling trend observed during the last decades in accordance with Díaz-Pacheco and García-Palomares (2014). S4 scenario shows the effect of the absence of strict land-use planning, i.e. sprawling urbanization patterns and a loss of natural and agricultural lands. S4 simulates future urban land in the four FUAs as a disorderly and non-sustainable expansion. Similar patterns also have been reported in previous studies that compare scenarios with and without planning (e.g. Zhou et al., 2020). Under S4 the allocation of industrial areas is influenced by transportation network location in the case of Valencia. Furthermore, the strong competition between industrial and urban land-uses generate unusual patterns as the increase of industrial areas close to cities. Some of the new developments simulated in S4 are found in isolated or small urban areas. These could be associated, as in the case of the central Ebro valley of Zaragoza FUA, with the existence of irregular buildings (buildings erected without valid permits). Similar urban growth patterns occurred in the 1970s in Barcelona FUA (Catalán et al., 2008), highlighting the importance of spatial planning to avoid discontinuous urban development.
Overall, spatial patterns show that medium and higher urban densities, which are associated with lower land consumption and compact urban forms that avoid the increase in sprawl and land degradation, could be reached when applying planning instruments to create sustainable urban environments. Explorative scenarios, especially when supported by cartographic outputs, assist policymakers to evaluate at low-cost the magnitude and spatial patterns of alternative development schemes and policy proposals (Zhang et al., 2011).

Driving factors and simulations
The land-use distribution patterns are generally influenced by socioeconomic drivers, distance to networks and terrain characteristics in agreement with previous studies (Kim, Newman, & Güneralp, 2020;Zhou et al., 2020). Probability models included between 2 up to 9 drivers out of 23 computed (Appendix A, Tables A2-A5). Population was the major social driving force, being correlated with urban land needed. Despite the fact that the increase of urban area in Spain has surpassed population needs during the last decades, shaped by legislative changes, competitiveness and private and political forces (Rivolin & Faludi, 2005;Varela-Candamio et al., 2019), population is still a key driving factor in depicting future land-use patterns in accordance with Kim et al. (2020). In this sense, the population is expected to continue growing for the analysed FUAs, as these constitute economic poles of attraction within the Spanish territory. Distance to transportation network is related to urban growth development as areas with higher accessibility are more easily urbanized (Zhou et al., 2020). A clear example is the growth of Valencia along the north-south and coast-interior axes related with the main highways (named AP-7 and A3). Furthermore, topographic drivers such as slope or elevation had a relevant importance to explain not only the urban development but also the location of natural and agricultural land-uses. In this sense, areas with smooth slopes favour urbanization (Price et al., 2015), and agricultural intensification (Pazúr & Bolliger, 2017) while steep slopes could limit the growth of forests land-uses. Probability model performance for individual land-uses was good with an average AUC value of 0.84, in line with Ornetsmüller et al. (2016) and Pazúr and Bolliger (2017). Overall model performance ranged from 93 to 99 %, while lower allocation agreement values were found in accordance with Pontius et al. (2008) and Varga, Pontius, Singh, and Szabó (2019). The complementary validation at patch level made it possible to determine similarity patterns between simulated and reference maps (Power et al., 2001), showing good agreement.
Our simulations show that integrating strict measures from zoning plans help maintaining compactness and urban sustainability, as was recently proposed by Zhou et al. (2020). Despite spatial plans are sometimes subject to some actualizations, that could be influenced by the increase power of interest groups or changes in planning governance (Janin Rivolin, 2017), the integration of spatial plans improves output reliability in land-change simulations (Li & Yeh, 2000). In this sense, the increasing laws to manage the use of land to be urbanized, as in the Spanish case the Land and Urban Rehabilitation law of 2015, provides new legal frameworks to restrict changes in boundaries and land-uses planned strengthening simulations validity.
The simulations results defined accurately potential patches of new development areas, though, the simulations present some limitations to allocate the exact pixel of change. These allocation disagreements at pixel level were previously reported by Pontius et al. (2008), 2018, and could be related with the low net change between the start and subsequent reference land-use maps and the short time simulation period (six years). In this sense, the use of EUA from 2006 was not considered to increase the time between the start of the simulation and the validation because a lower number of urban land-use classes are provided in the first EUA version. Despite harmonization of the datasets constitute an issue, finding consistent data for two historical periods is difficult (Ornetsmüller et al., 2016), and the EUA datasets constitutes a promising data source that can be used for validating land change simulations (Ornetsmüller et al., 2016) while enabling transferability of the work procedure for European urban areas.

Future research
Urban land-use change was predicted using CLUMondo framework based on a multi-scenario simulation that integrates zoning constrains. In the future, giving the increasing availability of digital plan data (Hersperger, Grȃdinaru, & Siedentop, 2020;Fertner et al., 2019), the use of this data need to be explored to model urban growth at the national level and for other European countries to better assess sustainable urban management at larger scales keeping urban compactness. Further work should be devoted to integrating other planning instruments, such as strategic spatial planning. In this sense, the use of methodologies to visualize planning intentions (Palka, Grȃdinaru, Jørgensen, & Hersperger, 2018) can help to incorporate strategic planning visions into land-change modelling. Furthermore, the increasing complexity of planning should embrace not only the integration of spatial plans as hard restrictions in modelling, but planning information related with land-use intensity that account for the spatial heterogeneity of planning zones. This study constitutes a first attempt to simulate changes integrating the spatial heterogeneity of zoning constraints based on maximum floor area derived from the plans. Land-use and urban morphology are key aspects of the built environment, which were traditionally analysed as a horizontal dimension while the exploration of urban 3D development is gaining momentum. Accordingly, the exploration of modelling frameworks that links land-use change with 3D urban morphology characteristics based on spatial planning constrains is a suitable research line to address sustainable urban development.

Conclusions
This study reveals the importance of integrating urban planning in urban land-change simulations to provide insights into future scenarios of Spanish urban area development. The CLUmondo modelling framework was used to simulate future land-uses for Madrid, Barcelona, Valencia and Zaragoza FUAs from 2012 to 2030 considering four scenarios. An innovative aspect of the study is the integration of digital plan data not as hard restrictions, but considering the heterogeneity within zoning areas along three intensity classes (low, medium and high). This softer integration of planning improved simulated real-patterns by adding the urban intensities based on maximum gross floor area defined in zoning plans. The 2030 simulations show that the greatest change occurs at the fringes of the main cities in the cases of Valencia and Zaragoza. In Madrid FUA, most changes are concentrated in the external rings close to Alcorcón, between Getafe and Coslada, and south of Alcobendas, while the greatest urbanization growth in Barcelona is located in sub-centers such as Terrasa, Rubí, and Sabadell. The simulated urban growth is inversely proportional to the intensity and constraints defined in the scenarios. This relationship leads to drastic changes in land-use patterns in the case of the unrestricted development scenario, which shows a disorderly increase of low density urbanization. Simulations shows that restricting urban growth to buildable areas with urbanization projects keep urban compactness and save almost 4200 ha of land-taken from grassland and cropland respect to unrestricted development. Despite populations is expected to continue growing in the FUAs at low rates, scenarios based on strict land-use planning constitutes a relevant output for planners to prevent non-sustainable trends and urbanization do not overpass populations demands, as previously occurred in the estate bubble of 2008. The proposed land-change multiscenario work procedure is transferable and could be implemented worldwide to account for different urban demands. Its application to Spain and European FUAs is straightforward, given the existence of the Urban Atlas, and considering the increasing availability of digital zoning plans that provides information about boundaries and intensity development. Resulting scenarios provide relevant information to support decision-makers and planners towards the spatial identification of vulnerable areas to urbanize, the definition of urban intensities according to population demands, the revision of existing zoning plans, and the design of new plans to develop regulations that restrict disorderly urbanization with the aim of reaching urban sustainability goals.

Declaration of Competing Interest
The authors declare there are no conflict of interest.

Acknowledgments
This research was funded by the Swiss National Science Foundation (Consolidator Grant number BSCGIO 157789). We would like to thank Dr. Vasco Diogo for his invaluable help during the simulation validation and for further suggestions. The authors thank the anonymous reviewers for their insightful comments and critical reviews, which helped to improve our manuscript.

Table A2
Logistic regression independent drivers by land use at Barcelona urban region.

Table A3
Logistic regression independent drivers by land use at Madrid urban region.

Table A5
Logistic regression independent drivers by land use at Zaragoza urban region.

Table A7
Land-use conversion priority matrix (lusconv) homogenised for the different urban regions and scenarios. Values indicate the order of supply by land use varying from the highest rank (5) to the lowest rank (1). The value 0 indicates that no supply is provided by a specific land-use to supply the demand.

Table A8
Allowed conversions between land-uses (allow) standardized for Scenarios 1-3. The matrix determines whether the conversion from one land-use to another is allowed (1), restricted (0), allowed after 10 timesteps (110) or allowed in specific locations (17 to 20). The values 17, 18 and 19 refer to areas in which low, medium or high urban density could be allocated according to zoning plans, respectively. The value 20 refers to areas to which industrial land-use could be allocated.

Table A9
Allowed conversions between land-uses (allow) standardized for Scenarios 4. The matrix determines whether the conversion from one land-use to another is allowed (1), restricted (0), allowed after 10 timesteps (110) or allowed in specific locations (17 to 18). The values 17, 18 refer to areas in which low or medium urban density could be allocated according to zoning plans, respectively.