A Novel Technique for Modeling Ecosystem Health Condition: A Case Study in Saudi Arabia

: The present paper proposes a novel fuzzy-VORS (vigor, organization, resilience, ecosystem services) model by integrating fuzzy logic and a VORS model to predict ecosystem health conditions in Abha city of Saudi Arabia from the past to the future. In this study, a support vector machine (SVM) classiﬁer was utilized to classify the land use land cover (LULC) maps for 1990, 2000, and 2018. The LULCs dynamics in 1990–2000, 2000–2018, and 1990–2018 were computed using delta ( ∆ ) change and Markovian transitional probability matrix. The future LULC map for 2028 was predicted using the artiﬁcial neural network-cellular automata model (ANN-CA). The machine learning algorithms, such as random forest (RF), classiﬁcation and regression tree (CART), and probability distribution function (PDF) were utilized to perform sensitivity analysis. Pearson’s correlation technique was used to explore the correlation between the predicted models and their driving variables. The ecosystem health conditions for 1990–2028 were predicted by integrating the fuzzy inference system with the VORS model. The results of LULC maps showed that urban areas increased by 334.4% between 1990 and 2018. Except for dense vegetation, all the natural resources and generated ecosystem services have been decreased signiﬁcantly due to the rapid and continuous urbanization process. A future LULC map (2028) showed that the built-up area would be 343.72 km 2 . The new urban area in 2028 would be 169 km 2 . All techniques for sensitivity analysis showed that proximity to urban areas, vegetation, and scrubland are highly sensitive to land suitability models to simulate and predict LULC maps of 2018 and 2028. Global sensitivity analysis showed that fragmentation or organization was the most sensitive parameter for ecosystem health conditions.


Introduction
With economic prosperity, demographic growth, and LULC transition, urbanization is the primary driver of social and ecological transformation [1]. More than half of the world's population now lives in towns, and the rate of urbanization is expected to rise to 68% by the mid-century period [2]. The flow of material, resources, structure, and function of ecosystems can all be dramatically altered by urbanization [3]. For sustainable socio-economic growth, the growing population and urbanization necessitate higher ecological services [4], which leads to the emergence of a commonly affecting nexus between urbanization and ecology [3]. For example, urban expansion has resulted in the conversion of land fragmentation and LULC patterns from initial to built land, which has a significant impact on ecosystems, such as habitat loss, decreases in agro-forest productivity, and the reduction of plants' climate-regulating role [5]. Consequently, continuous ecosystems monitoring and impact assessments on urban ecosystems have become a scientific theme for sustainable urban landscape planning and environmental policy implication.
Ecosystem health (EH) has been widely used in a variety of areas, including urban, forest, grassland, and marine habitats, as a critical state of eco-environmental development [6,7]. Rapport [8] and Costanza [7] presented the idea of EH, which was described as an ecosystem's ability to preserve its actual function and original structure under external pressure. This concept of EH facilitates ecological studies and is well accepted by many academics [9,10]. The definition is also based on a philosophy, namely urban ecosystem health (UEH), which supports urban environmental managers to connect environmental, socio-economic, and human health drivers [11]. UEH encompasses not just the health and integration of natural and built landscapes but also the health and well-being of urban dwellers and the whole population [12]. The UEH has been a topic of interest and a major target in environmental planning, demonstrating how people perceive satisfaction with the ecosystem's capacity to restore itself. Furthermore, proper EH monitoring and quantification of the impact of urbanization on ecosystems are critical for LULC planning and management. Our research topic was inspired by widespread public interest and decades of progress in EH research [12].
Despite the complexities of indirect EH measurement, attempts have been made to quantify EH using a common method of developing an indices system based on multiple indicators [9,10]. To measure the function and service of EH in the field of the urban ecosystem, single or multiple indices with an individual indicator, e.g., landscape matrices, three indicators, e.g., vigor-organization-resilience (VOR), or four indicators e.g., vigor-organization-resilience-ecosystem service (VORS) models are used [2,13]. However, previous studies have mostly focused on the overall performance of EH, with only a few studies effectively applying the LULC transition and the transformation of landscape distribution at the local and regional scales [14,15]. Furthermore, although efforts to demonstrate EH change are ongoing, a deeper understanding of the impact of urbanization is still unknown. Appraising and quantifying EH from the perspective of LULC transition and landscape fragmentation pattern can provide a spatial view of the dynamics and eco-environmental consequences of UEH growth under certain conditions. Many scholars have focused on the relationship between urbanization and EH. For example, Xiao et al. [6] measured mountainous EH in southwest China based on ecosystem service importance and a set of metrics and discovered that urbanization has a negative relationship with EH. Wang et al. [9] used a selected series of metrics to investigate the effect of urbanization on EH in Zhuhai, China and concluded that the city has experienced ecological degradation since 1999, with a negative correlation between urbanization and EH. Cheng et al. [16] evaluated EH in the Haihe River Basin, China, using several indicators and found that arable land, gross domestic product (GDP) per capita, and population density had inverse relationships with river EH. In a data-scarce environment, Van Niekerk et al. [17] used five indices to investigate the key stresses from urbanization and EH. Styers et al. [18] computed EH using landscape indices and urbanization parameters. It is crucial to maintain EH at the top level for adequate resources and ecological development. However, evaluation of EH in urban environments results in a lack of coupling and specific research on fast-growing urban cities, specifically on the spatial heterogeneity of EH and its sensitivity to multiple human-induced pressures. As a result, there is an urgent need to research how to map the spatial heterogeneity of EH and its sensitivity using integrated and quantitative techniques. The study area, Abha district in the Kingdom of Saudi Arabia (KSA), is being investigated because it has undergone sustained and rapid urbanization. It has been chosen as one of the most important cities in the KSA from an economic standpoint. The rapid economic development and rapid urbanization have led to significant changes in the functioning and composition of the ecosystems across Abha city, Saudi Arabia, which have resulted in the reduction of the EHs [19]. In recent decades, Saudi Arabia's urban population has grown at the expense of arable land. This city has experienced rapid urbanization, leading to a severe decline in eco-environmental sustainability and the loss of the natural environment. Nevertheless, so far, to the best of the author's knowledge, no quantitative appraisal has been undertaken to evaluate the effect of urbanization on the EH in Abha city, and also the sensitivity of the EH trends and their spatial patterns remain vague. Therefore, the specific objectives of our study are to (1) map LULC from 1990 to 2018 using the support vector machine classifier; (2) predict LULC for 2028 using the CA-ANN model and analyze the sensitivity using RF, CART, probability density function (PDF), and Pearson's correlation techniques; (3) develop the fuzzy-based VORS model for building ecosystem health conditions from 1990 to 2028; and (4) investigate the performance of the EH model using the Morris technique-based global sensitivity analysis.
As compared to previous research, our research is novel in two ways: technical and objective aspects. With regard to technical issues, the implementation of SVM for LULC mapping in Saudi Arabia and developing countries is not yet a widespread endeavor. The use of CA-ANN for forecasting LULC is relatively new in developing countries, although the use of machine learning algorithms in conjunction with PDF and correlation coefficient to investigate the sensitivity of CA-ANN models is unique in terms of measuring the effect of input parameters for CA-ANN models. Furthermore, we incorporated fuzzy logic knowledge with the VORS model to improve the accuracy of EH models, which can be considered the unique work in the domain. Finally, since ground truth evidence for validating the EH models is unavailable and difficult to obtain, the global sensitivity analysis (Morris method) was used to make the models reliable to stakeholders and the scientific community. This is the first research to use global sensitivity analysis to investigate the performance of EH models. As a matter of fact, in terms of technical aspects, this research is very unique. In terms of objective, on the other hand, the present work offers detailed analysis on past to present LULC mapping and dynamics using advanced methods, a new approach to present EH conditions for past to future period, and uncertainty estimation. As a result, to date, this is the first study that considers various factors at the local scale and provides benchmark references for future urban eco-environmental planning.

Study Area
Abha is a city in Saudi Arabia's southwestern Asir province ( Figure 1). This city's topography is rugged and heterogeneous, with elevations ranging from 1950 to 2982 m above sea level. The city is characterized by a cold and semi-arid climate. The average annual precipitation is 355 mm and occurs mainly from June to October. The average annual minimum and high temperatures were 19.3 • C and 29.70 • C, respectively. This city faces extreme difficulties due to land degradation caused by man-made activities, higher slopes, weaker geological formations, and unpredictable rainfall, which affects ecosystem services (ESs). Plant communities are thriving as a result of climate change and diverse topography. This city is one of the Asir Mountains' most diverse floristic regions. The urban hills are the most appealing tourism destination in this region, with the most diverse flora and fauna.

Materials
The Landsat 4-5 TM and 8 OLI data (path/row: 138/44 and 138/45, spatial resolution: 30 m) for 1990, 2000, and 2018 were obtained from the USGS Earth Explorer website (https://earthexplorer.usgs.gov (accessed on 21 November 2019)). To derive topography parameters, an ALOS PALSAR radiometrically terrain corrected (RTC) Digital Elevation Model (DEM) with 12.5 m spatial resolution was obtained from Earth Science Data Systems NASA. These satellite images were collected between October and January in order to achieve a credible and cloud-free dataset. As a basic pre-process, atmospheric correction was performed on the above datasets.

Methods for Image Classification and Validation
A non-parametric machine learning algorithm (MLA) called the support vector machine (SVM) classifies data based on statistical learning theory [20] and has been used in the present study for classifying the Landsat images for 1990, 2000, and 2018. In the field of remote sensing, the SVM model has been successfully used for data classification

Materials
The Landsat 4-5 TM and 8 OLI data (path/row: 138/44 and 138/45, spatial resolution: 30 m) for 1990, 2000, and 2018 were obtained from the USGS Earth Explorer website (https://earthexplorer.usgs.gov (accessed on 21 November 2019)). To derive topography parameters, an ALOS PALSAR radiometrically terrain corrected (RTC) Digital Elevation Model (DEM) with 12.5 m spatial resolution was obtained from Earth Science Data Systems NASA. These satellite images were collected between October and January in order to achieve a credible and cloud-free dataset. As a basic pre-process, atmospheric correction was performed on the above datasets.

Methods for Image Classification and Validation
A non-parametric machine learning algorithm (MLA) called the support vector machine (SVM) classifies data based on statistical learning theory [20] and has been used in the present study for classifying the Landsat images for 1990, 2000, and 2018. In the field of remote sensing, the SVM model has been successfully used for data classification [20]. In data mining, SVM is a common classification algorithm. It makes use of decision hyper-Remote Sens. 2021, 13, 2632 5 of 24 planes, which establish decision boundaries in the input space or high-dimensional feature space. As a result, it is more suited to binary categorization. However, most datasets in remote sensing include many classifications. To classify datasets with several classes, multi-class SVM is employed. The kernels and parameters used in SVM kernels have a significant influence on the accuracy of SVM classification. The appropriate selection of two parameters, C and γ, has a major impact on the SVM's performance. The cost of misclassification is denoted by C. Since just a small number of points are chosen as support vectors to define the decision boundary, a low value of C leads to under-fitting. A high value of C, on the other hand, leads to over-fitting, since more data samples are picked as support vectors, resulting in a wider decision boundary. The ideal value of C is the one for which the training error is the smallest and the SVM's generalization capacity is the greatest. The distance between the decision border and the nearest support vectors is represented by the parameter γ. Low γ values cover points that are far away from the decision border, whereas high γ values only cover the ones that are closest. As a result, for defining the decision boundary, an ideal value of is necessary. The grid search technique is the standard way for picking these parameters, although it is time intensive. Furthermore, the SVM's performance is influenced by the kernel selected.
We employed a radial basis function (RBF) for the classification of multi-class LULC datasets to address these limitations. As a result of its high generalisation capability, the RBF kernel is often utilized. Equation (1) is the formula that defines it.
where γ = 1/2σ 2 (σ is a positive parameter for regulating the radius). The distance between the hyperplane and adjacent data points is represented by γ. The mistake is minimized with a wider margin of γ error. Changing the width σ has an effect on γ. The RBF kernel is a local kernel in which data points adjacent to each other impact kernel points. Choosing the best values for γ and C is a tough task that is still being researched. When the RBF kernel is employed, the right values of γ and C improve the accuracy of the SVM classifier.
The SVM in the present study was run with the following parameters and kernel RBF: a kernel width gamma (γ) of 0.143, a penalty parameter (C) of 100, and a classification probability threshold of 0.3. For classifying the LULC for the years 1990, 2000, and 2018, we used five spectral bands, except for the thermal band and pan band for Landsat 4-5TM (years 1990, and 2000), and six spectral bands for Landsat 8OLI (except for band 1,[8][9][10][11]. The training and testing data for 2018 have been generated from field survey and Google Earth images; on the other hand, the training and testing data for the years 1990 and 2000 were collected from Google Earth and topographical maps. A total of 300 sample sites were selected to be training and testing samples. Then, we divided the whole dataset into training and testing datasets as per an 80:20 ratio.
Validation of the prepared LULC map is a must before it can be used for further research. As a result, in this study, LULC maps for 1990, 2000, and 2018 were validated using ground truth and secondary sources. The LULC map of 2018 was validated using a field survey, whereas the LULC maps of 1990 and 2000 were not validated using a field survey or a Google Earth image. As a result, we had to rely on a secondary source. Landsat 4-5TM blue, shortwave infrared, near-infrared bands, and topographical maps were used for testing between 1990 and 2000. Based on the testing samples, we extracted information of different land use features through visualization and compared it with topographical maps. Based on the collected samples, we compared the samples with prepared LULC maps of 1990 and 2000. In statistical validation of LULC maps, the Kappa coefficient and RMSE were used. On the other hand, the RMSE was calculated based on the correct LULC samples in the ground truth and classified LULC maps.

CA-ANN Model
In this study, a MOLUSCE QGIS plugin was used to predict future land coverage. The artificial neural network (ANN) method was used in the cellular automata (CA) model because it has been shown to produce better results than other techniques such as linear regression [21] (for more details, see the method section of the Supplementary File). In the present study, CA-ANN was used to predict the LULC map of 2028.

Sensitivity Analysis
Since future LULC cannot be validated, the model's performance should be tested for reliability and robustness. As a result, sensitivity analysis has been used in this research. We have used RF, CART, Pearson's coefficient of correlation, and probability distribution function techniques (for details of the methods, please see the method part in the Supplementary Section). For conducting sensitivity analysis, we extracted information from these variables and final models such as simulated LULC, predicted LULC of 2028, and four ecosystem health models, based on 3000 randomly chosen points using 'create random points' and 'extract multi values to points' tools in ArcGIS 10.5 software. Then, the sensitivity models have been used on these collected data.

Ecosystem Health Assessment Framework
To assess Abha City's EH, we used the VORS model, which includes four indicators: vigor, organization, resilience, and ecosystem services (ES). The ES is assessed by two criteria, namely the status of the ecosystem physical health, including three indicators [22,23] and the ES of Abha City. A high ES value indicates that the community is receiving sufficient material and resources, and the ecosystem is in a very stable and balanced condition. We hypothesized that the ES and EH levels are interrelated positively. In our study, EH values are classified into five categories: very good (80-100%), good (60-80%), moderate (40-60%), very poor (20-40%), and poor (0-20%) [2]. Equations (2) and (3) for calculating EH are as follows [2,24]: where H denotes the EH level of the units being appraised in Abha City and PH denotes the ecosystem's physical health. ES is the ecosystem service value. However, the input parameters were standardized according to the positive and negative characterization before modeling EH conditions. As a result, the standardized parameters were varied from 0 to 1. Nonetheless, researchers have used a variety of expertbased and computational methods [25]. However, we used a fuzzy membership function to make the standardization process more robust. The fuzzy membership function generates fuzzy crisp layers depending on the layer's characteristics. The MS small membership function was used for the fragmentation of organization parameters in this analysis, while the MS large membership function was used for V, R, and ES. On the other hand, the fragmentation parameters were calculated using a fuzzy membership function and fuzzy operators. The specifics were given in the ecosystem organization (O) section. The EH conditions were predicted using Equations (1) and (2) based on the fuzzified parameters.
In general, ecosystem vigor refers to the input, metabolism, or primary productivity of an ecosystem [7]. The NDVI was used to quantify V in our study, since it denotes the extent of vegetation dynamics and primary productivity of vegetation [26]. The NDVI value ranges from −1 to +1. The NDVI value close to zero indicates no vegetation growth, while the value close to one indicates the greatest concentration of green vegetation.
The number and diversity of relationships among the components of an ecosystem are referred to as ecosystem organization (O). Ecosystem O was used to assess the spatial heterogeneity of landscape patterns and landscape connectivity in this study. Spatial distribution patterns are well known to be the most influential factors in controlling Remote Sens. 2021, 13, 2632 7 of 24 ecological processes at the landscape level. Shannon's Diversity Index (SHDI), Patch Density (PD), and the Area-Weighted Mean Patch Fractal Dimension (AWMPFD) were used to assess landscape heterogeneity [27]. The Patch Cohesion Index (COHESION), Connectance Index (CONNECT), Integral Index of Connectivity (IIC), and Contagion Index (CONTAG) were used to measure landscape patterns. CONTAG, COHESION, Patch Density, and Edge Density were used in this analysis. We used fuzzy membership functions and fuzzy operators to integrate the four chosen fragmentation parameters. To build fuzzy crisp layers, we used the MS large membership function for all fragmentation parameters. Then, to create the final fragmented layers, we used fuzzy OR operators to combine all the parameters.
Ecosystem resilience (R) relates to an ecosystem's capacity to maintain its original function and structure increasingly under human pressure [28]. LULC contributes to ecosystem R in a variety of ways [29], and the following Equation (4) was used to calculate Abha's ecosystem R: The resilience was calculated in this study using the coefficients available for each land use type (Supplementary Table S1).
In accordance with Costanza's ESV model [24], the unit value of ES for specified regions has been widely designed to quantify ES. For 2000 and 2018, we used an ESV coefficient (Supplementary Table S2) to get the ES value. The following equation is used to calculate the ES value (5).
where A i denotes the area of LULC type i, n denotes the number of LULC types, and P i denotes the values per unit area of LULC type i.

Sensitivity Analysis of Ecosystem Health Model Using Morris Method
Morris (1991) developed the Morris one-at-a-time method (MOAT) for screening the parameters as a global sensitivity analysis method. The theoretical basis of this approach is the calculation of the overall effect and interaction effect of all input parameters using the mean m and standard deviation s of all gradients sampled from the r MOAT path. All input parameters have been modified by the same relative amount using this global sensitivity method. The Morris method considers modifying the variable in question between two model simulations, which separate it from the conventional OAT analysis approach [30].

LULC Mapping and Dynamics
An SVM classifier was used to build LULC maps of Abha, Saudi Arabia in 1990, 2000, and 2018 (Figure 2a-c). Built-up areas, water bodies, dense vegetation, sparse vegetation, agricultural land, scrubland, bare soil, and exposed rock were all classified as LULC types in the present study area. The areas under various LULC types were computed and are summarized in Table 1 for the term 1990-2018. In 1990, Abha city had a built-up area of 6246 ha, which had risen to 27,145 ha in 2018 (Table 1). Furthermore, the delta (∆) change in the built-up area between 1990 and 2000 was 4275 ha (3.34%), which rose to 16,624 ha (13%) (2000-2018) and 20,899 ha (16.34%) . From 1990 to 2018, the built-up area grew at a rate of 334.6%, indicating significant urban expansion over a 28-year period. Scrubland, one of the eight LULC types, experienced a significant decline over 28 years. Between 1990 and 2018, the area of scrubland declined by 14,528 ha. Figure 2 shows that the built-up area has converted the majority of the scrubland. Other land-use types, such as dense vegetation and exposed rocks, witnessed almost no change during this time span. The water body and sparse vegetation endured minor losses of approximately −0.07% and 0.77%, respectively (Table 1). Since 2000, the scrubland in the study area has seen a rapid conversion to other land uses, especially built-up areas. One of the key factors for this land-use type's higher conversion rate was its proximity to the main city as well as its topographic location, which has a lower elevation (2000 m) and lower slope (10 • ) than other land-use types.
2018, the built-up area grew at a rate of 334.6%, indicating significant urban expansion over a 28-year period. Scrubland, one of the eight LULC types, experienced a significant decline over 28 years. Between 1990 and 2018, the area of scrubland declined by 14,528 ha. Figure 2 shows that the built-up area has converted the majority of the scrubland. Other land-use types, such as dense vegetation and exposed rocks, witnessed almost no change during this time span. The water body and sparse vegetation endured minor losses of approximately −0.07% and 0.77%, respectively (Table 1). Since 2000, the scrubland in the study area has seen a rapid conversion to other land uses, especially built-up areas. One of the key factors for this land-use type's higher conversion rate was its proximity to the main city as well as its topographic location, which has a lower elevation (2000 m) and lower slope (10°) than other land-use types. A transitional probability matrix was used to compute the trend of LULC transformation from each type to other LULC types. It denotes the likelihood of each cell of land-use types being transformed into another land-use type. The transitional probability matrix was constructed in the current study using the MOLUSCE plugin of QGIS (version 2.18.33) and the Markovian approach for the periods 1990-2000, 2000-2018, and 1990-2018 (Supplementary Tables S3-S5). The transitional probability matrix between the 1990 and 2018 LULC maps is shown in Supplementary Table S5. The diagonal bold values indicated that the probabilities remained constant over the period. A value of zero indicates that the land is not being used. The built-up area in the study area was the most consistent and stable land-use type, with 83.6% of the largest transition probability. Sim-   Tables S3-S5). The transitional probability matrix between the 1990 and 2018 LULC maps is shown in Supplementary Table S5. The diagonal bold values indicated that the probabilities remained constant over the period. A value of zero indicates that the land is not being used. The built-up area in the study area was the most consistent and stable land-use type, with 83.6% of the largest transition probability. Similarly, exposed rocks are a consistent landuse type with a high transitional probability of staying unchanged by 75.08%. Meanwhile, agricultural land had the lowest transitional probability of remaining unchanged by 15.3% over a 28-year period. However, the built-up area was not drastically changed to other land uses. The amount of built-up area in scrubland and exposed has gained by 5.8% and 7% over 28 years. The study area has experienced rapid urbanization. Sparse vegetation, agricultural land, scrubland, bare soil, exposed land, and water bodies were converted into urban areas at a rate of 15.2%, 17.2%, 21.8%, 23.8%, 12.4%, and 19.5%, respectively. As a result, urban areas have absorbed all resource-producing land-use types in the name of development, which is extremely harmful for the environment and ecosystem.
The goal of accuracy evaluation is to quantify how well pixels were sampled into the right land cover groups. A Khat statistic (an estimate of kappa) is produced by kappa analysis and is a measure of agreement or accuracy between ground reality and classified LULC maps. A kappa coefficient of one indicates complete agreement, whereas a number close to 0 indicates agreement that is no better than would be anticipated by chance. The kappa and RMSE errors for all LULC images are presented in Table 2. The LULC map of 2018 had a kappa coefficient of 0.87, indicating that the ground truth and the LULC map were very good at each other. While there is very little RMSE for the 2018 LULC map. Thus, the performance of the LULC classifier is highly satisfactory in preparation of the 2018 LULC map. Furthermore, in the 1990 and 2000 LULC maps, the kappa coefficient was more than 0.8, which indicated good results. The overall classification accuracy for classified LULC maps from 1990, 2000, and 2018 is 88.64%, 87.60%, and 90.71%, respectively. The kappa coefficients for 1990, 2000, and 2018 classified LULC maps are 0.85, 0.84, and 0.87, respectively (Table 2). Furthermore, the measure of producer accuracy (sensitivity) shows the specific category's prediction accuracy. The user accuracy indicates the classification's dependability to the user. The accuracy of the user is a more important determinant of the classification's actual applicability in the field. In the present study, water bodies had the highest producer's accuracy (96.92%, 96.91%, and 98.01%), while sparse vegetation had the lowest producer accuracy (82.31%, 80.92%, and 84.01%) for the LULC of 1990, 2000, and 2018 respectively (Table 2). However, the accuracy above 80% is considered as good agreement between ground reality and classified images. On the other hand, water bodies had the highest user's accuracy or reliability (97.82%, 99.93%, and 98.71%), while built up and sparse vegetation had the lowest producer accuracy (81.88%, 75.1%, and 75.07%) for the LULC of 1990, 2000, and 2018 respectively ( Table 2). As a result of the statistical validation analysis, all LULC maps were found to be valid and can be used for further research.

LULC Prediction
The current study is designed to forecast the future LULC of 2028 based on historical LULC maps and conditioning factors. The first step in forecasting future LULC is to predict current or existing LULC maps. If the predicted model yields a satisfactory result, the same model can be used to forecast future LULC maps. Since it is impossible to evaluate the model's performance for future LULC, predicting an existing LULC map is a necessary task. Since we have no idea of future LULCs, therefore, the validation of future LULC is impossible. Hence, the model's performance has been evaluated by simulating the existing LULC through optimizing the model's parameters before it can be used to forecast future LULC maps. Otherwise, there would be no reliability of the future LULC map. Thus, the LULC map of 2018 was first simulated, and the model performance was assessed in the present study. The same model was applied for the 2028 LULC map after obtaining the good performance of the model. In the present study, seven LULC change conditioning factors were identified for 2018 and 2028, including elevation, slope, proximity to an urban area, agricultural land, scrubland, sparse vegetation, and water bodies (for detail of the role of factors for projecting future LULC, see Supplementary Materials, Supplementary  Figures S1 and S2). The ANN model was used to integrate the conditioning parameters for predicting the LULC transitional suitability map after obtaining change maps and transition probability between the LULC periods of 1990-2000 and 2000-2018 using QGIS. The QGIS MOLUSCE plug-in includes a system for gathering training and testing datasets. The back-propagation algorithm based multi-layer perceptron architecture was used to generate the LULC transitional suitability model. This step's input data consist of a change map (which was derived by change analysis between the initial layer (2000) and last year (2018) LULC) and are expected to influence the likelihood of changes occurring. The weights of the components selected by the model are the output values. As a result of tuning the model, each factor receives one or another weight, depending on its contribution to the probability of changes occurring. Weight can be either positive (a close relationship between the factor and the possibility of changes occurring) or negative (feedback-if the factor is present, then changes are unlikely). Then, the output of the ANN model was integrated with the CA model. Several model parameters have been optimized during the ANN model's training to achieve better results. We used a trial-and-error process to optimize the ANN model parameters. The optimized parameters are are follows: Iteration rate: 1000, Learning rate: 0.001, Momentum: 0.02, Neighborhood: 10 px, Hidden layer: 10, Activation function: numpy.tanh sigmoid function. The prediction of LULCs was carried out using CA simulation after obtaining the land-use probability or suitability model by ANN. Then, the CA model was optimized. We started by setting the model to one iteration (10 years), which represents the forecast for the next ten years.
The CA model was used to create the LULC maps for 2018 and 2028, which were based on the ANN model-based land suitability model. After generating the simulated LULC map of 2018, it was compared to the original LULC map of 2018 to judge or validate the ANN-CA model's prediction performance. For the LULC map of 2018, the correction of prediction, correlation coefficient, and kappa coefficient were 76.2%, 0.841, and 83%, respectively, indicating that the model performed well. Following the validation of the predicted LULC map for 2018, the LULC map for 2028 was created using the same optimized model that used the LULC change conditioning parameters for 2018 as well as the LULC maps for 2000 and 2018. During the 2028 modeling procedure, no model parameters have been changed. Figure 3 depicted the 2018 simulated LULC map and the 2028 predicted LULC map. Since the LULC of 2028 cannot be validated, a sensitivity analysis was conducted for the sake of the predicted results' reliability. rection of prediction, correlation coefficient, and kappa coefficient were 76.2%, 0.841, and 83%, respectively, indicating that the model performed well. Following the validation of the predicted LULC map for 2018, the LULC map for 2028 was created using the same optimized model that used the LULC change conditioning parameters for 2018 as well as the LULC maps for 2000 and 2018. During the 2028 modeling procedure, no model parameters have been changed. Figure 3 depicted the 2018 simulated LULC map and the 2028 predicted LULC map. Since the LULC of 2028 cannot be validated, a sensitivity analysis was conducted for the sake of the predicted results' reliability. The dynamics of LULC land use from 1990 to 2028 are presented in Table 3, revealing that the area being constructed was 34,372 ha, up by 28,134 ha (22.14%) and 7291 ha (5.74%) compared to the 1990 and 2018 LULC maps. The area of bare soil declined from 2018 to 1471 ha and from 1990 to 7246 ha. Since 2018, the exposed rock has fallen by 8648 ha and since 1990, it has fallen by 7922 ha. Similar to the 2018 LULC, other land-use types would remain in 2028. The built-up area in the study area was the most consistent and stable land-use type, with a 96.5% transitional probability of staying unchanged, but it gained a significant amount of area from other land-use types. Based on this analysis, it is possible to conclude that the urbanization process has accelerated in order to capture other land-use types in the process of urbanization and development. As a result, this period can be referred to as the period of urbanization in the study area. The dynamics of LULC land use from 1990 to 2028 are presented in Table 3, revealing that the area being constructed was 34,372 ha, up by 28,134 ha (22.14%) and 7291 ha (5.74%) compared to the 1990 and 2018 LULC maps. The area of bare soil declined from 2018 to 1471 ha and from 1990 to 7246 ha. Since 2018, the exposed rock has fallen by 8648 ha and since 1990, it has fallen by 7922 ha. Similar to the 2018 LULC, other land-use types would remain in 2028. The built-up area in the study area was the most consistent and stable land-use type, with a 96.5% transitional probability of staying unchanged, but it gained a significant amount of area from other land-use types. Based on this analysis, it is possible to conclude that the urbanization process has accelerated in order to capture other land-use types in the process of urbanization and development. As a result, this period can be referred to as the period of urbanization in the study area.

Sensitivity Analysis, and Correlation Analysis
In the present study, the sensitivity analysis was performed using several machine learning algorithms and statistical techniques, such as RF, CART, and probability distribution function, and for simulated LULC of 2018 and projected LULC of 2028 ( Figure 4). Pearson's correlation technique was used to explore the relationship of land use changing variables with simulated LULC of 2018 and predicted LULC of 2028. For the sensitivity analysis, the ANN-based land-use suitability or probability models (generated using MO-LUSCE plug-in of QGIS) for 2018 and 2028 were used. The land suitability models were used as the target variable, while seven LULC change conditioning variables were used as independent parameters. Before performing sensitivity analysis, we randomly collected data from target and independent variables based on 5000 points. Based on the extracted data, we performed a sensitivity analysis to explore the influence of the independent variables for modeling the target variable or land suitability models, based on which the LULC projection can be prepared.

Sensitivity Analysis, and Correlation Analysis
In the present study, the sensitivity analysis was performed using several machin learning algorithms and statistical techniques, such as RF, CART, and probability dis tribution function, and for simulated LULC of 2018 and projected LULC of 2028 (Figur 4). Pearson's correlation technique was used to explore the relationship of land us changing variables with simulated LULC of 2018 and predicted LULC of 2028. For th sensitivity analysis, the ANN-based land-use suitability or probability models (generated using MOLUSCE plug-in of QGIS) for 2018 and 2028 were used. The land suitabilit models were used as the target variable, while seven LULC change conditioning varia bles were used as independent parameters. Before performing sensitivity analysis, w randomly collected data from target and independent variables based on 5000 points Based on the extracted data, we performed a sensitivity analysis to explore the influenc of the independent variables for modeling the target variable or land suitability models based on which the LULC projection can be prepared. Multiple features are involved in the decision-making process utilizing RF and CART, and it is necessary to consider the relevance and implications of each feature, therefore assigning the relevant feature at the root node and traversing the node splitting downward. Moving in the downward direction reduces the degree of impurity and uncertainty, resulting in improved categorization or elite split at each node. Splitting metrics such as entropy, information gain, Gini index, and others are utilized to resolve the issue. However, the Gini index outperforms all others. As a result, the Gini index utilizing RF and CART was utilized to calculate the essential parameters in the current investigation. The Gini index, also known as Gini impurity, estimates the likelihood of a given characteristic being wrongly categorized when chosen at random. If all of the elements are connected to a single class, it is said to be pure. The Gini index ranges from 0 to 1, with 0 indicating purity of categorization and 1 indicating a random distribution of components across multiple classes. The Gini index value of 0.5 indicates an equal distribution of components across some classifications. When developing the CART and RF, the characteristics with the lowest Gini index value would be prioritized. However, the RF algorithm identified proximity to urban areas as the main indicator (Gini value: 5), which is responsible for LULC transformation based on the parameters of 2018 for the prediction LULC of 2028, followed by proximity to vegetation (Gini value: 15), proximity to agricultural land (Gini value: 12) (Figure 4a). Meanwhile, the CART model also predicted proximity to the urban area, proximity to vegetation, proximate to agricultural land, proximity to scrubland, and elevation as highly influencing parameters for LULC transformation (Figure 4b). For the sensitivity analysis of simulated LULC of 2018, we again applied RF and CART for exploring the influence of independent parameters for land-use suitability models. Figure 4c,d represented that proximate to urban and vegetation were identified as the most responsible indicators for a land suitability model of simulated LULC maps of 2018 based on the parameters of 2000 by RF and CART.
In the present study, we prepared PDFs and correlation coefficient for seven independent parameters and LULC probability or suitable models of 2000 and 2018 for simulating and predicting the LULC of 2018 and 2028 to explore the data pattern of the parameters and compare the shape of the PDFs to analyze the sensitivity (Supplementary Figures S3-S5). Based on the analysis of the data distribution pattern, it can be stated that proximity to urban, vegetation, and agricultural land can influence or control maximally the land suitable model.
We also computed the relation between the land suitability model and land-use changing drivers or variables using Pearson's correlation coefficient. Supplementary Figure S3 showed that proximity to the urban area had a higher correlation (r: 0.714) at the significance level of <0.01, while the elevation had a very low correlation with the land suitable model for LULC of simulated 2018. Figure 5 showed that proximity to the urban area had a higher correlation (r: 0.569) at the significance level of <0.01, while the elevation had a very low correlation with the land suitable model for LULC maps of 2028.

Changes in EH Indicators
The results of spatiotemporal analysis of vegetation showed that the vegetation ha significantly decreased from 1990 to 2028 (Supplementary Figures S6-S8). The vegetation has decreased by 148.4 km 2 for the period of 1990-2028. Therefore, it can be concluded that the declines in vegetation have led to a detrimental impact on the ecosystem health In addition, natural ecosystem services have been harmed.
On the other hand, we estimated the valuation of the ecosystem services using th PtW: proximity to water bodies, PtSL: proximity to scrub land, PtV: proximity to vegetation, PtU: proximity to urban areas, PtAL: proximity to agricultural land.

Changes in EH Indicators
The results of spatiotemporal analysis of vegetation showed that the vegetation has significantly decreased from 1990 to 2028 (Supplementary Figures S6-S8). The vegetation has decreased by 148.4 km 2 for the period of 1990-2028. Therefore, it can be concluded that the declines in vegetation have led to a detrimental impact on the ecosystem health. In addition, natural ecosystem services have been harmed.
On the other hand, we estimated the valuation of the ecosystem services using the coefficient of Costanza et al. (1997) from 1990 to 2028. The estimated ESV values for different LULC types are presented in Table 3. The results showed that the ESV for the urban area was zero (0) for the period of 1990-2028 (Table 3) because Costanza et al. (1997) reported that urban areas did not provide any ecosystem services. As the results showed that the areas of water bodies, sparse vegetation, and scrubland have been decreased from 1990 to 2028, therefore, the ESV from these LULC types would also experience the decreasing tendency accordingly. Table 3 showed that the ESV for water bodies were 1,155,728 US$ ha −1 yr −1 in 1990, while it was decreased to 650,624 US$ ha −1 yr −1 in 2018 as per the coefficient of 1997 of Costanza et al. (1997). On the other hand, ESV for scrubland was 1,978,960 US$ ha −1 yr −1 in 1997, whereas it decreased to 138,323,698 US$ ha −1 yr −1 in 2018. It was reduced significantly. However, the ESV from dense vegetation increased slightly. The overall ESV from the LULC was 14,591,612 US$ ha −1 yr −1 in 1990, while it was reduced by 3,617,680 US$ ha −1 yr −1 in 2018. On the other hand, the ESV value of LULC in the study area decreased to 3,194,714 US$ ha −1 yr −1 during the period 1990-2028 using the coefficient of Costanza et al. (1997). The huge reduction in ESV was caused by the very rapid urbanization process. Talukdar et al. [31] estimated the ecosystem services for different LULC maps in the charlands of the Ganges River and reported that due to urbanization, the overall ESV decreased significantly.
The resilience shows the ecosystem's ability to resist external pressures associated with LULC processing. Greater resilience increases the ecosystem's ability to resist or to cope with external pressure. Therefore, the results of this study showed that high resilience, dense vegetation, sparse vegetation, and scrubland for the period 1990 to 2028 were found in water bodies. During the study, a major process of urbanization was observed, capturing and converting the regions of natural resources to urban areas. As a result, the ecosystem's resilience has been impaired or reduced during this time. The health of the study area is affected by the ecosystem.
The preparation of fragmentation or organization (O) parameter is quite a difficult task, as there are no proper ways to derive it. Each fragmentation index represents the different organization of the study area. In the present study, we adopted the fuzzy logic model to derive the high-precision fragmentation parameter at a spatial scale. First, the fragmentation indices, such as cohesion, contag, patch density, and edge density were derived from the LULC maps of 1990, 2000, 2018, and 2028 at spatial scale using FRAGSTAT software. Then, we used fuzzy MS large membership function (high membership indicates high fragmentation and vice versa) to make the parameter standardized (value ranges from 0 to 1). Subsequently, we applied the fuzzy OR operator (to combine all parameters based on the high membership function) to integrate all parameters and construct the final fragmentation (O) parameter. In the present study, the results showed that fragmentation status has been increased over the periods (Supplementary Figures S10-S13). Whereas the natural resources, such as dense vegetation, sparse vegetation, and scrubland observed higher fragmentation value, which indicates that these natural resources have been undergone a huge transformation. Consequently, it hampers the ecosystem's health condition.

Ecosystem Health Conditions Models
In the present study, we integrated the knowledge of fuzzy logic with the VORS model to prepare the ecosystem health conditions model for the period of 1990-2028. Before proceeding with the VORS model, four parameters, such as V, O, R, and S were standardized by using the fuzzy membership function. We used the MS small function (low membership function of the parameters inversely affect the ecosystem health conditions) for V, R, and S parameters to make them fuzzy crisp and standardized. On the other hand, we used the MS large membership function (high membership highly influences the ecosystem health conditions) for fragmentation (O). Based on the standardized parameters, the VORS models were applied to integrate the parameters. Finally, the ecosystem health conditioning model for 1990, 2000, 2018, and 2028 was generated ( Figure 6). The ecosystem health condition maps were classified into five classes using interval techniques, such as very poor, poor, moderate, good, and very good ecos health conditions for the periods of 1990-2028 ( Figure 6). The area coverage of dif ecosystem health condition zones was computed and tabulated in Table 4 R showed that 234.63 km 2 and 450.61 km 2 areas are predicted as very-poor and poo system health conditions in 1990, while the areas of very-poor and poor health cond have been increased to 442.02 km 2 and 269 km 2 in 2028 ( Table 4). The very-good e tem health conditions have been increased slightly during the period (from 2.64 k 9.1 km 2 ). However, these findings do not mean that the environmental conditions been improved, because the huge area from moderate ecosystem health condition been shifted to very poor and poor ecosystem health condition zones. Meanwhi very slight area was shifted toward very good and good ecosystem health con zones. Therefore, these findings showed alarming environmental situations.  The ecosystem health condition maps were classified into five classes using equal interval techniques, such as very poor, poor, moderate, good, and very good ecosystem health conditions for the periods of 1990-2028 ( Figure 6). The area coverage of different ecosystem health condition zones was computed and tabulated in Table 4 Results showed that 234.63 km 2 and 450.61 km 2 areas are predicted as very-poor and poor ecosystem health conditions in 1990, while the areas of very-poor and poor health conditions have been increased to 442.02 km 2 and 269 km 2 in 2028 ( Table 4). The very-good ecosystem health conditions have been increased slightly during the period (from 2.64 km 2 to 9.1 km 2 ). However, these findings do not mean that the environmental conditions have been improved, because the huge area from moderate ecosystem health conditions has been shifted to very poor and poor ecosystem health condition zones. Meanwhile, the very slight area was shifted toward very good and good ecosystem health condition zones. Therefore, these findings showed alarming environmental situations.

Sensitivity Analysis
As the validation of ecosystem health models is not possible because of the absence of ground truth data, the reliability of the models can be judged and explored by using sensitivity analysis. The uncertainty of the input parameters and output was investigated using sensitivity analysis. In the present study, we employed the Morris method as the global sensitivity analysis to investigate the uncertainty and sensitivity of the models. We collected 1200 samples for performing the Morris model in R studio (package: sensitivity) for the period of 1990-2028. Based on the collected samples, the probability distribution function and correlation coefficient were performed for input parameters and output for the period of 1990-2028 (Figure 7; Supplementary Figures S18-S20). The results of the PDFs for 2028 showed that the data distribution pattern of the resilience is quite similar to the output or ecosystem health model, which is followed by ESV, NDVI, and fragmentation ( Figure 7). A high correlation was found between output and NDVI (r: 0.98), followed by resilience (r: 0.9), ESV, and fragmentation. However, a quite similar pattern of data distribution can be found in the case of the ecosystem health models of 1990, 2000, and 2018 ( Supplementary Figures S18-S20). Therefore, it can be stated that resilience, NDVI, ESV, and fragmentation affect the output model truly.

Sensitivity Analysis
As the validation of ecosystem health models is not possible because of the absence of ground truth data, the reliability of the models can be judged and explored by using sensitivity analysis. The uncertainty of the input parameters and output was investigated using sensitivity analysis. In the present study, we employed the Morris method as the global sensitivity analysis to investigate the uncertainty and sensitivity of the models. We collected 1200 samples for performing the Morris model in R studio (package: sensitivity) for the period of 1990-2028. Based on the collected samples, the probability distribution function and correlation coefficient were performed for input parameters and output for the period of 1990-2028 (Figure 7; Supplementary Figures S18-S20). The results of the PDFs for 2028 showed that the data distribution pattern of the resilience is quite similar to the output or ecosystem health model, which is followed by ESV, NDVI, and fragmentation (Figure 7). A high correlation was found between output and NDVI (r: 0.98), followed by resilience (r: 0.9), ESV, and fragmentation. However, a quite similar pattern of data distribution can be found in the case of the ecosystem health models of 1990, 2000, and 2018 ( Supplementary Figures S18-S20). Therefore, it can be stated that resilience, NDVI, ESV, and fragmentation affect the output model truly. The Morris results were evaluated by comparing the mean and standard deviation of the distribution function of each input parameter (Figure 8). The mean of the absolute value of the elementary effect generates the overall sensitivity of the ith input parameters on the ecosystem health model or output response. In this case, this is the accuracy of the ecosystem health model procedure. The higher the calculated value of μi*, the higher the sensitivity of the ecosystem health condition procedure to the ith input parameters. A standard deviation of elementary effect represents the likely interactions with other pa- The Morris results were evaluated by comparing the mean and standard deviation of the distribution function of each input parameter (Figure 8). The mean of the absolute value of the elementary effect generates the overall sensitivity of the ith input parameters on the ecosystem health model or output response. In this case, this is the accuracy of the ecosystem health model procedure. The higher the calculated value of µi*, the higher the sensitivity of the ecosystem health condition procedure to the ith input parameters. A standard deviation of elementary effect represents the likely interactions with other parameters and/or that the parameter has a nonlinear influence on the ecosystem health model or output [32]. rameters and output (Figure 8). The graphs of mean value and standard deviation showed that the most influential parameter of the ecosystem health model of 1990 is resilience, followed by NDVI, and fragmentation, while the least affected parameter is ESV (Figure 8a). On the other hand, resilience was the most influential parameter for the ecosystem health model of 2020, followed by ESV, fragmentation, and NDVI (Figure 8b). In the case of the ecosystem health model of 2018, fragmentation was computed as the most sensitive parameter, followed by ESV and NDVI, while the least sensitive parameter was resilience (Figure 8c). The ESV was found as the most sensitive parameter for the ecosystem model of the future, followed by fragmentation and resilience, while the least sensitive parameter was NDVI (Figure 8d).

Discussion
The EH condition of the study has reflected the impact of LULC changes on ecological health [33][34][35][36]. Evaluation of ecosystem health based on the variations in the LULC types of the ecosystem permits us to spatially understand the drastic variations in health conditions [33]. Therefore, LULC classification with robustness is the primary element to obtain better EH results. The SVM-RBF kernel-based model performs better in classifying Results showed the plot of the mean value and standard deviation of all input parameters and output (Figure 8). The graphs of mean value and standard deviation showed that the most influential parameter of the ecosystem health model of 1990 is resilience, followed by NDVI, and fragmentation, while the least affected parameter is ESV (Figure 8a). On the other hand, resilience was the most influential parameter for the ecosystem health model of 2020, followed by ESV, fragmentation, and NDVI (Figure 8b). In the case of the ecosystem health model of 2018, fragmentation was computed as the most sensitive parameter, followed by ESV and NDVI, while the least sensitive parameter was resilience (Figure 8c). The ESV was found as the most sensitive parameter for the ecosystem model of the future, followed by fragmentation and resilience, while the least sensitive parameter was NDVI (Figure 8d).

Discussion
The EH condition of the study has reflected the impact of LULC changes on ecological health [33][34][35][36]. Evaluation of ecosystem health based on the variations in the LULC types of the ecosystem permits us to spatially understand the drastic variations in health conditions [33]. Therefore, LULC classification with robustness is the primary element to obtain better EH results. The SVM-RBF kernel-based model performs better in classifying the LULC with high precision [37,38]. However, the climate change in recent decades has affected LULC severely in Abha city, similar to other cities in the world [39]. Furthermore, some factors, such as the expansion of industrial activities, combined with global climate change and water scarcity, have resulted in a deterioration of vegetation coverage. Thus, land infertility may happen due to socio-environmental drivers such as the temporal migration of humans [40]. This factor plays a crucial role in the reduction of the size of vegetation coverage and its conversion to built-up land in the study basin [41,42]. The built-up area has increased significantly during the study period due to rapid urbanization. Topography might have played a significant role in the expansion of Abha city's builtup area. This finding is in line with other studies where huge human interventions are the main reason for landscape fragmentation, such as widespread urban pressure and ecological disruption [43,44]. However, this work did not include rural planning and conservation strategies in the modeling procedure of LULC change, which could warrant further examination.
Several techniques, including EMM, PSR, and VORS, are widely used to assess ecosystem health (EH) conditions. Of these techniques, VOR has paid considerable attention to the EH structure studies [45]. The main advantage of this model is that it can evaluate health through the condition of the ecosystem and outside disruption, and it disregards the capability of the ecosystem to give proper services to humans [46,47]. However, in the case of EH, they must meet human demand while also ensuring sustainability. The VORS method couples the structure of the ecosystem and ecosystem services to human requirements [46], and it is easy to define in theory and develop in practice [48]. Hence, it is appropriate for the research on the LULC change in areas with severe constraints of anthropogenic inputs. The designated weights allocated using fuzzy logic to these parameters of the VORS method have an influential impact on the evaluation outcomes.
Evaluating EH includes human subjective decisions that may take into consideration human requirements of ecosystem service [48]. Using the subjective technique to estimate the particular weights to be allocated to the parameters focuses on the relative importance of the various parameters. In this study, we used fuzzy logic coupled with VORS based on objective judgment, which is a novel method. The VORS method has been proven to be scientifically robust and has been extensively used in the EH condition studies worldwide [33,46,49]. Thus, we depend on professional knowledge and practical experience, as well as the actual state of Abha city, which presents its robustness, to estimate its weight. However, in order to adopt this method in other cities, it must be adjusted suitably based on the actual states of the study region.
Data selection has a critical effect on the fuzzy logic with VORS model outcomes. Generally, the remote sensing dataset can be utilized to monitor and evaluate the spatial heterogeneity of EH at various time scales [50]. Hence, remote sensing data and GIS tools are adopted extensively in EH conditions [51]. Although RS large-scale data can evaluate the space-time pattern and evolution of EH in the urban city, large-scale datasets have several limitations in accurately reflecting the state of EH. Thus, the ground-based field dataset was adopted to verify the outcomes of large-scale datasets [52].
The areas with "very-poor to poor" EH is mainly located within the northern and eastern zones, with high ecosystem health, and it mostly remained stable from 1990 to 2018. The main reason is that those zones are close to Abha city and susceptible to tourist activities, which may reduce the EH [53]. Qian et al. [54] mentioned that low EH level zones in China were easily influenced by anthropogenic activities. Thus, intense urbanization in the area may be triggered by the low level of vegetation cover [55]. In contrast, the areas with "good" EH was distributed in the western zone. The reason is that the areas that received strict protection and vegetation coverage had a high level of EH conditions [56]. This finding is consistent with the findings of He et al. [47], which discovered that moisture and land utilization intensity were the root variables that estimated the EH [47], which is a common cause of high EH levels.
However, population growth and economic activity have increased the consumption of ES associated with human inputs, and the existing literature reveals that the effect of economic activity and population aggregation on the EH is presented in the expansion of the built-up area, which spoils the ecosystem's structure and function. That led to the ecosystem worsening [35,46]. Due to socio-economic transformation, a large number of infrastructure activities have been built in the study area, resulting in increased landscape fragmentation and decreased ecosystem organization [57]. In recent years, the intensity of human input into Abha's natural ecosystem has been increased steadily. Thus, it requires restoring vegetation growth and improving ecological carrying ability.
As a result of rapid economic and population growth, as well as the expansion of built-up areas, the likelihood of fragmentation between landscape matrices has increased, limiting sustainable development. Mo et al. [58] pointed out that expansion is the key cause of the declination of the ecosystem organization. Earlier literature has indicated that vegetation cover plays a pivotal role in improving environmental quality and the sustainable development of urban cities [47,59]. During the fast urbanization process, much ecological land (i.e., wetlands) has been changed into industrial land and residential areas, which decreases the resistance of the natural ecosystem to external interference, triggering some ecological problems [57]. On the similar note, Li et al. [36] reported that the disparity between the supply and demand of ecosystem services was the main reason for ecosystem degradation, which limited the survival of humans and the development of socio-economic activities. Thus, to meet the future demand for ecological services brought about by population growth, Abha city has made tremendous efforts.
The expansion of the built-up area in Abha city has increased the landscape's fragmentation and deteriorated the ecosystem's balance. Previous studies have indicated that the increase in landscape fragmentation will lead to many eco-environmental problems, such as urban heat islands (UHI), air pollution, and geo-environmental hazards [36,60]. The large-scale expansion of built-up areas has extremely distressed the ecosystem's balance, triggering a reduction in the resistance of the natural ecosystem to outside disturbances [57]. Since 1990, the built-up area has increased by 16.34%, while it was only 3.34% between 1990 and 2000 and 13% between 2000 and 2018, which has exceeded the carrying capacity of resources and the environment, resulting in ecosystem deterioration.
In this research, we established an integrated fuzzy logic-based VORS method for EH conditions based on the spatiotemporal variations in the landscape categories. However, this study has some shortcomings which cannot be ignored. For example, due to some inconsistencies, data from the study area could not be obtained spatially homogeneously during the field sampling procedure for confirming the EH model. In addition, we used the RS data of 1990, 2010, and 2018 to evaluate the EH of Abha city and applied the field dataset of 2018 to confirm the outcomes of EH in 2018. As the field sampling procedure was only done in 2018, we could not confirm the outcomes of EH in 1990 and 2010, which is another drawback of this study. Thus, further study should be conducted to verify that the sampling should be uniformly located and tourist activities should be taken into account within the framework for the evaluation. The findings presented here can help land planners understand the spatial and temporal patterns of ecosystem health changes and the factors that influence them. Finally, sincere efforts are needed to examine the factors that drive EH as well as their effects, which will provide a scientific basis for establishing appropriate protection strategies and formulations.
The use of moderate-resolution satellite images such as Landsat is a significant drawback of this research. High-resolution images, such as SPOT-5, KOMPSAT series, IKONOS, QUICKBIRD, etc., would be very beneficial in obtaining a very high precision result. The lack of socio-economic data, which often impacts land-use changes, may result in poor prediction. We utilized fewer explanatory factors to forecast the LULCs for 2018 and 2028. Several fragmentation parameters could be estimated instead of only five parameters for better explanation the organization. On the other hand, the net primary product might be more precise vigor parameters than NDVI. In the future, more advanced machine learning techniques such as deep learning could be used to replace the VORS model to get very high accuracy results.

Conclusions
This study provides in-depth analysis on the LULC's dynamics from 1990 to 2018, as well as a prediction of the future LULC. Furthermore, this study proposes the development of a fuzzy-VORS model for analyzing ecosystem health conditions from 1990 to 2018, including a future period (2028). The future LULC and EH models were evaluated using a global sensitivity analysis. The results revealed that the built-up area has increased by 16.34% since 1990, while it was only 3.34% between 1990 and 2000 and 13% between 2000 and 2018. The process of urbanization has gained acceleration after 2000. Furthermore, the LULC map of 2028 was predicted by using the ANN-CA model. The ANN-CA performed satisfactorily. The LULC's of 2028 showed that the built-up area increased to 343.72 km 2 . On the other hand, bare soil and exposed rocks were reduced significantly by 86.48 km 2 and 14.71 km 2 , respectively. Agricultural land showed a significant decrease in 2028. To the best of the authors' knowledge, this is the first study in which we used machine learning algorithms to perform sensitivity analysis, such as CART and RF, PDF, and correlation coefficient. The results of ecosystem health conditions revealed that from 1990 to 2018, very-poor and poor health conditions increased, while very-good and good conditions decreased. Global sensitivity analysis showed that ecosystem organization and resilience parameters were found to be the most sensitive parameters, while ESV and NDVI were the least sensitive parameters for all periods.
We developed the level of ecosystem health in Abha city, Saudi Arabia, and analyzed the characteristics of its space-time evolution, which will assist local authorities and environment scientists in formulating precise protection measures for various regions with varying levels of ecosystem health, as well as provide guiding significance for the protection of Abha's EH. In addition, we created an EHA framework, which may be used for other cities in the future.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/rs13132632/s1, Figure S1: LULC change conditioning factors, such as (a) elevation, (b) slope, (c) proximity to urban area, (d) proximity to agricultural land, (e) proximity to scrubland, (f) proximity to sparse vegetation, (g) proximity to water bodies for predicting 2018 LULC; Figure S2: LULC change conditioning factors, such as (a) elevation, (b) slope, (c) proximity to urban area, (d) proximity to agricultural land, (e) proximity to scrubland, (f) proximity to sparse vegetation, (g) proximity to water bodies for predicting 2028 LULC; Figure S3: Correlation between simulated LULC for 2018 and LULC change conditioning parameters; Figure Table S2: Coefficient for ESV estimation; Table S3: Transitional probability matrix  between 1990-2000 LULC maps; Table S4: Transitional probability matrix between 2000 and 2018  land use maps; Table S5: Transitional probability matrix between the LULC maps of 1990 and 2018; Table S6: Transitional probability matrix between the LULC maps of 1990 and 2028.

Conflicts of Interest:
The authors declare no conflict of interest. One-at-a-time method