Modelling the future vulnerability of urban green space for priority-based management and green prosperity strategy planning in Kolkata, India: a PSR-based analysis using AHP-FCE and ANN-Markov model

Abstract Changes in land-use and land-cover (LULC) in urban areas affect the natural environment, especially urban green spaces (UGS). The present study examines the loss of UGS due to LULC transformation at different periods to predict the future vulnerable zone of UGS, based on the 'Pressure-State-Response’ framework. To calculate the weight of each factor, a combined Analytical Hierarchical Process and Fuzzy Comprehensive Evaluation method have been used. An integrated multilayer perceptron based artificial neural network and Markov chain (MLP-ANN-MC) model has been employed to predict the UGS vulnerable area in Kolkata. Results indicated that growth rates of built-up area, land-use dynamic degree, change intensity index, and proximity factors are the major responsible for UGS vulnerability. Applying the MLP-ANN-MC model, future vulnerable zones were identified for management and conservation of UGS. The methodology developed and demonstrated in this study expands LULC change analysis and provide a new dimension for UGS vulnerability assessment.


Introduction
Urbanization is one of the most dramatic forms of land alteration, adversely influencing urban environments and ecosystems (Lin et al. 2020). Land use land cover (LULC) change in urban areas often results from rapid urban expansion and has emerged as a serious environmental threat (Banti et al. 2019). The analysis of LULC change have recently been a focus of the research community to understand various environmental changes, including transformation in vegetation cover, shrinkage of wetlands and water bodies, loss of open space, and alteration of urban micro-climates (Ali and Patnaik 2019;Kim and Kwon 2021). Urban green cover is often adversely affected by the process of rapid and unplanned urban expansion (Xu et al. 2018;Ghosh et al. 2021). A considerable amount of research has been conducted by the academic community to assess the impacts of LULC change and alteration of urban green space (UGS) (Zhou et al. 2011;Zhao et al. 2013; Bardhan et al. 2016;Dinda et al. 2021). UGS refers to the green areas of a city, including forests, trees, parks, residential gardens, playgrounds, grasslands, and any other green cover (Taylor and Hochuli 2017;Nath et al. 2018;Wang et al. 2018). UGS provides multiple ecosystem services that support biodiversity, minimize air pollution, absorb CO 2 , reduce urban heat islands, improve residential quality of life, and provides various health benefits (Chen 2015;Anguluri and Narayanan 2017). Modelling UGS dynamics through LULC change analysis provides a perceived knowledge regarding the future vulnerability of green space (Shishir and Tsuyuzaki 2018). Moreover, UGS management is essential to mitigate the undesirable effects of urbanization and to support sustainable urbanization (Gavrilidis et al. 2019;Hou et al. 2020).
UGS vulnerability denotes the external factors which endangering the green space's health, obstructing the ecological system, and hampering the ecological security of green space (Daniels et al. 2018). UGS is one of the sensitive elements of the natural environment; hence it affected by various factors, driven by rapid urbanization and its associated alterations (Weng and Lo 2001;Ibrahim et al. 2014;Dinda et al. 2021). Therefore, to measure the vulnerability of UGS, emphasis should be given to assess the role of those factors. Scholars have frequently used the 'Pressure-State-Response' (PSR) model to measure ecological vulnerability (Malekmohammadi and Jahanishakib 2017;Sun et al. 2019;Wu et al. 2020;Yang et al. 2021). Moreover, the expert-opinion-based PSR model is widely used by researchers because it is a systematic and very effective method for impact assessment (Sun et al. 2019). The model was developed by the Organization for Economic Co-operation and Development to analyze the various environmental issues in the development context (OECD 1993). The PSR approach conceptualizes causal relationships among the ecological vulnerability indices, based on three categories, viz. 'pressure', 'state', and 'response'. Assessment of UGS vulnerability is essentially a multi-criteria decisionmaking (MCDM) process. Analytic Hierarchy Process (AHP) is a popular multi-objective analysis and MCDM method (Taheri et al. 2015;Dinda et al. 2019;Ghosh et al. 2019;Lin and Kou 2021). Moreover, AHP can elaborate and decomposes the complex relationship among criteria and solve decision-making problems . However, the sensitivity analysis becoming an essential element of the MCDM model (Banerjee et al. 2020). The Fuzzy Comprehensive Evaluation (FCE) is now a widely used method to dealing with interdependencies problems among the indicators and increase the robustness of the decision-making method (Chiu et al. 2014;Shao et al. 2020;Khan et al. 2021).
Likewise, the land-use change results can help land-use planning by identifying individual components of LULC, their patterns, and the direction of change. But, simulationbased LULC models can provide insights into the areas of future expansion associated with urban growth (Musa et al. 2017;Dahal et al. 2018;Abass et al. 2019). Many LULC models are available, including statistical methods, mathematical simulations, expert system methods, cellular automata models, and machine learning-based methods (Bihamta et al. 2015;Shishir and Tsuyuzaki 2018;Nasehi et al. 2019;Gaur et al. 2020). Various hybrid models have also been used to measure urban spatial pattern and future growth such as cellular automata (CA) and Markov chain (CA-Markov) models, CA-logistic regression (CA-LR) models, and Artificial Neural Network-Markov chain (ANN-MC) models (Feng et al. 2011;Behera et al. 2012;Gidey et al. 2017;Xu et al. 2019;Mohamed and Worku 2020). Among the different possible land simulation models, the ANN-MC model is a widely used model because it takes into consideration spatio-temporal elements of LULC, based on complex non-linear relationship among the variables (Tong and Liu 2005;Chakraborti et al. 2018;Xu et al. 2019). Compared to other predictive models, ANN-MC has a few advantages, such as a simple calibration method, self-adaptive, high prediction accuracy, and its ability to simulates multiple land covers in a complex multiland use urban environment (Paliwal and Kumar 2009). By combining both ANN-MC models, it may be possible to produce more accurate results in land predictions despite the complications of inappropriately parameterizing and optimally configuring of model (Gaur et al. 2020). The ANN-MC model is thus one of the most accepted and widely used methods for modelling land change simulation (Moghadam and Helbich 2013;Liu et al. 2018). Thus, in this study, we used an ANN-MC approach for LULC modelling to highlight the possible future state of UGS. Although, satellite-based remote sensing data have become popular for monitoring urban growth and associated issues (Tripathy and Kumar 2019;Vigneshwaran and Vasantha Kumar 2021), yet field investigation remains important to understand LULC change, which helps to visualize various factors or drivers of land-use change at the microscale and thus enhance the accuracy of modelling. In this study, both satellite data and field investigation data are used to identify the factors associated with LULC change and the loss of UGS.
Many developed countries have addressed the need for green space in urban areas and evaluated the vulnerability of UGS in response to urbanization (Gan et al. 2014;Nath et al. 2018). However, the situation is different in most developing countries where UGS has gradually been encroached by rapid and unplanned urbanization (Zhou et al. 2011). In some cases, the percentage of urban green spaces is below the minimum required standard of the World Health Organization (Asian Green City Index 2011; WHO 2012). In India, the conservation priority of UGS is not fully incorporated neither into urban planning nor in urban renewal programs, due to lack of city-wise detailed assessment (Gupta et al. 2012). The Ministry of Urban Development, Government of India has recently, emphasized the importance of green space conservation and development of the green growth framework through different urban renewal programs (MoUD 2014).
Considering the above methodological challenges, and the pressing need for UGS conservation in developing regions, the present study demonstrates an analysis of the future vulnerability of UGS, based on PSR-based indicators and an ANN-Markov model in an Indian city. The study responds to important research questions regarding the loss of urban green space of the city of Kolkata, including (a) what is the pattern of land-use change in the city; (b) how much green area has been lost; (c) which PSR indicators are the most responsible possible vulnerability of UGS and (d) what will be future UGS vulnerability in different land change scenarios? To address these research questions, the study examines LULC dynamics in the study area from 1980 to 2018, measures spatiotemporal changes in UGS in response to LULC change, predicts the state of LULC patterns and availability of UGS, and assesses the future vulnerability of UGS based on PSR indicators in the study area.

Study area
Kolkata is the capital city of the state of West Bengal and is situated along the east bank of the Hoogly river in the north-south direction (Figure 1). It covers an area of 185.52 km 2 with 141 wards i.e. smaller administrative units (KMC 2020). The city is under the jurisdiction of the Kolkata Metropolitan Corporation (KMC). Kolkata is the only city in India that has surpassed a complex colonial history and established itself as one of the largest cities in India. The city of Kolkata sits within the lower Gangetic delta with an average elevation of 1.5-9 m. The soil is alluvium in nature. The city lies under the tropical wet and a dry climatic environment with an annual mean temperature of 26.8 C and an annual average rainfall of about 185 cm.
Cities in India are incessantly growing by the process of population growth and urban expansion (Bhatta 2009;Ghosh et al. 2021); the megacity of Kolkata is no exception. Changes in LULC are rapid and impinging natural landscapes. Kolkata is one of the densest cities in India, with 24,306 persons per km 2 (Census 2011). Although population growth has slowed down and even became negative (-1.88%) between 2001 and 2011, the built-up area has kept expanding and UGS shrinking due to urban expansion and related infrastructure development. This trend is more likely to continue in the near future (Mondal et al. 2017;Dinda et al. 2021). Sharma et al. (2015) report that in Kolkata the built-up area has increased from 22.18% to 48.16% between 1989 and 2010, and, subsequently, the green cover has been reduced from 43.47% to 26.76%. In this viewpoint, the city is a suitable site for the assessment of LULC change and potential impacts on UGS.

Database
Multiple thematic data layers were used in the study, including satellite imagery, socioeconomic and environmental attributes, population data, location attributes, and environmental attributes, collected from different authentic sources. Multi-temporal Landsat satellite images were used for estimation of LULC change and UGS mapping because, Landsat provides a continuous series of data, it covers most of the area of the earth and is extensively used by researchers (Zhu et al. 2019). Images were acquired from the United States Geological Survey (USGS 2016). In comparison to the Sentinel data, although it provides finer resolution but data available since 2015. Detailed information about the satellite data used in the study is given in Table S1. Socio-economic data and population data were collected from the Census of India. Land value was collected from 'wbregistration.gov.in' website. Roads, railways, and other transport network data were collected from various published maps and recognized open data sources, including Open Street Map. The database and methodology outline are shown in Figure 2.

Methodology
2.3.1. Systematic approach for image pre-processing and LULC classification All images were pre-processed with ERDAS Imagine 2015 and other subsequent analyses were done with ArcGIS 10.7 software. The radiometric correction of the selected images was done by using ENVI 5.3 software. The Top of the Atmosphere process was adapted to convert digital numbers to radiance and reflectance (Hadjimitsis et al. 2010), which helps to reduce atmospheric attenuation. Suitable atmospheric correction was done before LULC mapping by using the Fast Line-of-Sight Atmospheric Analysis of Hypercubes described by Vermote et al. (1997). Geometric correction of the images was done by using 245 Ground Control Points. Among the various land classification methods, the supervised classification technique is a widely-used method for urban LULC mapping (Otukei and Blaschke 2010;Ghosh et al. 2019). In this study, the Maximum Likelihood Classifier (MLC) algorithm was used for the supervised image analysis and LULC mapping. The MLC decision rule was applied to reorganize seven LULC classes, including built-up area, vegetation and plantation, grassland and open space, agricultural land, water bodies, wetland and fallow land. Detail definition of these LULC classes was given in Table 1. In the present study, all types of urban green area, including urban forest, plantations, residential green, parks, organized green area, roadside trees and aquatic vegetation and grassland are defined together as urban green space (UGS).

Land transformation estimation using trend surface analysis
The post-classification comparison approach is very efficient for land transformation estimation of different LULC types and change detection analysis (Gaur et al. 2020). Spatial transformation of LULC maps were measured using Trend Surface Analysis (TSA) method. In order to assess the transition of each LULC class, the coordinates of the 'n' number of points (X i , Y i ) are taken as independent variables. It can be expressed as: Where Z is the variables of land transition between LULC types, X and Y represents the locational coordinates, and a, b, . . . , j are the polynomial coefficients. To produced TSA, the weight of pixels has been assigning to '1 0 for the transition scenario and '0 0 for no transition scenario. Moreover, to accurately identify the land-transition scenario, nineorder polynomial trend was employed (V aclav ık and Rogan 2009).

Accuracy assessment
Accuracy assessment was done to measure the scale of difference between classified images and reference data (Stehman 1996). Ground truth samples were selected from field survey sites and high-resolution Google Earth imagery for validation of the classification. Moreover, topographical map and maps collected from the KMC were used for the validation of classified images of previous periods. To evaluate the accuracy of the classified images, 245 random sample points were created in ArcGIS. Those points were selected based on a decision that each LULC class should cover at least 10% of samples (Congalton and Green 2009). The validity of selected points was checked by using GPS values during the field survey and corresponding Google Earth images. Error matrices were generated in the ArcGIS 10.7 software to calculate user accuracy (UA), producer accuracy (PA), overall accuracy (OA), and Kappa statistics. ANNs are widely used self-organizing and self-adapting machine learning technique for land-simulation (Azari et al. 2016;Rahman and Esha 2020). The ANN is based on forwarding structure black-box model, which is proficient by a back propagation (BP) algorithm (Xu et al. 2019). In this study, multilayer perceptron (MLP) based ANN architecture was employed for LULC simulation. MLP can adjust the weights between the variables based on mathematical relationship (Gaur et al. 2020). MLP-ANN is a forward neural network model with three layers, i.e. input, hidden and output .
The transition probability of each sub-model was constructed using 12 input neurons, 7 hidden neurons and 2 output neurons using the iteration limits to 10,000 and accepted error is 0.01. the optimal accuracy of each transition sub-model was above 85%, selected though trial and error process. The strength of association among the explanatory variables and urban expansion was measured using Cramer's V value. Generally, the Cramer's V value of 0.15 to 0.4 are good and the value greater than 0.4 are considered as excellent (Hamdy et al. 2017).

Calibration stage
Calibration stage helps to identify the location of change of LULC layers at base year (t) and (t )þ 1, using explanatory variables. In the present study, we have used four categories of variables, including proximity, socio-economic factors, neighbouring factors and resistance factors. All variables were selected based on previous studies ( (Table S2). After assessing the strength of the variables, simulation has been done.

Simulation using Markov chain
The Markov chain (MC) model is a stochastic process that efficiently predicts the conversion rates and transition probability from one land class to another (Mondal et al. 2017). The transition probability means a probability of future LULC change from one class to another which is represented by the 'transition probability matrix' (Nasehi et al. 2019;Dinda et al. 2021). Based on the Markovian process and Bayes theorem probability formula, the MC prediction function can express by: Where qij represents the transition probability of LULC change from class i to class j at the time of t and t þ 1, respectively; and n is the total number of land use classes in concern. We used the MC model to simulate LULC changes for six different periods, viz. 1980-1991, 1991-2001, 2001-2011, 2011-2018, 2018-2025 and 2025-2035. Finally, MLP-ANN model was combined with MC (MLP-ANN-MC) for LULC simulation and prediction analysis. The MLP-ANN-MC has higher ability to LULC simulate, shown in previous studies (Moghadam and Helbich 2013;Chakraborti et al. 2018).

Scenario development
To minimize the uncertainties of land-simulation factors, scenario development was important (Lu et al. 2016;Saeidi et al. 2018). The scenario set helps the LULC simulation model to be dynamic and recommend alternative land policies based on changes in different socio-economic factors (Goodarzi et al. 2017;Kim and Kwon 2021). In this study, two future scenarios were considered for the LULC simulation. The first scenario was the current land-change scenario, which was termed as the rapid-growth scenario, in which the present LULC change rate will be continued, and the population growth rate will persist as same as the present (Census 2011;Ghosh et al. 2021). The second scenario was the conservation scenario, in which the probability of UGS encroachment will be reduced by the promotion of green space conservation policies and the management of urban growth (MoUD 2014). Apart from this, it is also assumed that the population growth rate will reduce about 1% as compared to the present growth rate (UN 2019).

Model validation
The Markovian transition estimator tab was used for the assessment of model performance. The model calibration was performed by the comparative analysis between satellitegenerated actual LULC in 2018 and simulated LULC in 2018. Initially, after visual comparison analysis, chi-square (v2) statistics were performed for classification similarity assessment. To assess the areal distribution of the simulated land-use map, a more sophisticated Kappa index was used to validate LULC predictions (Pontius and Millones 2011). Kappa standard (K standard ), Kappa for no information (K no ), Kappa for location (K location ), and Kappa for quantity (K quantity ) was calculated. However, to measure classification agreement/disagreement, the values equal to '1 0 are considered perfect and imperfect, if the values equal '0 0 (Pontius and Millones 2011;Singh et al. 2015). Values more than 0.80 (>80% accuracy) indicate the reliability and efficiency of the simulation model ).
2.5. Pressure-state-response (PSR) framework development and derivation of weights 2.5.1. Developing of the index system Comprehensive, scientific, and feasible methods were adapted for the development of the indicators. Initially, 21 preliminary indicators have been selected based on a literature review (Chen et al. 2010;Gao et al. 2015;Bardhan et al. 2016;Li et al. 2016;Lu et al. 2016;Chakraborti et al. 2018;Munthali et al. 2020;Wu et al. 2020;Dinda et al. 2021;Ghosh et al. 2021) and finally, 15 indicators have been selected by proper screening method through field survey and experts' opinion. In order to understand the relationship between urban green space and associated risk factors, a panel of experts with relevant disciplines was engaged. Among them, seven members with significant knowledge of urban ecology, urban planning, and local urban biodiversity were invited. The average weight of each indicators was used for further processing. Indicators were further optimized/categorized under the PSR framework. The 'Pressure' (P) means the stresses originating from an adjoining area of UGS and it includes eight indicators. The 'State' (S) can be defined as the current status of the environment, including the quality and quantity of urban green space. 'Response' (R) indicates the results of the impacts and the kind of actions taken for management. Six indicators were selected for 'Pressure' and three indicators were for 'Response' parameters. The selected indicators, their significance, and their sources are mentioned in Table 2.

Determination of indicator weights using Analytic Hierarchy Process (AHP)
Analytic Hierarchy Process (AHP) was used to calculate the weights of the selected indicators which are responsible for UGS vulnerability. AHP is a widely used multi-criteria decision-making method (Ghosh et al. 2018;Sun et al. 2019;Banerjee et al. 2020). It consists of three steps, including the establishment of hierarchical structure, constructs a pairwise comparison matrix, and consistency check. The relative importance of indicators was selected based on a hierarchical structure with an appropriate scale ranging from 1 to 9 and their reciprocal (Saaty 2008). Relative importance, i.e. the score of the indicators was assigned according to the experts' opinion, invited from the field of environment, urban planning departments. The pairwise comparison matrix was formed using the mean score, denoted as A: Each entry (a ij ) of the comparison matrix (a ij > 0Þ, expresses the experts' score of the comparison of the criteria i-th with the criteria j-th. It is worth expressing that a ij ¼ 1, and a ij ¼ 1=a ji , whenever, i ¼ j and a ij 6 ¼ 0: The eigenvector is designated as w, matching the largest eigenvalue, lambda max (k max ) of the matrix is the derived solutions of the model. It can be expressed as (Banerjee et al. 2020): In another form, Here eigenvector should be P n i¼1 w i ¼ 1 and k max ¼ n, under the ideal condition. n denotes the number of indicators. The reliability of the AHP should be assessed because the relative importance of different indicators is varying (Ghosh et al. 2018;Yu et al. 2019). Consistency ratio (CR) determines the reliability of the model, which is calculated from the ratio of consistency index (CI) and random consistency index (RI). It can be expressed as: Where CI ¼ ðk max À nÞ=ðn À 1Þ and RI value was obtained from Saaty (1980). CI value of less than 0.1 is acceptable for an AHP model. In this study, Expert Choice software was used to conduct the test of the model.

Fuzzy Comprehensive Evaluation (FCE) model
FCE model is a combined model of fuzzy set theory and mathematical modelling (Chiu et al. 2014). It is a very effective model to manage the inter-relationship and interdependence problem of multi-criteria decision-making methods (Liu et al. 2014;Shao et al. 2020). FCE consists of three steps: build the factor set, allocate an evaluation set, and finally developed the FCE set (Sun et al. 2019). First, the parameter index (I) was used to construct the factor set. It can be represented as, I ¼ ði 1 , i 2 , . . . , . . . , i n Þ, where i is the original value of the indicators. Then the indicators of the parameter were divided into five equal classes from 0 to 1 with an interval of 0.2 (Zhu et al. 2017;Sun et al. 2019). It can be expressed as C ¼ ðc 1 , c 2 , c 3 , c 4 , c 5 Þ, where C is the allocated evaluation sets of five selected classes representing, 'very high', 'high', 'medium', 'low', and 'very low'. Therefore, the result of the AHP model was evaluated under FCE evaluation using the following matrix (Sun et al. 2019 Where bx are the FCE sets, b is the weighting coefficient; a 1 , a 2 , . . . , . . . , a m are the AHP derived weights, x is the fuzzy transformation matrix; t m1 , t m2 , t mn are the score of allocated evaluation sets of five equal grades, and b 1 , b 2 , . . . , . . . , b n are the FCE values, respectively. The score of allocated evaluation sets was normalized using the following equation: Where l is the normalized value, x i is the expert opinion-based score of i-th indicator; x min and x max are the lower and upper limit of each indicator, respectively.

Future vulnerable zones assessment of UGS
Future vulnerable zones aim to insert proper management strategies to conserve UGS and to formulate policy accordingly. Based on PSR indicators and CA-Markov predicted LULC data, we have developed the present and future vulnerable zones of UGS for Kolkata. The weighted linear combination (WLC) model was adopted to produce future vulnerable zones of UGS using a raster layer of risk factors and their criteria weights derived from AHP-FCE. The 'union' operator function of ArcGIS 10.7 was used for this operation (Nath et al. 2020). The future vulnerable zones can be expressed by the following formula: Where FVZ represents the future vulnerable zones of UGS, f n represents the UGS vulnerability indicators and b n are the FCE normalized weights of the indicators.

Spatio-temporal change analysis and accuracy assessment (1980-2018)
The built-up area increased rapidly from 92.71 km 2 in 1980 to 148.30 km 2 in 2018 and occupied more than 80% of the total area of the city (Figure 3). Being a metropolitan and capital city, Kolkata has experienced ceaseless urban expansion (Bardhan et al. 2016;Sahana et al. 2018;Ghosh et al. 2019). The results revealed that the built-up areas are rapidly expanding towards the outskirt of the city, following the major transportation axes. Built-up land expanded faster since 2001, indicating an accelerated urban sprawl than identified in previous studies (Bhatta 2009;Mondal et al. 2017;Dinda et al. 2021). The percentage change of built-up land was 17.23% between 1980 and 1991, but gradually increased to 21.13% between 2011 and 2018. In contrast, other LULC areas are gradually decreasing due to built-up expansion and encroachment ). The overall change in the vegetation cover and grassland area from 1980 to 2018 was (-) 2.48% and (-) 1.84%, respectively. The total area of the vegetation cover decreased from 45.02 km 2 to 13.75 km 2 from 1980 to 2018. All classified LULC maps were verified using Kappa statistics, user accuracy, producer accuracy, and overall accuracy (Figure 4). Overall accuracies of the classified maps were >85% in all the base year. The Kappa values of the classified maps were >0.85, reached at a satisfactory level (!0.85) for further analysis (Monserud and Leemans 1992).

Validation analysis and LULC change simulation (2018-2035)
The association among the explanatory variables and land-use change was measured using the Cramer's V value for all selected time intervals (Table 3). To validate the LULC simulation model, the simulated LULC map of 2018 was compared with the actual LULC map of 2018. Initially, after the visual similarity assessment, the chi-square test was performed to compare and validate the area of each LULC type. For this purpose, we considered that the LULC areas of the simulated map and actual map are to be the observed (O) and expected (E) LULC areas, respectively (Figure 5a). The ROC curve of the model was shown in Figure 5b, represents the level of efficiency of the simulation model. It was found that the tabulated value (v 2 0:05, ð8Þ ¼ 15:5) was greater than the value of v 2 (0.18); hence, we reject the alternative hypothesis. Results indicate that the adapted model was a good fit to run the land simulation and future prediction in the study area (Table S3). All Kappa coefficient values ranged from 0.86 to 0.92, which can be described as 'almost perfect agreement' (Congalton and Green 2009). The overall agreement between actual and simulated maps of 2018 is 0.91, while the error due to quantity (disagree quantity) is only 3.06% and the error due to allocation (disagree grid cell) is only 3.86% (Table S4). Therefore, from these indices, it can be inferred that the model is efficient to simulate LULC for 2025 and 2035. The simulation result shows that the built-up area will persist to expend up to 2.00% to 2.34% and 1.16% to 3.77% between the period of 2018-2025 and 2025-2035, respectively under the current land-change scenario. The simulated results indicate that built-up land area will be expanded up to 155.67 km 2 of total land cover (185.52 km 2 ). It was also found that a major decrease will be observed in four LULC classes, including vegetation cover, grassland, water bodies, and, wetlands. The main objective of this study was to assess UGS changes and likewise, the green spaces that will be severely reduced by the process of urban expansion. In the current land-change scenario (rapid growth scenario),  the rate of decline of vegetation cover will be up to (-) 11.78% and (-) 17.42% in the period 2018-2025 and 2025-2035, with a yearly decreasing rate of (-) 1.68% and (-) 17.43%, respectively. However, in the conservation scenario, little development will be observed in vegetation and grassland with a diminutive sharing in the area in comparison to the total study area. From the transition probability matrix of 2018 to 2025, it was found that vegetation and grassland have more prospects to change, around 63.79% (p ¼ 0.63) and 53.15% (p ¼ 0.46), respectively. Higher rates of change were seen for bare land, cropland, and water bodies. Similarly, in the scenario for 2018 to 2035 land transition, vegetation and grassland will be gradually reduced and converted to other land-use classes with a probability of 87.09% (p ¼ 0.12) and 76.01% (p ¼ 0.23), respectively.

Overall trend analysis of LULC change
The LULC change analysis and simulation results predict an exceedingly rapid increase in the proportion of land under built-up area and a fast decrease in the case of the vegetation cover and grassland. Change of different LULC areas and their spatio-temporal trend was depicted with the help of Trend Surface Analysis (Figure 6). The growth rate in built-up land was 31.77% and 95.49% from 1980 to 2018. In comparison to the built-up area, the other categories of LULC classes are expected to deteriorate. The overall change in UGS from 1980 to 2035 was observed to be (-) 79.47%. The UGS occupied 45.02 km 2 in 1980, which decreased to 13.75 km 2 in 2018 and is expected to further decrease to 10.21 km 2 in 2035. Hence, the peripheral area of the city, particularly the southern and south-eastern part of the city are shown an increasing trend of built-up area.

Weights of the PSR indicators and analysis of UGS vulnerability factors
Using the AHP-FCE method, weights of the green space vulnerable indicators were estimated. Assessment of the PSR indicators was based on base-level data, i.e. 2018, and the simulated LULC data, i.e. 2035 which were further classified as 'rapid growth scenario' and 'conservation scenario' through the scenario setting. The consistency ratio of the model was 0.006, 0.010, and 0.011 of 2018, 2035-with 'rapid growth scenario', and 2035with 'conservation scenario', respectively (Table 4). Therefore, the models were within acceptable limits which allow for further analysis. The overall weights of 'pressure', 'state', and 'response' indicators in 2018 were 0.307, 0.197, and 0.102, respectively, which was predicted to be increased by 0.411, 0.309, and 0.195, respectively under the rapid growth scenario in 2035 (Table 5). Thus, it was found that 'pressure' indicators had a higher impact on the vulnerability of UGS. Among them, proximity of built-up area (0.112), proximity of major and minor roads (0.109), distance from the city core area (0.099), proximity of metro stations (0.098) had the higher weights. Similarly, these same factors were largely responsible for the analysis of the UGS vulnerability in 2035 under the rapid growth scenario. The weights of indicators, including proximity of built-up area (0.277), proximity of major and minor roads (0.199), distance from the city core area (0.156), proximity of railway (0.109), and the population density (0.188) had increased. Spatial layers of the PSR indicators were shown in Figure 7.
The built-up growth rate (0.109), land value (0.071), and proportion of open space (0.078) were the most important factors under the 'state' category for further urban expansion and responsible for green space vulnerability. Similarly, the weight of the 'state' indicators had increased during the assessment of the AHP model in 2035 under the rapid growth scenario. Land use dynamic degree and green space change intensity index were the most weighted indicator of the 'response' category. It is further observed that the weight of these indicators will increase further in 2035 under the scenario of rapid growth. However, it was found that the weights of the vulnerability indicators will be reduced in the conservation scenario of 2035, as compared to the rapid growth scenario of 2035. The green space management plans indicator in the conservation scenario of 2035 showed the largest weight of all the 'response' indicators. The spatio-temporal dynamics of 'Pressure', 'State', and 'Response' indicators are shown in Figure 8.

Future vulnerable zone identification and analysis
As stated above, the total available green area of the city was only 15.58% in 2018 and it was found to likely further reduce by 29.11% in 2035. Hence, the simulation of green space vulnerable zones is a helpful instrument for the management of UGS Kim and Kwon 2021) The future vulnerable zones can contribute to understanding the undergoing stress in the green areas and the risk of future encroachment. Accordingly, we have combined the AHP-FCE method with the MLP-ANN-MC produced LULC data to generate three future vulnerability zones of UGS in Kolkata. The study exhibits that the vulnerable areas will be gradually expanded in the period 2018 to 2035. However, the map shows that in the case of the conservation scenario, the areas covered by the 'very high' and 'high' categories will be slightly reduced. The MLP-ANN-MC based simulated results reflect the rapid growth of the low-density built-up area due to urban sprawl, particularly in the urban fringe area of the city Ghosh et al. 2019). Hence, the green space in the outskirt areas will be more vulnerable due to the high growth rates of urban growth. Figure 9 shows that the southern, south-eastern, and south-western part of the city will be more prone to vulnerability if the rate of urban expansion continues.

Analysis of the role of green space vulnerability indicators
Understanding the future Spatio-temporal pattern of LULC is one of the effective pathways to evidence-based management of natural environments (Munthali et al. 2020). In this study, the the main objective was to assess the vulnerability assessment of green space due to changes in LULC using PSR framework-based indicators and highlighting its present and future scenarios. Kolkata city has experienced a high rate of urbanization since the colonial period. Thus, urban land-use drastically changed and resultant modification has been found in the natural landscape . The present study was based on a combined GIS-based method and field survey to analyze LULC change. Initially, LULC classification was done using Landsat satellite imageries. The MLP-ANN-MC was adopted for LULC simulation in 2025 and 2035 using a series of LULC data, driving factors, and suitability layers. The result shows that the built-up areas increased gradually by encroachment towards green areas, i.e. vegetation, and grassland area. Between 1980 and 2018, the total loss of UGS was more than 69.47% of its total area. It mainly occurred due to the encroachment of the built-up area at a tremendous rate (95.49% between 1980 and 2018) . This study presented a multi-dimensional PSR based analysis using an integrated AHP and FCE model to estimate the weights of the UGS vulnerability indicators. The outcome of the model revealed that factors, including nearness of built-up area, roads, railway, city core area, and population density belong to the 'pressure' category are mainly responsible for UGS vulnerability in Kolkata. Ghosh et al. (2021) have mentioned all these factors in determining the urban ecological security of Kolkata. The land use dynamic degree index of the built-up area and green space change intensity index was negatively correlated for Kolkata, reported by Dinda et al. (2021). It is worth noting that, dynamic degree index of the built-up area progressively increased from 1.57% to 2.02% between 1980-1991 and 2011-2018 and it was simulated that the same will be increased up to 3.17% in 2035 ( Figure 10). The value of the green space change intensity index was decreased gradually from (-) 0.14% to (-) 0.74% between 1980-1991 and 2011-2018. The projected value of the change intensity index also exemplifies a negative change in UGS due to the percentage increase in other LULC classes. Hence, the AHP-FCE model has produced higher weights of these two indicators.
The integrated MLP-ANN-MC model was found to perform well in this study in comparison to previous research. Mondal et al. (2017) applied the CA-Markov model to predict the land suitability surface of Kolkata, reaching 81% accuracy, whereas the accuracy achieved by the present study is about 91%. LULC statistics exhibit that vegetation and grassland were two important land components that share about 41.15% of the total area in 1980, which reduced 15.58% of the total area in 2018 by the expansion and encroachment of built-up areas. We found several factors, including proximity to roads, railway, metro, availability of open land, land price, and other infrastructure opportunities to be major factors contributing to the growth of built-up land and urban sprawl. Previous studies have also found that vegetation cover, grassland area, and open land provide a suitable/potential surface for built-up growth (Bardhan et al. 2016;Sahana et al. 2018). The analysis of the dynamic degree index and change intensity index provides a state of the relationship between built-up area growth and associated changes in green space ). However, unabated urbanization not only reduces the green areas of a city, but can also deteriorate ecologically sensitive areas (Wu et al. 2020). Therefore, the future vulnerable zones of UGS will be evidence of green space stress and level of ecological sensitivity, because criteria-based vulnerability modelling are useful to replicate spatial characteristics of UGS changes under feasible land management practices (WEF 2015;Ghosh et al. 2021) (Figure 11).

Advantages and limitations of the model, and scope of future work
Kolkata has experienced rapid urban expansion, as mentioned in previous studies (Sharma et al. 2015;Mondal et al. 2017;Sahana et al. 2018;Ghosh et al. 2019;Dinda et al. 2021). In this study, we have used a land-use simulation model to assess the UGS vulnerable zones with two different scenarios, i.e. rapid growth scenario and conservation scenario. Moreover, the status of UGS was clearly demonstrated in the LULC change scenario and TSA. Besides, the present study has focussed on the PSR-based vulnerability of UGS, rather than only land-use change or patterns. It is worth mentioning that, AHP-FCE has been widely used in various types of environmental research, which is more efficient and more robust (Sun et al. 2019;Banerjee et al. 2020;Shao et al. 2020). Furthermore, GIS-based maps were produced using the extracted indicators' weight to  comprehensively analyse the vulnerability scenario. In addition, this study reinforced the fact that rapid growth of the built-up area, land-use dynamic degree, proximity of roads, and railway are the major responsible for the UGS vulnerability.
However, the present study does have a few limitations. First, the analysis of LULC maps was scale-dependent as the derived remote sensing database had a 30 m spatial resolution. Second, the AHP is a purely expert opinion-based model; therefore, the opinion can be varying. To minimize this error, the assigned weights by the experts were normalized into non-dimensional values using the FCE method. As a result, the weights were more reliable for the entire analysis. Third, the assimilation of various landscape matrices would be better able to measure the distribution of UGS and vulnerability levels. Therefore, the application of landscape matrices and the quality of UGS requires further research. Yet, the current study is probably the first initiative for Kolkata city to assess the  vulnerability of UGS. The future vulnerable zone of UGS will help to accentuate the green space sensitivity area of the city. This study might provide new dimensions to assess the UGS and, in particular, help redesign city greening policies.

Conclusion
This study has presented a straightforward and replicable framework for the assessment of UGS vulnerability that can assist policymakers in land management in a city. We simulated LULC changes using the integrated MLP-ANN-MC model and analyzed associated loss of UGS. The simulated LULC changes allowed the assessment of impacts of rapid urban growth in UGS availability and showed its future states. The result shows that the landscape of Kolkata was highly modified during different periods and green space was the most significantly altered land class by the process of urban development. The loss of UGS poses major challenges to urban planners and governments in the achievement of sustainable urban growth. Results indicate a decreasing trend of UGS from 1980 to 2018 and likely further decline in the future. The assessment model performed well in terms of accuracy assessment. The vulnerable area of UGS was assessed by future vulnerable zone mapping, based on PSR-based indicators. It can help to adapt and implement management policies. The findings of this study can help to understand the state of UGS and its vulnerability, which is one of the important tasks for the urban local body (Anguluri and Narayanan 2017).
Contemporary environmental problems such as air pollution and urban heat island, which cause various health issues, have increased the requirements of green space in a city (Kim and Kwon 2021;Dutta et al. 2020). Dinda et al. (2021) have reported that percapita UGS of Kolkata was only 7.2 meter 2 per person in 2018, which is below the minimum requirement for UGS, as recommended by the World Health Organization. The findings of the study also suggest that the increase of UGS vulnerable areas will be reduced in the conservation urban growth scenario, as compared to the rapid growth scenario. Thus, policies should be implemented by following the 'Urban and Regional Development Plans Formulation and Implementation Guidelines' (MoUD 2014) to protect the existing green space and also took a plan to increase the areas of UGS to fulfil the minimum requirement. Moreover, future research is needed to investigate the UGS suitability zones by considering the demand for various land uses and feasible land management practices.
However, the academic contributions of the present study are: (a) responsible factors for the UGS vulnerability were selected by incorporating the literature review, expert opinion, and field survey, and categorized under PSR framework, (b) incorporating expert opinion based AHP with FCE method, make the model more robust, and (c) prediction of UGS vulnerability with controlled and uncontrolled land-growth scenarios has provided a comprehensive framework for assessment. To conclude, this study will help the town planners and policymakers to understand the possible risk of UGS due to the continuous urban growth process, and develop policies, accordingly.

Declarations
There has no specific funding particularly for this work but Dinda, S, and Ghosh, S is willing to thank the University Grant Commission (UGC), India for providing the research fellowship.