Delimiting Future Urban Sprawl Boundaries Using a GIS-based Model for Ecological Sensitivity Index Assessment and Optimization Techniques. The case of Mytilene (Lesvos Island, Greece)

Delimiting urban sprawl boundaries have been generally regarded as a regulatory policy measure for controlling chaotic and sparse urban expansion and for protecting ecological areas. The conservation of ecologically sensitive areas plays a key role in environmental protection; so, harmonizing urban sprawling with nature conservation can be viewed as a binary compatibility planning problem. This study, aims to employ a geographical allocation model, based on the minimization of the environmental cost in order to apply complex spatial clustering techniques. Firstly, five ecological sensitivity factors affecting the ecological footprint of the study area are modeled through Geographic Information Systems (GIS) and Analytic Hierarchy Process (AHP) method in order to evaluate the Ecological Sensitivity Index. Then, several spatial objectives and constraints such as continuity and compactness are applied and finally, the most optimal areas are extracted for future urban sprawl. Spatial regulations, sitting rules considerations and scenarios based on the parameters of the spatial clusters outputs are tested to the commune of Mytilene, located on Lesvos Island, Greece, where strong land use changes have been recorded by the urban sprawl over the last three decades.


INTRODUCTION
Nowadays, it is imperative to manage environmental protection in order to protect natural resources. We live in an era where the terms of sustainability and risk expressing the relationship between human -nature is an integral part of our daily lives and our duty is to ensure a balance between these two systems.
A crucial question is how can we measure the balance between human activities and environment in a sustainable manner? The evaluation of ecological importance in a regional scale is to emphasize on the harmonious development between production space, living space, and ecological space. Ecological risk assessment is conducted in order to transform scientific data into meaningful information about the risk of human activities to the environment. The purpose is to enable planners, risk managers, and stakeholders to make informed environmental decisions. A set of new technologies contributes in this direction, such as Geographic Information Systems (GIS) and Spatial Optimization models, which in combination with a broader set of disciplines can be a useful management tool to address the problem.
Goal programming approaches present many advantages compared to specific algorithms: simplicity of implementation if a programming software is available, reliability of the method, computational speed, exact or guaranteed approximate solution of the problem, and finally, possibility of easily modifying the model. Kao 2012;Beyer et al 2016), 2) heuristic methods based on a) genetic algorithms (Haque & Asami, 2011;Cao et al 2012;Datta et al 2012;Stewart and Janssen 2014;Mohammadi et al 2015;), b) simulated annealing (Stewart et al., 2004;Aerts et al 2005;Sante΄-Riveira et al 2007) or c) particle swarm optimization Masoomi et al 2012;Sahebgharani 2016). It is important to highlight that the adoption of integer linear programming methods in land use management has been slow as long as early trials, tests and models proved unsatisfactory. The main reason was the lack of computing power when the size of the problem is spatially expanding (Williams and ReVelle 1998;Shirable 2005). Therefore, multiprocessing methods were applied to enforce IP model's capability through parallel computing environment (Liu and Kao 2012) and generative land-use modeling techniques enhanced with additional algorithms in order to solve the objectives faster (Ligmann-Zielinska et al. 2008). Contrary, heuristic algorithms are based on optimizing the objectives simultaneously in multi-objective mode focusing on Pareto front. Although, these methods do not generate an exact optimal solution, the generated solutions are meaningful and significant in many complex problems and case studies (Aerts et al. 2005).
Another important characteristic that characterizes MOLU is the structure of the input data and the ability the models have to interact with GIS. Many efforts have applied in recent years to combine spatial optimization procedures with GIS interfaces (Church 2002;Sante Riveira et al. (2008);Chen et al. 2010;Cao et al. 2012, Sanussi et al. 2014Stewart and Janssen 2014;Mohammadi et al. 2015;Ligmann-Zielinska 2017). GIS make use of two types of data: grid data (raster format) and attribute data (vector format). Data in a raster model are stored in a twodimensional matrix of uniform cells on a regular grid. By their nature raster data are substantially easier to include in mathematical representations of the world for purposes of optimization. Using a grid-based representation of a planning region, Stewart et al. (2004) and Janssen et al. (2008) showed that it was possible to formulate a spatial planning problem in mathematical terms and apply MOLU to generate optimal solutions interactively.

Compactness, contiguity and spatial constraints
Compactness and as a result contiguity is an issue that belongs to optimization spatial analysis which in turn applies diverse analytic and computational techniques in order to find optimal or near to optimal solutions (Vanegas P. et al. 2010). Land-allocation models in terms of the method used to encourage compactness and contiguity are separated within two types of constructs: solution-based and explicit constraint-based. The explicit constraint approach may be further divided into adjacency-based clustering, perimeter-based compactness, and block aggregation (Ligmann-Zielinska et al. 2008).
Contiguity can be explicitly structured in spatial optimization models or implicitly accounted for in a solution algorithm. Most explicit approaches are based on graph theory imposing network connectivity (Williams 2002;Shirabe 2005;Datta et al. 2012; (Billionnet 2013). Williams (2002) defined necessary and sufficient conditions for spatial connectivity. Rather than utilizing paths and spanning trees, Shirabe (2005) formulated contiguity constraints based on network flows. This work was extended by Duque et al. (2011) to include multiple network flows, one for each region. Datta et al. (2012) formalized their model as a multiobjective partitioning network problem using an Integer coded Genetic Algorithm.

Urban sprawl management
The increasing demand for land resources due to growth in population, urban areas, and economy has posed great challenges to rural and urban sustainable development in regional or national scale worldwide. (Haque and Asami 2014).
In recent years, great efforts have been succeeded on sustainable urban development and sprawl handling on the horizon of more compact and ecological friendly cities (Handayando et al. 2017). In particular, Kumar et al. (2016) with a Mixed Integer Quadratic Program (MIQP) consider just a single objective, maximizing suitability value of land to limit urban sprawl, based on two spatial constraints. Gabriel et al. (2006) on the contrary, take a multi-objective approach to controlling sprawl in land development by considering objectives from the perspective of the government, planners, environmentalists, conservationists, and land developers. Masoomi et al. The most comprehensive framework, according to the literature review, to quantify and measure sprawl is implemented more than a decade ago by Ewing et al. (2002). Based on causeeffect relationships of urban sprawl and quality of life, it includes numerous measures that are broadly divided into four categories, which are residential density, neighborhood mixture of homes, jobs, and services, strength of centers, such as business districts, and accessibility to the street network. Finally, Ligmann-Zielinska et al. (2008) developed probably the most robust optimization model for efficient utilization of urban space.

CONTRIBUTION OF THIS RESEARCH
Highlighting the most important factors in optimization procedures for land use management and ecological conservation the main purposes and contribution of this study are (1) to create a novel and integrated spatial index for evaluating and expressing the ecological sensitivity of a study area, (2) to demonstrate how the combination of theoretical models and GIS tools can contribute in environmental protection and land use management, (3) to develop a multiobjective land use allocation model for future urban sprawl boundaries expressed in a simple manner, (4) to achieve relevant computational time to obtain high quality and optimal compact sites through powerful programming techniques and finally (5) to incorporate these models into a GIS platform in order to be able to achieve multiple future scenarios and filtering out solutions as an interactive operation during a planning process.

Study Area and Data
The study area is the greater area of Mytilene, which is located in Lesvos Island in North-East Aegean Sea in Mediterranean Sea (Fig.1). The municipality covers an area of 107.5 km 2 with population of 37,890 inhabitants, according to the census of 2011. During the last three decades the land use changes regarding chaotic urban development lead to numerous impacts in the ecological and environmental status of the area. Data of land cover and ecosystem types were derived from the National Web Portal of Geospatial Data (http://www.geodata.gov.gr/). An extra spatial information added in Corine LC 2000 through digitalization as long as a lot of land use changes observed from 2000 to 2011. In order to create a detailed Digital Elevation Model (DEM) all the topographic elements (contours with contour interval 4 meters, heights, roads and drainage network) were digitized from maps of 1:5000 scale of the National and Geographical Military service.

Ecological sensitivity model
According to the related study of Troumbis (1995), the development of a mapping methodology assessing ecological sensitivity in natural systems, is an important condition to clarify the physical meaning of the term "ecological sensitivity". The term of "sensitivity" is not a specific variable and does not describe any physical property of a living system. Also, it can't be characterized by either as a single variable or as a complex descriptive parameter or factor which brings an etymological or conceptual dimension of the term "sensitivity". On the other hand, it is more accurate if we say that it is described as a pressure factor on natural systems, derived from human activities and is defined as a quality attribute that refers to the changes in the productivity and diversity of a natural system (mainly due to the land use changes) (Fig. 2.)

Combining ecological sensitivity model with GIS
The distinct steps of the framework are depicted in Fig 2 and according to these steps the model of relative quantification of ecological sensitivity arises from the combination of the following parameters rated as presented in Table 1. and consist of: 1) Vegetation, 2) Drainage basin system, 3) Heterogeneity of ecological landscape, 4) Anthropogenic disturbance sources (e.g. settlements, transport networks) and 5) Topographic slope. Each land-use type is the basic of the model regarding the expression of the productivity and the diversity of the system. The graded hierarchy of the different types of vegetation in relation to the productivity and diversity is achieved by using specific weighted values which reflect the system's relative productivity and diversity amounts.
The other criteria are used to change -"affect" the initial value of vegetation's specific gravity through the addition or subtraction of weight's values. Therefore, the highest values correspond to the most diverse vegetation types which have a reliable water supply and are distant from sources of human disturbances. The lowest sensitivity values correspond to intensive agricultural areas and settlements (Troumbis 1995). These values from higher to lower mean high ecological sensitive areas (negative values), ecological sensitive (close to zero) and none ecological sensitive (positive values).

Evaluation of ecological sensitivity indicators weighting
As ecological sensitivity is a multi-attribute index, it often requires a method combining qualitative and quantitative analysis to evaluate. The analytical hierarchy process (AHP) is a decision-making method and a theory of measurement through pairwise comparisons and relies on the judgements of experts to derive priority scales. In the AHP method, obtaining the weights or priority vector of the alternatives or the criteria is required. The decision-making

Vegetation Maps
Proximity to Streams

Productivity and Diversity
Landscape Heterogeneity Proximity to disturbance sources Accessibility

Hydrographic Network
Land Use Slope Figure 2. Geodatabase structure and procedural steps to determine the ecological sensitivity through GIS process starts with defining and dividing the problem into issues-criteria, which may optionally be divided further to form a hierarchy of issues-elements (Saaty 2008).

Productivity -Diversity
Natural fracture type Accessibility Slope

Proximity to Road Network
Distance 0-25m 2

Regional roads network
Distance 0-200m -4 Distance 200-500m -3 Gravel roads Distance 0-100m -2 Forest roads Distance 0-50m -1 Optimization model Formulating an optimization problem generally contains three basic steps: defining decision variable(s), formulating objective function(s) and defining problem constraint(s). The model is non-linear if the objective function and/or some of the constraints are non-linear. Moreover, Gabriel et al. (2006) highlight in their model the class of quadratic problems as part of nonlinear problems, however, the relaxed version of these problems are simply convex, quadratic programs with linear constraints and thus represent a reasonable computational burden given the state of the art in optimization solvers. In order to solve multiobjective optimization problem two main processes can be applied in general, (1) the weighting method and (2) the constrained method (Gabriel et al. 2006). Our algorithm is based on an Integer Quadratically Constrained optimization problem and Brunch and Bound algorithm upon linear constraints and weighted method approach, implemented through Gurobi optimization package (Gurobi Optimization 2017).

Size, data and complexity of the study area
In our model, as we referred to the literature, raster based formatted files are used (e.g. TIFF images) which are transformed into arrays as input for the optimization algorithm. Each raster cell represents ecological sensitivity index value. The dataset used as input for the model, in raster format, contains 287.028 cells (603 rows * 476 columns) with a cell size of 30 meters.
Compactness calculation for each cell is affected by neighbored cell's state related to ecological sensitivity index. For our model, we express as neighborhood of a cell the twodimensional square lattice composed of a central cell and its four adjacent cells (Von Neumman Neighborhood) (Fotakis and Sidiropoulos 2012). Among the total number of cells, NoData values of the raster file are not included but the model is capable to handle them in spite of the fact that these values do not participate in the final optimal solution. These can include borders, backgrounds, the sea or other data considered to not have valid values. Therefore, we have to customize the model to be able to search for an optimal solution above or near the boundary cells of the study area considering that the initial input raster is expressed as C i,j for i as the total number of rows and j as the total number of columns. In order to be able to account cross defined neighbors for the cells above the coastline we must add two additional rows and columns forming a final raster dataset X i,j for i = i + 2 and j = j + 2. Finally, the IQCP model contains 115.345 discrete variables (binary) and 232.414 constraints.

Decision variables, Objectives and Constraints
Minimum ecological sensitivity index: Ecological sensitivity index, and therefore, suitability of a land cell for urban growth is expressed as the status of location condition as it has formulated by MCA from previous steps. We apply multiple factor model of weighted summation to do the superposition calculation, ecological sensitivity index can be expressed by the following formula: where ES is the ecological sensitivity index, w i is the weight of each criteria i and f i is the rating of criteria i. Maximum compactness of urban growth patterns: A customized approach is applied for calculating compactness and continuity of cells based on Kao and Lin's (1996) model. In this study, compactness index is expressed as the total perimeter of the non-connected edges of the total cells leading the model to make urban patches regular and more rectangular. The continuity of the selected cells of the solution to the above model is guaranteed because the model seeks the smallest perimeter. Multiple factors optimization: For a problem with consideration of multiple factors, the objective function of the model should be modified adding the weight factor of each objective and can be expressed as: where: where I and J are the total number of rows and columns of cells; w u is the weight for compactness and w s is the weight for ecological sensitivity. S i,j represents the value of ES (Eq. 1) for each cell , V i,j is used to record the length of each cells perimeter; x i,j , an [0,1] indicator variable, is defined to represent whether cell i,j belongs to a considered site or not. Optimization constraints: The two fundamental constraints for our model refer indirectly to (Eq. 4) compactness calculation and neighbors accounting, NoData cell values handling and (Eq. 5) number of total amount of selected cells by the user, in order to calculate the objectives. N is the required size (in numbers of cells) of the desired site;

RESULTS
Sub -criteria The visualization output of each sub-criterion is shown in Fig.4 where the final maps of land use, road and drainage network proximity, slope and heterogeneity index are demonstrated.
Productivity and Diversity index: The most downgraded land units are observed in the central part of the study area as also to the South -East coastal part where urban sprawl areas and other individual industries are located. In contrary, the most productive areas are observed in North -West and South -East part of the study area where different types of forests, wetlands and sclerophyll vegetation exist. Proximity to disturbance sources index: The higher disturbance values exist because of the small highways that connect the city of Mytilene with the entire island as also the regional road network of the municipality that connects the different smaller settlements and covers the biggest area. Proximity to streams index: The drainage network was divided in different classes according to Stahler's method of designating stream orders (Strahler 1957). The biggest the stream order and the proximity is, the higher values of the index are observed. Landscape Heterogeneity index: Landscape diversity calculated to study the local landscape heterogeneity by drawing out a grid of 500m cell size and dimensions of 38*29 cells. Heterogeneity index reflects the diversity of the amounts of landscape elements, apart from areas with continuous urban development and the variation of their proportion. The higher the combination of different vegetation types is observed the less vulnerable is an ecosystem as long as the ecological footprint increases. Accessibility index: Higher slope values reduce ecosystem's regenerative potential as long as the capability of humans to intervene (e.g. forest fires) is decreased. The most vulnerable areas are observed in the South -East part of the municipality as the mountain area of Amali is located.

Weight determination and Ecological Sensitivity index
According to the calculation results derived from AHP, if CR < 0.1, where CR is the consistency ratio, the matrix is deemed to have the satisfactory consistency otherwise it should be adjusted. In our case, CR = 0,059 so the consistency value is accepted for number of criteria n = 5, λmax is the biggest eigenvalue, random inconsistency RI = 1,12 and consistency index CI = 0,067 (Malczewski 2008). The last step is to transform the ecological sensitivity single factor layers into the final raster map. The final results are presented in Fig.5 (no.1) where the ecological sensitivity index is calculated for the study area. Scale of the cell values ranges between -9 and 16, as a result, most downgraded areas are expressed with red and orange colors. Contrary, the most ecologically upgraded and stable areas are visualized with different green patterns. Proximity to disturbance sources, (4) Proximity to streams, (5) Landscape heterogeneity, (6) Accessibility Urban Sprawl Boundaries results Multiple scenarios of future urban growth patterns are presented in Fig.5 for number of total selected cells N = 100, 500 and 1000 in order to validate model's capability to allocate continuous and compact patterns related to the given weighting factors for ecological sensitivity and compactness. For small patterns (N=100) model returns absolutely squared areas irrespective of each criteria's weight. Searching for larger areas (N=500 and 1000) it is highlighted that as long as compactness weighting factor is maximized (>=0.6) contiguity and compactness are guaranteed. Although, is noted that for large clusters as compactness weight reduces (<=0.6), final optimal areas might be fragmented (e.g. in two separate clusters), something that seems normal considering that contiguity and compactness are expressed through the objective function and not as constraint factors.

CONCLUSION AND DISCUSSION
This study developed a goal programming-based framework with the integration of a multicriteria ecological assessment and external planning interventions, which aims to generate future urban sprawl land sites that minimize ecological sensitivity and compactness. The main modifications include: (1) include a different measure of ecological sensitivity; (2) implement a novel approach of formatting and calculating compactness (3) implement all these through a fast and easily handled model based on Integer Non-Linear programming.
Therefore, an integer quadratic constrained programming (IQCP) model for urban sprawl boundaries allocation was formulated. One of the main advantages of our approach is that it produces solutions that reflect precisely the default preferences of the user. Furthermore, it relies on a rigorous and systematic mathematical approach that avoids falling in sub-optimal solutions leading to optimal contiguous and compact sites in a promising amount of time considering the size of the spatial units. Future research consists of multiple clustering options (return more than one cluster), additional sitting objectives and constraints and a fully compatible toolbox with GIS platforms.