Predicting Land Use Changes in Philadelphia Following Green Infrastructure Policies

Urbanization is a rapid global trend, leading to consequences such as urban heat islands and local flooding. Imminent climate change is predicted to intensify these consequences, forcing cities to rethink common infrastructure practices. One popular method of adaptation is green infrastructure implementation, which has been found to reduce local temperatures and alleviate excess runoff when installed effectively. As cities continue to change and adapt, land use/landcover modeling becomes an important tool for city officials in planning future land usage. This study uses a combination of cellular automata, machine learning, and Markov chain analysis to predict high resolution land use/landcover changes in Philadelphia, PA, USA for the year 2036. The 2036 landcover model assumes full implementation of Philadelphia’s green infrastructure program and past temporal trends of urbanization. The methodology used to create the 2036 model was validated by creating an intermediate prediction of a 2015 landcover that was then compared to an existing 2015 landcover. The accuracy of the validation was determined using Kappa statistics and disagreement scores. The 2036 model successfully met Philadelphia’s green infrastructure goals. A variety of landscape metrics demonstrated an overall decrease in fragmentation throughout the landscape due to increases in urban landcover.


Introduction
On a global scale, urbanization is predicted to continually harm important ecosystem services far into the future [1], causing continuous challenges for governments, policymakers, and urban planners in resource reallocation [2].Many communities are combating the consequences of urbanization through policies focused on nature-based solutions (NBS).NBS policies encourage actions that help societies address a variety of environmental, social, and economic challenges in sustainable ways [3], using the science-policy-practice interface [4].Green infrastructure (GI) is a common example of NBS.The concept of GI describes the interdependence of land conservation and land development, and refers to a contiguous, interconnected green network consisting of a range of natural environments [5].Green infrastructure can enhance ecosystem function in fewer, larger areas compared to numerous, small patches [5], but connectivity cannot be achieved by solely enlarging the total area of GI.
Modeling potential urban scenarios and solutions has emerged as a useful tool to explore uncertain futures in complex urban systems and to further understand the impacts from land use/landcover changes (LULCCs).Although scenarios usually lack quantified probabilities [6,7], they instead function as alternative narratives that present important possibilities about the future [6][7][8].Using satellite images from the past and present via remote sensing techniques [9,10], researchers can calculate patterns of urbanization, apply drivers of change, and extrapolate LULCC trends into the future.To assess the best course of action, different land use policies can be applied to models by adding constraints and/or incentives.
To evaluate the consequences of urbanization and the validity of possible NBS, social and environmental scientists are increasingly using highly detailed LULCC models [11,12].Landcover models have been used to address general questions of landcover change and urbanization around the world [2,9,10,[13][14][15][16][17][18][19][20][21][22][23][24][25][26][27]; however, only one other study models LULCCs under GI policies [28].To predict precise landcover transitions and to answer specific questions of policy, future LULCCs need to be modeled at finer scales.Urban modeling studies are conducted at a variety of resolutions depending on the satellite imagery available.Most landcover models are created at a 30 m resolution using Landsat imagery [2,9,13,[15][16][17][18]27]; yet, small landcover features, like GI, require modeling at a much higher resolution, as GI projects can be smaller than 30 m.Similarly, urban models have been created at different levels of detail with varying numbers of landcover classes.Some studies present a broad overview of urbanization with only two landcover classes [10,19,[21][22][23][24][25][26][27], usually "urban" and "nature" or "nonurban."Other studies present more realistic models with seven to ten landcover classes representing many of the features in the urban system [13,15,16,18], such as buildings, roads, trees, and grass.A specific landcover class, such as GI, must be modeled with a larger number of landcover classes, as to accurately represent the landscape and the specific variables that effect GI's location.
A variety of different LULCC modeling tools exists today, all allowing for the prediction of socio-environmental changes in a study area over time and projecting these changes into the future in a way that it relates to measured land change [29].Modeling LULCCs involves historical estimates of landcover combined with biophysical and socioeconomic information to create estimates of future change [30].Selecting the method to model LULCC is an important first step, based on a study's purpose and available data [27].This study uses a hybrid modeling approach, utilizing two different tools to model future GI growth and continued urbanization.
One common method in land change prediction is cellular automata.Cellular automata (CA) models are based on the interaction of several components: the grid space can be a one, two, or multidimensional space; and the "cell" or the "automaton" is a discrete variable that represents the structural units of the grid [28].The "cell state" describes the characteristics of the cell which are subject to change [28].The change occurs according to specified transition rules, which are mathematical expressions that govern changes of the cell state [28].Cellular automata methods have been widely used in various modeling tools, such as SLEUTH [31], the CA-Markov module in TerrSet [32], and GEOMOD [33].GEOMOD was chosen in this study to model future growth in GI, as GEOMOD allows the user to define the quantity of change over time.
Machine learning is a newer approach that is gaining popularity in LULCC studies.Machine learning describes the automated procedures with which the knowledge can be acquired [34].This study utilizes the machine learning process of the multilayer perceptron (MLP) neural network in the Land Change Modeler (LCM) tool in TerrSet 18.3 [32].The MLP uses a back-propagation learning algorithm, one of the most widely used neural network models [35], to calculate transition potentials over time.The transition potential models are then combined in a Markov chain process to determine the overall quantities of change over time [35].LULCC models created in the LCM have been found to be more accurate than LULCC models created in other tools, such as SLEUTH and FUTURES [27].The LCM was used to model continued urbanization in this study.
Using detailed landcover models at a 1 m resolution, this study aims to model GI policies in Philadelphia, Pennsylvania, USA, and the resulting LULCCs from new GI and continued urbanization for the year 2036.A hybridization of GEOMOD and the LCM was used to model future changes in GI and urbanization, respectively.The model is driven by patterns of historical change, conservation and development constraints, the physical landscape, and distance variables.The model is validated using Kappa statistics and disagreement scores, and the landcovers are compared using a variety of spatial metrics.

Study Area
The city of Philadelphia, Pennsylvania, USA is located towards the eastern coast of the United States at 39.95 • N, 75.17 • W. The city experiences all four seasons with well-distributed precipitation throughout the year.The city is highly urbanized with an average population density of 4491 people per km 2 [36], and a number of universities, parks, and vacant lots for green space.Despite stated efforts to enhance and increase green space in the city [37], from 2008 to 2015, the city increased its urban areas by 11%, while its natural areas decreased by approximately 15% (Table 1).In 2009, the Philadelphia Water Department (PWD) initiated their latest plan to reduce combined sewer overflow events, Green City, Clean Waters (GCCW), to meet the regulations of the Pennsylvania Department of Environmental Protection and comply with the federal Clean Water Act [38].GCCW involves the use of GI as a way to alleviate the amount of runoff that flows into storm drains and, eventually, into the Schuylkill and Delaware Rivers and their tributaries.PWD measures the program's success using the concept of greened acres (GAs), or enough GI to manage one inch of stormwater from one acre of drainage area; approximately 27,158 gallons (103 cubic meters) [39].Greened acres are calculated with the following formula (Equation ( 1)): where IC is the impervious cover transformed into GI in acres.This quantity can include the area of the stormwater management feature itself, as well as the area that drains to it.Wd is the depth of water over the impervious surface that can be physically infiltrated into the ground in inches [39].This program is the first in the United States that prioritizes GI over traditional grey infrastructure to moderate stormwater runoff.Construction of GCCW projects officially began in 2011 and will continue to be implemented until 2036; however, some GI was built before 2011.GI projects are funded through credits to private developers, grants, and public works projects.In an effort to create more GAs, GCCW puts forth eight different best management practices to reduce the amount of impermeable surfaces within the city: 1) Green Streets, 2) Green Schools, 3) Green Public Facilities, 4) Green Parking, 5) Green Open Space, 6) Green Industry, Business, Commerce, and Institutions, 7) Green Alleys, Driveways, and Walkways, and 8) Green Homes [38].PWD implements a variety of GI practices throughout the city, including downspout planters, green roofs, rain barrels, tree trenches, bump-outs, stormwater planters, pervious pavement, wetlands, and rain gardens [38].By the year 2036, GCCW will have concluded with at least 9564 GAs, reducing stormwater pollution by 85% from its 2009 levels [38].

Data and Preprocessing
To create and validate future landcover distribution, at least three past and current landcover datasets are required.We acquired a 2008 Philadelphia landcover dataset [40], and generated landcover datasets for 2010 and 2015 (Table 2).City boundaries were set using city limits spatial data from the City of Philadelphia [41].In 2011, the University of Vermont Spatial Laboratory created a 1 m resolution 2008 landcover of Philadelphia for the city government [37].This landcover is accepted and highly used by the city government.The 2008 landcover was created using object-based image analysis techniques (OBIA) to extract landcover information from a 2008 orthophotography and 2008 LiDAR LAS data [40].Ancillary data sets were stacked on top of the OBIA data, which included shapefiles of building footprints, roads and railroads, and hydrography provided by the City of Philadelphia [40].For the purpose of this study, 2008 GI spatial data [42] was stacked onto the 2008 landcover (Figure 1).
To capture the impact of Philadelphia's GI polices, 1 m resolution 2010 and 2015 landcovers of Philadelphia were created by the Kremer lab at Villanova University.Both of the 2010 and 2015 landcovers were created using a supervised classification on 1 m aerial imagery from the United States' National Agriculture Imagery Program (NAIP) for 2010 and 2015 [43,44], respectively, in ESRI ArcGIS 10.5 [45].Similar to the 2008 landcover, ancillary data sources were stacked on top of each supervised classification, which included shapefiles of building footprints, roads, impervious surfaces, railroads, and hydrography (Table 2).After the landcovers were checked for accuracy, GI spatial data provided by City of Philadelphia [42] were stacked into each landcover (Figure 1).Overall, eight landcover categories were used: tree canopy, grass/shrubs, bare soil, water, buildings, roads/railroads, other paved surfaces, and green infrastructure.Although it was not ideal to use three different landcovers created with two different processes, the OBIA technique used for the 2008 landcover was not replicable within the scope of this project, in part due to over 30,700 manual edits it took to create the landcover [40]; however, the 2010 and 2015 landcovers created with the supervised classification achieved similar results and accuracy (Table 3).The 2010 and 2015 landcover datasets were validated using the accuracy assessment suite in ArcGIS 10.5 [45].The Accuracy Assessment suite uses ground truthing points and a confusion matrix analysis to calculate user's accuracy, producer's accuracy, overall accuracy, and Kappa coefficient scores.Kappa coefficient is a measure of the proportional improvement by the classifier over a purely random assignment to classes [20].The user's accuracy measures the proportion of each landcover class that is correct, whereas the producer's accuracy measures the proportion of the land base that is correctly classified [20].The classified images were assessed for accuracy based on a stratified sample of 105 reference points for each time period, which were visually evaluated using the NAIP imagery.The overall accuracies of the 2010 and 2015 classifications were, respectively, found to be 80% and 90.5%, with Kappa coefficients of 0.767 and 0.889 (Table 3).
Additional datasets used in this study are listed in Table 2. Yearly slope datasets were derived from their respective digital elevation models (DEM) using the slope tool in TerrSet 18.3 [32].Distance from roads and rivers were calculated using the Euclidean distance tool in ArcGIS 10.5.All of the data was projected into the UTM 18N projection.Although it was not ideal to use three different landcovers created with two different processes, the OBIA technique used for the 2008 landcover was not replicable within the scope of this project, in part due to over 30,700 manual edits it took to create the landcover [40]; however, the 2010 and 2015 landcovers created with the supervised classification achieved similar results and accuracy (Table 3).The 2010 and 2015 landcover datasets were validated using the accuracy assessment suite in ArcGIS 10.5 [45].The Accuracy Assessment suite uses ground truthing points and a confusion matrix analysis to calculate user's accuracy, producer's accuracy, overall accuracy, and Kappa coefficient scores.Kappa coefficient is a measure of the proportional improvement by the classifier over a purely random assignment to classes [20].The user's accuracy measures the proportion of each landcover class that is correct, whereas the producer's accuracy measures the proportion of the land base that is correctly classified [20].The classified images were assessed for accuracy based on a stratified sample of 105 reference points for each time period, which were visually evaluated using the NAIP imagery.The overall accuracies of the 2010 and 2015 classifications were, respectively, found to be 80% and 90.5%, with Kappa coefficients of 0.767 and 0.889 (Table 3).
Additional datasets used in this study are listed in Table 2. Yearly slope datasets were derived from their respective digital elevation models (DEM) using the slope tool in TerrSet 18.3 [32].Distance from roads and rivers were calculated using the Euclidean distance tool in ArcGIS 10.5.All of the data was projected into the UTM 18N projection.

LULCC Model and Validation
The 2036 LULCC model was created using the Land Change Modeler (LCM) in TerrSet 18.3 [32].The LCM was used to model continued urbanization based on spatial patterns from 2008 to 2015; specifically, landcover transitions to buildings, roads/railroads, and other paved surfaces.The LCM procedure involves change analysis, determining drivers of change, applying rules and restrictions, Markov chain-based transition predictions, and validation of the model (Figure 2).

Change Analysis
The landcover maps of Philadelphia in 2008 and 2015 (Figure 1) were analyzed for patterns of change, and exact quantities of transition between different landcover classes were calculated.Only transitions representing continued urbanization (i.e., change from trees, grass, and/or soil to buildings, roads/railroads, and/or other paved surfaces) were analyzed, as they were the transitions of interest.

Drivers of Change
Driver variables in this analysis included DEMs, slope, distance to existing roads and rivers, and evidence likelihood rasters.Distance to roads and distance to rivers were set as dynamic factors to be recalculated over time because as roadways and rivers change over time, so do the distances to these features.Elevation and slope have been documented to influence the location of urban growth [2,17,54].Evidence likelihood rasters were created for each LULCC that occurred.The evidence likelihood tool in the LCM transforms categorical variables, such as change from one landcover class to another, into numerical values so that they can be used in the modeling procedure [35].
Urban transitions were grouped into respective submodels by what they transitioned into and then used to compute transition potentials.A transition submodel consists of a group of transitions that share the same underlying driver variables [10], and so they can be modeled at once.The modeling of transition potentials is necessary for determining spatial change [55].The output of this step generates a series of transition potential maps that each correspond to a landcover transition based on the previous change analysis [55].The transition maps consider the suitability of pixels that have transformed into urban pixels based on a number of driving factors used for modeling processes of historical change.The transition potentials are created using a multilayer perceptron (MLP) neural network, which allows for transition submodels to be modeled at once [55].The MLP neural network is a feedforward neural network in which data flows in one direction from an input layer to an output layer through a number of hidden layers in between, as set by the user [56].A small number of hidden layers were used in this study, which expressed the common underlying themes in the variables [35] and resulted in higher accuracy levels.The computing elements (nodes) are grouped into layers, and each node receives an input signal from other nodes.After processing the signals locally through a transfer function, it outputs a transformed signal to other nodes for the final result [57].Each signal feeding into a node in a subsequent layer has the original input multiplied by a weight with a threshold added, and is then passed through an activation function [57] that, in this study, was non-linear.The weights are determined in the automatic training process before the network can be used for prediction purposes, aiming at changing the weights as to minimize the error between the observed and the predicted outcomes [57].Due to the non-linearity of the data, a sigmoid factor of 0.5 was applied to the weighted sum of inputs before the signal passed to the next layer.

Rules and Restrictions
Rules and restrictions were also applied to the model.A Boolean layer expressing areas of conservation and restricted development were added to the model so that urbanization would not occur in areas under conservation and the local Philadelphia airports.This spatial layer includes local urban parks, federal and state conservation areas, and the two airports in Philadelphia.

Transition Predictions Using Markov Chain Analysis
In the final change prediction, the LCM uses the change rates calculated in the change analysis step, as well as the transition potential maps to predict a future scenario for 2036.This step is responsible for determining the quantity of change to urban landcover in each transition using Markov chain analysis [55].The Markovian process is a method in which a predicted system can be estimated by finding its previous state and the probability of conversion from one state to another [18].By utilizing the 2008 and 2015 landcover maps, the Markov chain analysis determines, precisely, how much land would be anticipated to transition from 2015 to 2036, on the basis of projection of the transition potentials.A hard prediction was created in this study which yielded a projected map of 2036, where each pixel is assigned one landcover class-the class that it has been calculated to most likely become.

Validation
Validation was used to ascertain the quality of the predicted 2036 map.In order to validate the methods used to create the 2036 landcover map, the methods were replicated to create a predicted map of 2015 and then validated against the actual 2015 map.The 2008 and 2010 landcovers were used to predict changes in 2015.The validation module in TerrSet was used to statistically assess the quality of the 2015 predicted map against the 2015 reference map.The Kappa variation techniques were used in the validate module as the statistical validation procedure.Three variations of Kappa were calculated: Kappa for no information, Kappa for location, and Kappa standard.The Kappa for no information (K no ) signifies overall accuracy obtained in the simulation run, while Kappa for location (K location ) measures agreement level in a location [58].Considering the difficulty in interpretation encountered with the Kappa for correctly assigned proportion against the proportion of incorrectly assigned by change (K standard ) [59], the K standard was not used in this study; however, K location was useful for the validation in the absence of K standard [60].
However, Pontius and Millones (2011) have presented that Kappa scores can be useless, misleading, and/or flawed for practical applications in GIS and remote sensing [60].Instead, Pontius has developed other summary parameters calculated in TerrSet's validate module, two of which are presented in this study: disagreement at the grid cell level and disagreement due to quantity.Disagreement at the grid cell level is defined as the amount of disagreement associated with the fact that the comparison map fails to specify perfectly the correct locations of categories at grid cells within strata [35].Disagreement due to quantity is defined as the amount of disagreement associated with the fact that the comparison map fails to specify perfectly the correct quantity of each category according to the reference map [35].

Modeling Green Infrastructure
In order to assess the future implementation of GI, a cellular automata optimization approach was utilized.The GEOMOD module in TerrSet was used to assess the landcover change between GI and a "non-GI" category, which combined all other landcover categories (Figure 3).For the purpose of this analysis, the 2008 and 2015 landcovers were reclassified to GI and "non-GI," as GEOMOD only allows for the transition from "state one" to "state two." Green infrastructure was modeled separately from all other LULCC in GEOMOD because GEOMOD allows the user to specify the forecasted quantity of change.Philadelphia's GI plan indicates the amount of future GI by 2036, broken down by watershed [39], which was incorporated into the model.Future scenarios of GI are identified using driver variables, suitability maps, and calculating the amount of change in each Philadelphia watershed.Currently, there are no spatial plans for the future distribution of GI, as the city develops GI projects as they are needed with target quantities each year, and so it is challenging to evaluate the accuracy of this GI model, and it is just one possible scenario out of many possible scenarios.

Drivers of Change
GEOMOD uses an optimization algorithm that allocates cells to the most suitable locations based on the user-provided site suitability data [27].Similar to the LCM analysis, DEM, slope, distance, and evidence likelihood rasters were used in GEOMOD as the site suitability data.The suitability image is created by computing a weighted sum of all the driver images for each grid cell [33].

Change Allocation
Locations of land change are selected based on three decision parameters [33].The first parameter only allows for one directionality of change to occur.Specifically, in this study, all other land categories can change to GI, but GI pixels cannot change to anything else.The second parameter allows for regional stratification, where land change can be simulated differently for a series of user

Drivers of Change
GEOMOD uses an optimization algorithm that allocates cells to the most suitable locations based on the user-provided site suitability data [27].Similar to the LCM analysis, DEM, slope, distance, and evidence likelihood rasters were used in GEOMOD as the site suitability data.The suitability image is created by computing a weighted sum of all the driver images for each grid cell [33].

Change Allocation
Locations of land change are selected based on three decision parameters [33].The first parameter only allows for one directionality of change to occur.Specifically, in this study, all other land categories can change to GI, but GI pixels cannot change to anything else.The second parameter allows for regional stratification, where land change can be simulated differently for a series of user defined boundaries [27].Finally, neighborhood constraints that restrict growth within any time step to only edge cells between two classes can be set.Unconstrained growth was used in this study as different GI projects are not usually constructed right next to one another.

Amount of Change
The nature of the green infrastructure policy in Philadelphia requires a model restriction that specifies how much change will occur between state 1 (non-GI) and state 2 (GI) over a given amount of time.User-defined quantity of change was necessary for this analysis since it can be predicted how much GI will exist by 2036 according to GCCW.The prediction, in this case, focuses on the most likely location of the GI pixels.Information on 2016 green acreage from PWD's five year review, Green City, Clean Waters-Evaluation and Adaptation Plan [39], were used to calculate the amount of GI needed by 2036.PWD calculated the amount of green acreage they had accomplished in each of the four main watersheds in Philadelphia.The percent of GAs of the total 837 GAs was calculated for each watershed in 2016 (Equation ( 2)).This percentage was assumed to remain the same for 2036 predictions.Using GIS and the GI shapefiles provided by PWD [42], a ratio of physical acres of GI to GAs was then calculated for each watershed (Equation (3)).Predicted physical acres of GI in 2036 were then calculated using the ratio for each watershed (Equations ( 4) and ( 5)).
Physical acres GI/GA= GI ratio , per watershed in 2016, (3) 9564 GAs•GA% = GA 2036 , per watershed, (4) where GA 2016, WS is the GA in 2016 for each watershed, 837 is the total amount of greened acres in 2016 [39], GA% is the percentage of greened acres per watershed, GI ratio is the ratio of physical acres of GI to greened acres, and 9564 GAs is the target amount of greened acres in 2036 according to GI policies [39].Three smaller watersheds were not counted towards PWD's GAs goal.To account for the growth of GI in each of the smaller watersheds over time, the rate of change for each region was calculated and then extrapolated to 2036.The calculations for each watershed are displayed in Table 4.

Validation
The GI predictions were validated using the same approach as the LCM validation process.An intermediate GI prediction of 2015 was created in GEOMOD based on the 2008 landcover using the aforementioned methods.Again, the validate module in TerrSet was used to statistically assess the 2015 GI prediction against the 2015 GI reference image using Kappa statistics and disagreement scores.

Final 2036 LULC
The results of the GI-specific 2036 prediction were integrated into the 2036 predicted urbanization dataset.The end result is a predicted landcover dataset that includes the spatial distribution of predicted GI locations and continued intra-urbanization (Figure 4).* Data from Green City, Clean Waters-Evaluation and Adaptation Plan [45].† Data calculated using GIS and the PWD GI polygons [48].

Landscape Metrics
To understand the changes in the connectivity of green space and urban areas, each landcover was recategorized into two classes: green space (tree canopy, grass/shrubs, bare soil, and GI) and urban (buildings, roads/railroads, and other paved surfaces).To quantify the changes in the landscape over time, landscape metrics were calculated using Patch Analyst 5.0 [61], which integrates the spatial metrics found in FRAGSTATS [62] into ArcGIS.The number of patches, mean patch size, largest patch index (LPI), and patch cohesion index (PCI) were calculated at the class level.Patch cohesion index (PCI) and Shannon's diversity index (SHDI) were calculated for each landcover at the landscape level.FRAGSTATS 4.2 [62] was used to calculate the LPI and PCI for each landcover.Due to restrictions in FRAGSTATS 4.2, the landcover datasets had to be resampled to a 5 m resolution.A description of each metric can be found in Table 5.The physical connectedness of the corresponding patch type.PCI approaches 0 as the proportion of the landscape comprised of the focal class decreases and becomes increasingly subdivided and less physically connected, and vice versa.Dimensionless 0 < PCI < 100 Shannon's Diversity Index (SHDI) Relationship between the number of classes, the total number of patches, and the relative abundance of patches in each class.It has a value of 0 when no diversity is present and increases as the landscape becomes more fragmented.

Validation Results
The K location score (Table 6) was used for assessing accuracy, as it is difficult to interpret the Kappa for correctly assigned proportion against the proportion of incorrectly assigned by change (K standard ) [59].The LULCC prediction created with the LCM received a K location score of 78% (Table 6).The GI prediction created with GEOMOD received a K location score of almost 100% (Table 6).Most of the cells in the landcovers used in the GEOMOD analysis were the "non-GI" class, contributing to the high accuracy of the GI prediction.Additionally, in both the urbanization and the GI predictions, the disagreement scores represent very little disagreement between the models and their respective reference maps (Table 6).

Landcover Changes
The final 2036 landcover is featured in Figure 4. Overall, buildings are predicted to exhibit the largest net increase, followed by roads/railroads by 2036 (Figure 5).Grass/shrubs are predicted to exhibit the largest net decrease, followed by other paved surfaces by 2036 (Figure 5).The final 2036 landcover is featured in Figure 4. Overall, buildings are predicted to exhibit the largest net increase, followed by roads/railroads by 2036 (Figure 5).Grass/shrubs are predicted to exhibit the largest net decrease, followed by other paved surfaces by 2036 (Figure 5).The model was able to successfully meet Philadelphia's GI goal of at least 9564 GAs by 2036.The 2036 prediction allocated enough GI for 9945 GAs in 2036.A large portion of the new GI is predicted to be located in the middle and south Philadelphia (Figure 6).The new GI is predicted to mostly replace other paved surfaces, followed by grass/shrubs (Figure 7).The model was able to successfully meet Philadelphia's GI goal of at least 9564 GAs by 2036.The 2036 prediction allocated enough GI for 9945 GAs in 2036.A large portion of the new GI is predicted to be located in the middle and south Philadelphia (Figure 6).The new GI is predicted to mostly replace other paved surfaces, followed by grass/shrubs (Figure 7).
Land 2019, 8, 28 14 of 21 The final 2036 landcover is featured in Figure 4. Overall, buildings are predicted to exhibit the largest net increase, followed by roads/railroads by 2036 (Figure 5).Grass/shrubs are predicted to exhibit the largest net decrease, followed by other paved surfaces by 2036 (Figure 5).The model was able to successfully meet Philadelphia's GI goal of at least 9564 GAs by 2036.The 2036 prediction allocated enough GI for 9945 GAs in 2036.The new GI is predicted to mostly replace other paved surfaces, followed by grass/shrubs (Figure 7).A large portion of the new GI is predicted to be located in the middle and south Philadelphia (Figure 6).Buildings demonstrate a large net increase of almost 73 km 2 (Figure 5).A small reduction in buildings is due to new GI (Figure 7).Similar to increases in GI, most of the gains in building area are predicted to result from a loss in other paved surfaces, followed by grass/shrubs and tree canopy (Figure 8 right).It is forecasted that the new buildings will be well distributed throughout the city, with some concentration in south Philadelphia (Figure 8    Buildings demonstrate a large net increase of almost 73 km 2 (Figure 5).A small reduction in buildings is due to new GI (Figure 7).Similar to increases in GI, most of the gains in building area are predicted to result from a loss in other paved surfaces, followed by grass/shrubs and tree canopy (Figure 8 right).It is forecasted that the new buildings will be well distributed throughout the city, with some concentration in south Philadelphia (Figure 8

left).
Land 2019, 8, 28 15 of 21 Buildings demonstrate a large net increase of almost 73 km 2 (Figure 5).A small reduction in buildings is due to new GI (Figure 7).Similar to increases in GI, most of the gains in building area are predicted to result from a loss in other paved surfaces, followed by grass/shrubs and tree canopy (Figure 8 right).It is forecasted that the new buildings will be well distributed throughout the city, with some concentration in south Philadelphia (Figure 8

Spatial Metrics
SHDI is forecasted to reach its lowest levels of fragmentation throughout the entire landscape by 2036.The number of patches decreased in both categories from 2015.Mean patch size increased for both classes from 2015.However, the largest patch index decreased for green space over time but increased for urban areas.The patch cohesion index decreased over time for green space, increased slightly over time for urban areas, and increased to its highest levels for the entire Philadelphia landscape (Table 7).

Data Resolution
This research aimed to model fine-scale, multiclass prediction of intra-urban landcover change under GI policies.This study is unique in that it utilized one-meter resolution aerial imagery from the United States' NAIP to create the base landcover datasets.Most LULCC modeling uses Landsat data at a 30 m resolution or larger; however, GI can easily be less than 30 m in size and so higher resolution data is needed to map and model GI.
Remote sensing techniques and the availability of free to low cost satellite imagery and their temporal frequency has greatly enhanced the monitoring of urban growth and land use dynamics around the world [9].As technology continues to improve, the resolution of satellite images will also be enhanced.Developing and testing LULCC models with fine-scale data, as we do here, allows for more detailed models of the dynamics of change within urban environments.Advanced models can aid policymakers and planners in analyzing the effects of smaller landscape features, such as GI.

Green Infrastructure
This study spatially modeled the increase of GI in Philadelphia under the city's GI program for the year 2036.Many studies model future urban growth in areas around the world; however, we found only one other study that also modeled LULCC under GI policies [28].Contrasting to Mitsova et al. (2011), Philadelphia's 2036 model predicts an overall increase in cohesion throughout the entire landscape under GI policies (Table 6), which can be attributed to the decreased number of patches from 2015.Green space patches decreased by 920,970 patches from 2015 to 2036, with the mean patch size also increasing slightly (Table 6).
Although the number of patches decreased and mean patch size increased, the patch cohesion metric for green space still decreased over time (Table 6).GI serves as a nature-based solution to stormwater management in Philadelphia; nevertheless, to maximize the potential of GI's ecosystem services, additional attention is needed regarding the type of GI implemented and to increase its overall connectivity with other green spaces.GI and green spaces enhance ecosystem function in fewer, larger areas compared to numerous, small patches [5], but connectivity cannot be achieved by solely enlarging the total area of GI.
The model predicts that the new GI to be added by 2036 will be fairly well dispersed throughout the city, with the middle of the city and the southern portion of Philadelphia gaining some of the largest increases (Figure 6 top).A majority of the new GI will replace other paved surfaces (Figure 6), which includes urban features such as parking lots, alleys, and sidewalks.Replacing a majority of other paved surfaces with GI is a plausible prediction, as many of the GI strategies that the city currently uses does replace paved surfaces.These GI strategies include trenches, bump-outs, planters, and pervious pavement [63].A small portion of the expected GI is predicted to replace tree canopy, which is unlikely as the cities would probably not cut down trees for GI development.This prediction could be possible if tree trenches are added to these areas, so the tree could be preserved but, also, the trench underneath it would collect water for later use.
This study only used geophysical driver variables, whereas, in reality, other socioeconomic drivers have an additional influence on the spatial distribution of GI projects.To improve the GI model, GI experts at the City of Philadelphia should be consulted to identify other variables and constraints that may influence GI location.

Urbanization Land Use/Landcover Changes
Unlike the results presented by Mitsova et al. (2011), where urban growth is forecasted to stall, urban areas in Philadelphia are still predicted to increase, even under GI policy.Continued urbanization in the form of new buildings can be seen in south Philadelphia around the Delaware and Schuylkill rivers (Figure 8).Philadelphia's population is newly re-growing, as it has increased steadily over the past seven years, and is predicted to continue to grow [64], driven by urban regeneration and new employment opportunities in the city.By 2036, urban spaces are predicted to decrease by 189,436 patches across Philadelphia, but the mean patch size and largest patch index still increase from the 2015 measurements (Table 6), indicating further urban development in areas that are already urbanized.Specifically, buildings will experience the largest growth of the urban classes (Figure 5).With the addition of new urban areas throughout the city, the patch cohesion index for urban spaces continues to increase by 2036.
Road/railroad area also increases in 2036, approximately doubling from the 2008 area (Figure 5).New roads and railroads will be needed to meet the transportation needs of the new buildings, residents, and businesses.Other paved surfaces, such as sidewalks and parking lots, decrease in area over time as they are lost to new GI and buildings (Figure 5).
The amount of urban growth forecasted by the model should be analyzed further.The Delaware Valley Regional Planning Commission predicts that the population of Philadelphia will increase by only 6.37% from 2015 to 2035 [65].Under this population growth, the model most likely overestimates urban growth, as buildings are projected to have a net increase of almost 25%, and roads/railroads to have a net increase of approximately 23%.The demand for new buildings with only a 6.73% increase in population will not meet the rate projected by the model.Many socioeconomic factors affect urban and population growth rates, and they should be added to future versions of this model to accurately assess urban land change.

Future Research
There is little research on landcover changes in Philadelphia [66], and only one study that applies GI policies to LULCC using similar methods [28].This study aims to fill those gaps, as Philadelphia is a complex urban system and one of the few cities in the United States that prioritizes GI over traditional grey infrastructure.As the city continues to develop green infrastructure, it may act as a model for other cities in need of implementing nature-based solutions.
The methodology and the resulting landcover models used in this study can serve as a base for further exploration into how GI will function in the future.As a result of urban expansion, cities around the world are finding alterations of energy budgets with modifications to climatic, hydrologic, and biogeochemical cycles, and habitat fragmentation leading to a reduction in biodiversity [67].In Philadelphia, GI can serve as a solution to these problems.The degree to which GI can resolve these issues can be assessed by using future landcover models with predicted GI distribution as the basis of ecosystem service analyses.Additionally, further analysis of GI distribution based on socioeconomic equity should be studied in the future, but it was beyond the scope of this paper.This is integral to assessing future risk and vulnerability to environmental phenomena, such as climate change [17].

Conclusions
This study finds that LULCCs from GI policies can successfully be modeled in a heterogeneous intra-urban environment using fine-scale data.As Philadelphia continues to grow, GI will be implemented throughout the city to meet its GI program goals.Philadelphia city planners should consider GI cohesion to expand the many ecosystem services that GI contributes to the city.
Additionally, this study highlights the usefulness of LULCC modeling as a method for planning for the future.Availability of data and imagery will only improve in resolution and temporality in the future, allowing for more accurate landcover models that can capture smaller features, such as GI.As policymakers and city planners begin to plan for increased populations and the aggravated consequences of climate change, nature-based solutions will receive greater implementation as they allow for inexpensive and creative methods for city adaptation.
The 2036 LULCC model was created using the Land Change Modeler (LCM) inTerrSet 18.3 [32].The LCM was used to model continued urbanization based on spatial patterns from 2008 to 2015; specifically, landcover transitions to buildings, roads/railroads, and other paved surfaces.The LCM procedure involves change analysis, determining drivers of change, applying rules and restrictions, Markov chain-based transition predictions, and validation of the model (Figure2).

Figure 2 .
Figure 2. Processes involved in the Land Change Modeler (LCM) to create and validate the 2036 LULCC of Philadelphia.
Supervised Classification in ArcGIS +City of Philadelphia Spatial Data

Figure 2 .
Figure 2. Processes involved in the Land Change Modeler (LCM) to create and validate the 2036 LULCC of Philadelphia.

Figure 3 .
Figure 3. Process to create and validate the 2036 GI prediction using GEOMOD model based on Philadelphia's GI plans outlined in Green City, Clean Waters [38].

Figure 3 .
Figure 3. Process to create and validate the 2036 GI prediction using GEOMOD model based on Philadelphia's GI plans outlined in Green City, Clean Waters [38].

Figure 4 .
Figure 4.The predicted landcover of Philadelphia in 2036 under green infrastructure policies.

Figure 4 .
Figure 4.The predicted landcover of Philadelphia in 2036 under green infrastructure policies.

Figure 5 .
Figure 5. Net change for each land use category from 2008 to 2036.

Figure 6 .
Figure 6.Map of the increase in GI from 2008 to 2036 broken down by category.The middle and southern portions of Philadelphia are highlighted ("Other to GI" includes buildings, soil, and roads/railroads transitioning to GI).

Figure 5 .
Figure 5. Net change for each land use category from 2008 to 2036.

Figure 5 .
Figure 5. Net change for each land use category from 2008 to 2036.

Figure 6 .
Figure 6.Map of the increase in GI from 2008 to 2036 broken down by category.The middle and southern portions of Philadelphia are highlighted ("Other to GI" includes buildings and roads/railroads transitioning to GI).

Figure 6 .
Figure 6.Map of the increase in GI from 2008 to 2036 broken down by category.The middle and southern portions of Philadelphia are highlighted ("Other to GI" includes buildings, soil, and roads/railroads transitioning to GI).

Figure 7 .
Figure 7. Contributions to the increase in GI by area (km 2 ) from 2008 to 2036. left).

Figure 8 .
Figure 8. Changes in buildings from 2008 to 2036.Left: Map of change in buildings.Right: Contributions to increases in buildings.

Figure 7 .
Figure 7. Contributions to the increase in GI by area (km 2 ) from 2008 to 2036.

Figure 7 .
Figure 7. Contributions to the increase in GI by area (km 2 ) from 2008 to 2036. left).

Figure 8 .
Figure 8. Changes in buildings from 2008 to 2036.Left: Map of change in buildings.Right: Contributions to increases in buildings.

Figure 8 .
Figure 8. Changes in buildings from 2008 to 2036.Left: Map of change in buildings.Right: Contributions to increases in buildings.

Table 1 .
Land use/landcover changes (LULCCs) in Philadelphia from 2008 to 2015 in km2.Nature includes tree canopy, grass, and bare soil.Urban includes roads, other paved surfaces, and buildings.

Table 2 .
Datasets used in the study to develop the future projection of Philadelphia in 2036.

Table 3 .
Accuracy scores of the 2010 and 2015 supervised classifications of Philadelphia.

Accuracy * Tree Canopy Grass/ Shrubs Bare Soil Water Buildings Roads/ Railroads Paved Surfaces Overall Accuracy Kappa
* U Acc. = User's Accuracy; P Acc.= Producer's Accuracy.

Table 5 .
Description of spatial metrics.

Table 6 .
The Kappa scores and disagreement scores used to validate the predicted 2015 landcover against the reference 2015 landcover, and the 2015 GI against the reference 2015 GI.

Table 7 .
Spatial metrics of the Philadelphia landscape, urban areas, and green space: Shannon's diversity index (SHDI), number of patches (NP), mean patch size (MPS), largest patch index (LPI), and patch cohesion index (PCI).