Land subsidence hazard assessment based on novel hybrid approach: BWM, weighted overlay index (WOI), and support vector machine (SVM)

Land subsidence is a morphological phenomenon, which causes negative environmental and economic consequences for human societies. Therefore, identifying the areas prone to subsidence can be one of the necessary measures for reducing the potential risks. This study aims to evaluate the efficiency of the support vector machine (SVM) algorithm and weighted overlay index (WOI) models in zoning the rate of land subsidence hazard in Hashtgerd plain, Iran. First, the 19 criteria include groundwater depletion, groundwater extraction, aquifer thickness, alluvium thickness, aquifer recharge, well density, drainage density, groundwater depth, lithology, bedrock depth, average annual precipitation, average annual temperature, climate type, agricultural use, urban use, industrial use, distance from rivers and streams, distance from roads, distance from faults were considered. Then, the layers were weighed based on the Best–Worst Method (BWM). The results of BWM indicated that the factors of groundwater extraction (0.219), lithology (0.157), and groundwater depletion (0.079) have a greater effect on the potential for subsidence hazard. Moreover, the results of validation by performing ROC curve showed that the accuracy of RBF-SVM, LN-SVM, SIG-SVM, PL-SVM, and WOI were 95.7, 94.3, 94.9, 93.2, and 90%, respectively. Based on the ROC results, all of the models for preparing the subsidence hazard map in Hashtgerd plain exhibit excellent accuracy. Therefore, all of the models used here can predict the areas vulnerable to subsidence properly. In this study, the five land subsidence hazard maps were used as new input factors and integrated using fuzzy gamma-ensemble methods to make better hazard maps. The results of the ensemble model indicated that 19.3% of Hashtgerd plain is in the zone of high to very high sensitivity. The results of this study can help planners in managing and reducing the possible hazards of subsidence.


Introduction
Subsidence is one of the hazards, which human have faced in recent decades, especially in alluvial plains (Fabris et al. 2021;Shi et al. 2021). Subsidence is a type of surface deformation, which is associated with vertical deformation, sudden sinking, or gradual settling of the surface of the earth due to the removal or displacement of subsurface material (Shi et al. 2016;Galloway et al. 2018). This phenomenon is affected by human activities such as overexploitation of aquifers, changing the landuse, constructing and loading the engineering structures, drainage of organic soils, subsurface mining or pumping the fluids from the ground such as oil and gas (Lee and Park 2013;Rudiarto et al. 2018;Guzy and Malinowska 2020;Gido et al. 2020;Fabris et al. 2021), as well as natural geological phenomena including subsurface dissolution and karst collapse, deposit density, and movements tectonic faults (Modoni et al. 2013;Pradhan et al. 2014;Sopata et al. 2020).
Subsidence caused by natural factors often occurs gradually and over a long period of time. However, subsidence caused by human factors usually occurs suddenly and leads to severe environmental and economic consequences such as morphological irregularities (Meilianda et al. 2010), damage to man-made facilities (Ibrahim et al. 2018), disturbance in the pattern of hydrological flows (Machowski and Rzetala 2018;Meldebekova et al. 2020), and destruction of subsurface facilities (Sayyaf et al. 2014;Sahu et al. 2015) in urban and coastal areas within less time. Therefore, recognizing the effective factors in creating subsidence and determining the degree of risk of areas is one of the basic steps for preventing the damage to vital buildings and infrastructure, planning for sustainable urban development, awareness of the temporal and spatial distribution of land surface deformation, and reducing the risk rate (Fiaschi et al. 2018;Fabris et al. 2021). Therefore, various researchers around the world studied various aspects of this issue, most of which included monitoring the subsidence phenomenon (Aditiya et al. 2017;Liu et al. 2020;Orhan 2021), how it is formed, and the mechanism of its formation (Sowers 1976;Nie et al. 2013;Yu et al. 2018). Further, various researchers modeled and zoned the hazard of subsidence occurrence by using different semi-quantitative and quantitative methods such as artificial neural network (ANN) (Jung et al. 2005;Kim and Lee 2020), weight of evidence (WOE) (Oh and Lee 2010), support vector machine (SVM) (Zhi-xiang et al. 2009;Sui et al. 2020), evidential belief function (EBF) (Pradhan et al. 2014), adaptive neuro-fuzzy inference system (ANFIS) , analytical hierarchy process (AHP) (Ibrahim et al. 2018;Rezaei et al. 2020), random forest (RF) algorithm (Mohammady et al. 2019a) and multiple criteria decision making (MCDM) (Ghorbanzadeh et al. 2018;Arabameri et al. 2021a). Implementing these models, along with the use of geographical information system (GIS), are suitable tools for estimating the potential of land subsidence based on achieving high-resolution data, analyzing the factors affecting this phenomenon, and reducing operating costs (Pradhan et al. 2014;Tien Bui et al. 2018).
Due to the pervasive environmental risk of land subsidence in more than 300 plains of Iran, the need to study about this hazard is important (Mirzadeh et al. 2021). Hashtgerd plain is one of these vulnerable plains in Iran. The studies of Geological Survey and Mineral Exploration of Iran (2007), which were performed by using the interferometric synthetic aperture radar (InSAR) method, indicated that land subsidence occurred with a maximum rate of 16 and the average rate of 8.4 cm per year in a wide zone of this plain. This plain has been seriously faced with a drought crisis and decreased groundwater levels in recent years. It seems that the groundwater levels in this area will decrease significantly in the future due to climate change and continuing the existing conditions in the exploitation of water resources. In fact, changes in agricultural pattern, reduced precipitation, and occurrence of continuous droughts lead to the unplanned and unprincipled use of groundwater resources and declining groundwater levels in Hashtgerd plain, which causes the conditions for the occurrence and expansion of the land subsidence. Therefore, providing a suitable model and preparing a land subsidence hazard zonation map can be a significant help to local people and competent centers for protecting the human and natural resources and reducing the damage caused by land subsidence. In this regard, this study aims to identify and classify the sensitive areas, and determine the most important factors affecting the occurrence of subsidence based on the weighted overlay index (WOI) model and support vector machine (SVM) algorithms, and evaluate the accuracy of these models in Hashtgerd plain. The most significant innovation includes: • The appropriate and comprehensive combination of 5 criteria and 19 sub-criteria effective in assessing the subsidence in Hashtgerd plain based on the opinion of experts and review of extensive sources. • Making critical and basic maps such as alluvial thickness, aquifer thickness, aquifer recharge, and bedrock depth for the first time in Hashtgerd plain for future researches. • Combining the WOI method and the best-worst method (BWM) in order to create a linear optimization problem of criteria weights for the first time. • Comparing the efficiency of the WOI method and 4 kernels of the vector machine algorithm including linear, polynomial, sigmoid, and radial to prepare the subsidence map in Hashtgerd plain.

Study area
Hashtgerd plain is located in Alborz Province, Iran. The area of this plain is about 1168.8 Km 2 , which is located between the latitude 35° 47′ 45" and 36° 03′ 05″ north and longitude 50° 29′ 05″ and 50° 54′ 28″ east. The climate in Hashtgerd plain is considered as dry and cold based on Emberger and dry on Dumarten climate classifications. The aforementioned area experiences cold and wet winters and dry summers, which is among the characteristics of a dry and cold climate. The annual average temperature in the study area equals 14 °C and its absolute minimum and maximum temperature equal − 16.1 and 38 °C, respectively. In addition, the 33-year average precipitation in Hashtgerd plain about 373 mm. The climate of the area in the northern parts is semi-humid, which gradually changes to semi-arid toward the south of the plain due to decreasing altitude (Fig. 1). Geological studies of the area are necessary for knowing the subsidence phenomenon. Based on the hydraulic conductivity, the geological formations of this area are divided into four main groups, which are described as follows. The first group includes the formations with very low permeability. This group includes the formations consisting of ancient and very fine-grained sediments of shale, siltstone, and sandstone such as Barut, Zagun, Dorud, and Shemshak formations. The low permeability formations belong to the second group. The group includes the formations related to the late Paleozoic to early Paleogene, including Ruteh, Elika, and Dalichai formations. These formations are predominantly carbonate and do not have high hydraulic conductivity. Moderate permeability formations are known as the third group, which mainly includes Hezar Dareh alluvial formation and somewhat Neogene 1 3 sediments. Hezar Dareh Formation is in fact Quaternary alluvial terraces consisting of rubble and sand with weak silty-clayey cement. This formation reduces the lateral recharge of the aquifer due to its small pores and low permeability. Formations with high to very high permeability are considered as the fourth group. This group includes young alluvial sediments and river alluvium. In general, the main aquifer of the plain is composed of high permeability sediments, which have high hydraulic conductivity.
Further, the tectonic structure of the studied area has been affected by Alborz tectonic activities. Therefore, most folds and faults have an east-west trend. Tectonic activities in the area have caused more faults than folds. Important faults in this area include Abyek, south Taleghan, and Valian faults. Moreover, faults of the area played a major role in forming the south mountain ranges of Taleghan as a horst and Hashtgerd plain as a graben (Geological Survey and Mineral Exploration of Iran 2007). Its lithological properties are shown in Fig. 2 and Table 1.

Research framework
This study was performed in five general steps including identifying and preparing the effective environmental criteria, weighting environmental criteria based on the Best-Worst Method (BWM), implementing the Weighted Overlay Index (WOI) model and Support Vector Machine (SVM) algorithms, preparing the map of the Subsidence Potential Index (SPI), and validating the models by using the Relative Operating Characteristic (ROC) curve. Figure 3 shows the general steps of the research.

Identifying and preparing the effective environmental criteria
Based on the conditions of Hashtgerd plain, reviewing the literature and opinion of experts, the effective criteria in the land subsidence zoning were compiled in the two physical-chemical and economic-social dimensions, five groups of criteria including hydrological, geological, climate, land use, and buffer, and 19 sub-criteria including groundwater depletion, groundwater extraction, aquifer thickness, alluvium thickness, aquifer recharge, well density, drainage density, groundwater depth, lithology, bedrock depth, average annual precipitation, average annual temperature, climate type, agricultural use, urban use, industrial use, distance from rivers and streams, distance from roads, distance from faults. In this study, Inverse Distance Weighted (IDW) method was used to prepare five maps including temperature and precipitation maps, groundwater depletion, groundwater extraction, and alluvium thickness. Based on the IDW method, the degree of correlation and similarity between neighbors is proportional to the distance, which can be defined as a function

Fig. 3
Steps of land subsidence zoning process in Hashtgerd plain with the inverse of the distance of each point from the neighboring ones (Uyan and Cay 2010). It is possible to control the significance of known points upon the interpolated values, based on their distance from the output point with IDW. The weights for samples in IDW decrease with an increase in distance between the known samples and the estimated points. These weights are controlled by weighting powers, so that greater powers reduce the effect of farther estimated points and smaller powers distribute the weights more uniformly among the neighbors' points (Ouabo et al. 2020). The IDW method is considered as appropriate when the sample points exhibit enough distribution at the local scale levels (Li et al. 2014), so the IDW method was utilized according to the distribution of the samples. In order to check the efficiency of the IDW method, the data were removed at one point and the actual value was compared with the calculated one based on the root mean square error (RMSE) method, indicating that the above-mentioned method is regarded as appropriate for preparing the maps intended for this area.
In this study, a land subsidence inventory map was acquired from the Geological Survey and Mineral Exploration of Iran (2020). A total of 293 LS events have been recorded in the study area. Some of the field photographs in this study area are shown in Fig. 1. In order to create land subsidence susceptibility maps by statistical, probabilistic, and soft computing approaches, land subsidence training and testing datasets are required (Pradhan 2010;Mokhtari and Abedian 2019). In this study, the models were developed and validated by using training and testing dataset, respectively. For this purpose, the land subsidence locations were randomly divided into two parts: 70% of these land subsidence (205 land subsidence locations) were used to train the model, and the remaining 30% (88 land subsidence locations) were used to validate the performance of the model.
A description of the criteria is given in Table 2. Each criterion should be indicated as a map layer in the GIS-based database since a set of criteria was specified. The description of the criteria is as follows.
-Distance from rivers and streams Watercourse is one of the important factors in inducing instability in the region (Ghorbanzadeh et al. 2018;Ranjgar et al. 2021). Therefore, this map was extracted from the 1:50,000 topographic map of the studied area. Then, a map of the distance from rivers and streams was prepared by using the Distance tool in ArcGIS software, the continuous values of which are between 0 and 4826 m. This layer was classified into 9 classes with equal intervals. Score 1 was assigned to the first class, which represents the closest distance from rivers and streams, and score 9 to the ninth class. The other classes were assigned the appropriate score on a scale of 1 to 9.
-Distance from fault Faults are considered as one of the important criteria in land subsidence events. Vibration related to the fault formation can cause an inelastic compaction due to increased effective stress in the soil, which leads to subsidence (Galloway et al. 1999;Pradhan et al. 2014;Navarro-Hernández et al. 2020). In this study, the faults of the region were extracted from the geological map with a scale of 1:100,000. After that, the distance from the fault map was prepared by using the Distance tool in ArcGIS software. In the continuous map, the values were determined between 0 and 15,030 m. This layer was classified into 9 classes with equal intervals. Score 9 was assigned to the first class, which represents the closest Table 2 Criteria affecting on land subsidence  Ebrahimy et al. (2020) distance from faults, and score 1 to the ninth class. The other classes were assigned the appropriate score on a scale of 1 to 9.
-Groundwater depletion and groundwater extraction One of the causes of land subsidence is related to its natural reaction against groundwater overexploitation and consequently depletion of the groundwater level due to reduced hydraulic pressure and increased stress on sediments (Hu et al. 2009;Huang et al. 2012). Groundwater depletion is influenced by factors such as the reduction in rainfall, long-term droughts, especially in arid and semi-arid areas such as the study area, and geological conditions, as well as human factors such as indiscriminate groundwater extraction. In fact, groundwater extraction is among the critical human factors in reducing the quantity and quality of underground water and the occurrence of subsidence phenomenon. In this study, the data related to water level changes during 2005-2020 (15 years) obtained from Alborz Regional Water Company (2020) were applied to prepare a map of groundwater depletion in the study area. Then, the groundwater depletion was prepared by using IDW interpolation method in 5 classes, that class 1 represents the lowest and class 5 represents the highest groundwater depletion. Then, scores 1, 3, 5, 7, and 9 were assigned to classes 1 to 5 of the map, respectively. In addition, the groundwater extraction rate was prepared by using the discharge data of the exploitation wells. After that, a continuous map of groundwater extraction was prepared by using IDW interpolation method, the values of which were determined between 7 and 1756 cubic meters. This layer was classified into 9 classes with equal intervals. Score 1 was assigned to the first class, which represents the lowest groundwater extraction, and score 9 to the ninth class. The other classes were assigned the appropriate score on a scale of 1 to 9.

-Aquifer recharge
Recharge is one of the effective hydrological factors in subsidence. The higher rate of recharge caused the higher hydraulic pressure and greater distance between the grains, which reduces the effective stress leading to the lower possibility of subsidence (Nadiri et al. 2018;Sadeghfam et al. 2020). This study used the Piscopo method (Piscopo 2001) for preparing the aquifer discharge layer. In the Piscopo method, the slope, precipitation, and surface permeability of the soil (depth of 0.5-2 m of surface soil) layers were converted to the raster layer and the continuous map of aquifer recharge was prepared by weighting and overlaying the layers. The values of which were determined between 0.39 and 1. Score 9 was assigned to the first class, which represents the lowest aquifer recharge, and score 1 to the ninth class. The other classes were assigned the appropriate score on a scale of 1 to 9.
-Alluvium thickness Alluvium thickness includes the distance between the land surface and the bedrock. In general, subsidence mostly occurs in the layers with higher alluvium thickness since much water can be extracted by increasing the alluvium thickness of the area. To prepare this layer, the amount of alluvium thickness in each well of the region was calculated by using the data of the surface geophysics and well logging obtained from Regional Water Company of Alborz (2020). Then, these points were interpolated based on the IDW method and converted to a raster file. Therefore, a continuous layer of aquifer thickness was obtained. The values of which were determined between 60 and 329 m. This layer was classified into 9 classes with equal intervals. Score 1 was assigned to the first class, which represents the lowest alluvium thickness, and score 9 to the ninth class. The other classes were assigned the appropriate score on a scale of 1 to 9.
-Climatic parameters Climatic parameters, especially temperature and precipitation, are one of the effective factors in the subsidence phenomenon due to their effect on groundwater (Zheng et al. 2017;Tafreshi et al. 2019). Therefore, in this study, tabular data such as location and information of meteorological stations related to the average annual precipitation and temperature of the studied area were imported as the point into the ArcGIS software, which were converted to the annual isohyetal and isothermal maps by using IDW interpolation methods. Information about climatic stations is presented in Table 3. Values of isohyetal and isothermal maps were determined between 199 and 850 mm per year and 0-14 °C, respectively. These layers were classified into 9 classes with equal intervals. In the average annual precipitation, score 9 was assigned to the first class, which represents the lowest average annual precipitation, and a score 1 to the ninth class. In addition, in the average annual temperature, score 1 was assigned to the first class, which represents the lowest average annual temperature, and a score 9 to the ninth class. The other classes were assigned the appropriate score on a scale of 1 to 9. Moreover, subsidence usually occurs in arid climates where groundwater extraction is higher and aquifer recharge is poorer. The map of climate type in the studied area was received from Alborz Meteorological Office. Based on the map, the climate of the area includes six classes, the first and sixth of which indicate the driest and the wettest climate class, respectively. Scores 9, 8, 6, 4, 2, and 1 were assigned to classes 1 to 6 of the map.
-Land use Land is used for different purposes and the effect of each of these factors on subsidence can be varied. Agricultural and urban uses are one of the main factors in increasing groundwater extraction and occurring subsidence. Further, the weight of these industrial structures causes the compaction of soil layers and land subsidence Zhou et al. 2020). Therefore, the land use map at the scale of 1:100,000 was obtained from the Forests, Range, and Watershed Management Organization of Iran. The agricultural, residential, and industrial areas were extracted from it. In this study, industrial, urban, agricultural, and other land uses were assigned a score of 7, 7, 9, and 1, respectively.

-Bedrock depth
In an area where the bedrock is located at a shallow depth, drilling a well is not possible due to low thickness of alluvium. Thus, groundwater is transported to the thicker areas, which leads to an increase in the subsidence and horizontal displacement of layers in the area (WRI 2014). In this study, the layer of the bedrock depth was prepared from the global map of depth to bedrock (https:// data. isric. org/ geone twork/ srv/ api/ recor ds/ f36117ea-9be5-4afd-bb7d-7a3e77bf392a). The values of which are continuous and determined between 0 and 231 m. This layer was classified into 9 classes with equal intervals. Score 1 was assigned to the first class, which represents the lowest depth of the bedrock, and score 9 to the ninth class. The other classes were assigned the appropriate score on a scale of 1 to 9.
-Distance from roads Communication axes increase the land subsidence and sinkholes by causing the waves to the sides and the depth of the earth and collapsing the soil particles in places where moisture decreased, which is related to much passages (Hakim et al. 2020;Arabameri et al. 2021b). Road maps were extracted from the 1:50,000 digitized topographic map of the studied area. Then, the map of distance from roads was prepared. In the continuous map, the values were determined between 0 and 8424 m. This layer was classified into 9 classes with equal intervals. Score 9 was assigned to the first class, which represents the closest distance from roads, and score 1 to the ninth class. The other classes were assigned the appropriate score on a scale of 1 to 9.
-Well density Density of wells and qanats directly related to the overexploitation of groundwater and amount of subsidence (Jeanne et al. 2019). In this study, the location of wells as x and y was received from Regional Water Company of Alborz (2020). The density of wells was calculated based on the number of wells per square kilometer after entering the ArcGIS software. In the map, its values were determined between 0 and 223 wells per square kilometer. This layer was classified into 9 classes with equal intervals. Areas with high well density, which are more prone to subsidence, were ranked higher and areas with low well density were ranked lower. -Lithology Rocks show different resistances against the external forces due to differences in the type of constituent sediments, conditions, and period of formation (Conforti et al. 2014; Mokhtari and Abedian 2019). Therefore, a digital lithological map of the studied area with a scale of 1:100,000 was received from the Geological Survey and Mineral Exploration of Iran (2020), and the characteristics of the lithological units of Hashtgerd plain were extracted. Hashtgerd area has 19 lithological formation classes, and these formations were ranked on a scale of 1 to 9 based on the type of constituent materials and with the help of review sources. Fine-grained formations, which are more prone to subsidence, were ranked higher and coarse-grained formations were ranked lower. The scoring of lithological formations is presented in Table 4.
-Drainage density The amount of drainage density is directly related to the subsidence events (Hakim et al. 2020). Drainage density, which is the ratio of the length of the rivers and streams to the studied area, was calculated by using a digital height map and ArcHydro extension in ArcGIS software. In the continuous map, the values were determined between 0.04 and 2.14 km/Km 2 . This layer was classified into 9 classes with equal intervals. Areas with high drainage density, which are more prone to subsidence, were ranked higher and areas with low drainage density were ranked lower.

-Aquifer thickness
The alluvium thickness is the distance between the ground surface and the bedrock, and the aquifer thickness is the distance between the water table and the bedrock. Therefore, the aquifer thickness can be easily estimated by subtracting the water table depth from the alluvium thickness. In other words, the part of the alluvium thickness saturated with water can be considered as the aquifer thickness. In this study, due to the relatively accurate map of alluvium thickness and water table depth, this method was used. In the continuous map of aquifer thickness, the values were determined between 0 and 285 m. This layer was classified into 9 classes with equal intervals. Score 1 was assigned to the first class, which represents the lowest aquifer thickness, and score 9 to the ninth class. The other classes were assigned the appropriate score on a scale of 1 to 9.
-Groundwater depth There is a direct relationship between groundwater depth and subsidence. The contour map of piezometers was obtained from Regional Water Company of Alborz (2020) and its values were determined between 10 and 145 m. The continuous map of depth of groundwater level was prepared by interpolating the contour map of piezometers based on the IDW method. The values of which are continuous and determined between 10 and 145 m. This Table 4 The scoring of lithological formations in land subsidence events No. layer was classified into 9 classes with equal intervals. Score 1 was assigned to the first class, which represents the lowest depth of the groundwater level, and score 9 to the ninth class. The other classes were assigned the appropriate score on a scale of 1 to 9.

Weighted overlay index (WOI)
Identifying and assessing the area prone to subsidence is one of the most important issues in assessment and spatial development programs (Abidin et al. 2015). Multi-criteria evaluation is a spatial decision support tool to identify the best place or the best pixels based on their ranking by evaluating some main criteria. The steps of this method included defining the problem, setting goals, selecting the criteria, standardizing, determining the weights, using one of the methods of multi-criteria evaluation, and analyzing the results (Saaty 2008;Song and Chen 2018). The Weighted Overlay Index (WOI) is one of the most common methods in the multi-criteria spatial decision making, which is widely used in the process of assessing and determining the vulnerable zones (Hailegebriel et al. 2007;Vázquez-Quintero et al. 2020). In this method, the decision maker directly assigns weights to the criteria based on the relative importance of each criterion. In this research, the best-worst method (BWM) has been used to determine the importance of criteria and sub-criteria.

Best-worst method (BWM)
Based on the BWM, best and worst criteria are determined by decision maker and a pairwise comparison is made between each of the two criteria (best and worst) and other criteria (Rezaei 2015). Next, a max-min problem is formulated and solved to determine the weight of various criteria (Foroozesh et al. 2022). Moreover, a formula is considered for calculating the incompatibility rate to check the validity of the comparisons. The steps of the method are as follows (Rezaei 2015).
Step 1 A set of decision criteria is determined. The set of criteria is defined as C 1 C 2 … C n , which is needed for making a decision.
Step 2 The best (most important, most desirable) and worst (least important and least desirable) criteria are identified. In this step, the decision maker defines the best and worst criteria in general and no comparison is made.
Step 3 The preference of the best criterion over the other criteria is determined by using the numbers between 1 and 9. The preference vector of the best criterion over the other criteria is displayed as A B = a B1 .a B2 … .a Bn . In the mentioned vector, a Bj indicates the preference of the best criterion B over the criterion j . It is clear that a BB = 1.
Step 4 The preference of all the criteria over the worst criteria is determined by using the numbers 1-9. The preference vector of the other criteria over the worst criterion is displayed as A W = a W1 .a W2 … .a Wn T , where a Wj indicates the preference of the worst criterion w over the criterion j . Thus, we have a WW = 1.
Step 5 The optimal values of the weights are determined. The pairs Further, the problem can be rewritten as follows, where W B indicates the weight of the most important criterion, W w shows the weight of the least important criterion, W j is considered as the weight of the criterion j, and a ik × a kj = a ij , and i and j represent a perfectly consistent pairwise comparison matrix.
In these cases, instead of minimizing the maximum value in the set of a Bj × a jW = a BW and a BW ∈ {1.2. … .9} . The consistency ratio can be calculated by using the consistency criteria given in Table 5 and the presented formula.
The results are more consistent when the values of consistency ratio are closer to zero. In the next step, the sum of the product of multiplying the relative weight in the value of the criterion leads to obtaining a potential of the map for the desired goal through Eq. (5).
where SPI indicates the amount of subsidence potential, W i is considered as the weight of each sub-criteria, and X i shows standardized value of each sub-criteria.

Support vector machine (SVM) algorithm
SVM model is based on the statistical learning theory. Based on the theory, the learning machine error rate bound for unclassified data can be considered as a generalized error rate (Fig. 4). These bounds are considered as a function of the total of training error rates, which indicate the complexity rate of the classifiers. The rate of training error and complexity of classifier should be reduced for minimizing the generalized error rate, which can be performed by maximizing the margin separating. Therefore, this classifier has an appropriate performance since it does not depend on the dimension of the input data. SVM model, which has been used extensively in recent decades, and is based on the nonlinear transformation synchronously with a high-dimensional feature space (Vapnik 2013).
To provide the subsidence zoning map using SVM, a remote sensing software called ENVI 4.3 was used. The SVM classifier in ENVI 4.3 software provides four types of kernels, including linear, polynomial, radial basis function (RBF), and sigmoid. The mathematical representations of each kernel (linear, polynomial, radial base function, and sigmoid) are listed below, respectively (Tien Bui et al. 2012): (4) Consistency ratio = * Consistency criteria where γ is the term gamma in the kernel function for all types of kernels except linear kernels. d is the term polynomial degree in a kernel function for a polynomial kernel. r is the term bias in the kernel function for polynomial and sigmoid kernels. γ, d, and r are usercontrolled parameters, because their correct determination increases the accuracy of the SVM solution.

Ensemble of SVM algorithm maps using fuzzy operator
In the next step, in order to improve the prediction accuracy based on multiple outputs of vector machine functions, a set of kernals were combined using a fuzzy operator. The fuzzy set theory was used to integrate the calculated subsidence hazard indices through each model in this study. The concept of fuzzy logic is to consider the spatial objects on a map as members of a set. In the fuzzy set theory, membership can take on any value between 0 and 1 reflecting the degree of certainty of membership (Park et al. 2014). Bonham-Carter (1994) introduced five operators that were useful for combining exploratory datasets: five operators, namely the fuzzy and, fuzzy or, fuzzy algebraic product,

Fig. 4 SVM principles.
A n-dimensional hyperplane differentiating the two classes by the maximum gap: B non-separable case and the slack variables ξ (Yao et al. 2008) fuzzy algebraic sum, and fuzzy gamma operator. This study applied the fuzzy gamma operator for combining the fuzzy membership functions. The Fuzzy Algebraic Product is defined as: where µ i is the fuzzy membership function for the ith map, and i = 1, 2,..., n maps are to be combined. The fuzzy algebraic sum is complementary to the fuzzy algebraic product, is defined as: The gamma operation is defined in terms of the fuzzy algebraic product and the fuzzy algebraic sum by: where λ is a parameter chosen in the range (0,1), and the fuzzy algebraic sum and fuzzy algebraic product are calculated using Eqs. (10) and (11), respectively. In the fuzzy gamma operation, when λ is 1 the combination is the same as the fuzzy algebraic sum, and when λ is 0 the combination equals the fuzzy algebraic product. The judicious choice of λ produces output values that ensure a flexible compromise between the 'increase' tendencies of the fuzzy algebraic sum and the 'decrease' effects of the fuzzy algebraic product (Park et al. 2014).

Evaluating the models
To evaluate the models, the Relative Operating Characteristic (ROC) curve was used for assessing the predictability of the models. The ROC curve is a graphical plot of the balance between the negative and positive error rates for each possible amount of cut-off. The area under the curve indicates the rate of system prediction by describing its ability to accurately estimate what occurred and what did not occur. The most ideal model has the greatest Area Under the Curve (AUC). The AUC values range from 0.5 to 1. The classification of the area under the curve includes excellent (0.9-1), highly appropriate (0.8-0.9), appropriate (0.7-0.8), average (0.6-0.7), and weak (0.5-0.6). The best accuracy of the prepared zonation maps is obtained when the amount of area under the curve is closer to 1 (Agarwal and Garg 2016;Chen et al. 2018).

Results
In this study, two physico-chemical and socio-economic dimensions including five criteria and 19 sub-criteria were compiled for determining the zones prone to subsidence. Figure 5 shows the criteria map for determining the zones prone to subsidence.

Results of a Weighted overlay index (WOI)
The WOI method was used for assessing the potential of subsidence in the area. In the used method, after preparing the 19 sub-criteria as the normalized raster layer, the next step related to determining the relative importance of each effective component and assigning the appropriate weight to each of the sub-criteria was performed. The weighting of each criterion was obtained based on the opinion of each expert separately by solving a nonlinear model in Lingo software. Then, the significance coefficient or weight of each criterion was obtained by averaging the opinion of all the experts (Table 6). * indicates the consistency of a comparison, the value of which for all the experts in this study is close to zero. Therefore, the results were very well consistent and acceptable.
At this stage, the layers were combined using the weight overlap index method to prepare a potential zones of subsidence. High values of the map indicate the high sensitivity to the subsidence occurrence and increasing the probability of subsidence occurrence.

Results of Support Vector Machine (SVM) algorithm
Determining the appropriate values for each of the four kernels is necessary for achieving the maximum efficiency in the SVM algorithm process. The setting parameter (C) should be determined for LN-SVM. However, C and the gamma parameter (γ) should be determined for RBF-SVM. Three parameters including C, γ, and bias (r) should be determined in the case of SIG-SVM. Four parameters C, γ, r, and degree of kernel (d) are needed for PL-SVM. Table 7 indicates the four parameters required for optimizing the algorithm for zoning the subsidence. Four output maps of the sensitivity of the basin to subsidence were extracted based on this algorithm (Fig. 7). The potential of subsidence was calculated between 0 and 1. The greater sensitivity to subsidence occurrence and the higher possibility of the occurrence were observed when its value is closer to 1.

Result of land subsidence vulnerability zoning
Finally, the subsidence potential resulting from the WOI and four kernal SVM algorithm was classified into five classes including very low, low, moderate, high, and very high based on the Jenks natural breaks classification method (Figs. 6 and 7). Figure 8 showed the percentage of land subsidence vulnerability in Hashtgerd plain based on the WOI and four kernal SVM algorithm. As shown in Fig. 8 16, 34.1, 18,87, 7.94, and 5.93%, respectively. In this study, a slight difference was observed in susceptibility classes, which could probably be related to the four functions, their modeling process, and parameters (C, γ, r, and d)) required to optimize the SVM algorithm in subsidence zoning. In addition, Fig. 9 shows the distribution of the five susceptibility classes of the ensemble of four SVMs-WOI models. The areas of very low to very high land subsidence susceptibility areas in the ensemble of four 31.9,10.5,9.7,and 9.5%, respectively. Areas with very high subsidence vulnerability zones are found in the center parts of Hashtgerd plain. Comparatively, eastern and northern parts have low and very low vulnerability areas.

Results of validation based on the ROC curve
The percentage of area under the ROC curve was used for validating the accuracy of the WOI method and the functions of the SVM algorithm (Fig. 10). The results of validating the performance of the mentioned models by implementing the ROC curve indicated that the accuracy of WOI, RBF-SVM, LN-SVM, SIG-SVM, and PL-SVM models were equal to 90, 95.7, 94.3, 94.9, and 93.2%, respectively. Based on the ROC results, all of the models for preparing the subsidence sensitivity map in Hashtgerd plain exhibit excellent accuracy. Therefore, all of the models used here can predict the areas vulnerable to subsidence properly. Also, the ROC curve for the ensemble model was 96%. Based on the ensemble of four SVMs-WOI models, which has the highest ability to predict the vulnerable areas to subsidence occurrence, 19.3% of Hashtgerd plain is in the zone of high to very high sensitivity. The output map of the SVMs-WOI models indicated that most zone of the Hashtgerd plain basin is very low to low sensitivity to subsidence occurrence. Based on the fuzzy gamma-ensemble model, the central part of Hashtgerd plain is more sensitive to subsidence and the possibility of this hazard occurrence is very high in this part of the basin. Moreover, these conditions are prevailing in the northwest of the plain, and the possibility of subsidence is high there. The overlap of this map with the criterion maps showed that most subsidence occurs in the areas near the faults, with the highest groundwater extraction, as well as in agricultural lands, old alluvial fans, and old alluvial terraces. Further, based on the BWM method and the opinion of the experts, factors including groundwater extraction (0.219), lithology (0.157), groundwater depletion (0.079), and groundwater depth (0.078) have greater importance on the potential of the subsidence hazard in Hashtgerd plain. Abdollahi et al. (2019) reported that the piezometric data (groundwater), NDVI, and altitude were the most effective factors in the occurrence Fig. 6 Zoning of subsidence vulnerability based on WOI method 1 3 of land subsidence in Kerman province, Iran. In addition, Motagh et al. (2008) demonstrated that extraction from the aquifers is one of the most important factors for the subsidence occurrence observed in recent years in the plains of Iran. Mohammady et al. (2019b) indicated that groundwater extraction, geological structure, distance from fault, and slope percentage played an effective role in subsidence events. Based on the results of weighting in the researches, selecting the criteria and their weighting was a function of environmental conditions and the location of the studied area. In fact, knowing the studied area was not needed if such conditions did not exist, in which case the results of the research would be completely wrong.
Based on the ROC model, the efficiency of SVM classification with different kernel functions was compared with the WOI classification method. Based on the results, the SVM method was regarded as superior to the WOI one in terms of accuracy. The high accuracy of SVM should be attributed to its ability to identify the optimal separating hyperplane. Statistically, the optimal separating hyperplane identified by the SVM algorithm should be generalized with the least error compared to other ones by other classifications. The SVM method is considered as superior to other existing techniques due to its valuable advantages. For example, the SVM does not face the problem of local optima in its training, builds the classifier with maximum generalization probability, determines its structure and topology optimally, and forms nonlinear functions easily with low calculations utilizing nonlinear kernels and the concept of inner product in Hilbert space. In addition, SVM searches the optimal hyperplane between classes and requires less training data for classification. Further, SVM applies the kernel method to transfer the data to a high-dimensional feature space where the data can be separated linearly when such a process cannot be implemented. Furthermore, SVM uses optimization rules to locate the optimal boundaries between the classes and is regarded as an appropriate alternative to other  Fig. 9 Zoning of subsidence vulnerability based on the ensemble of four SVMs-WOI models common classification algorithms. The ROC rate of the four support vector method kernels RBF-SVM, LN-SVM, SIG-SVM, and PL-SVM was 90, 95.7, 94.3, and 94.9, respectively, indicating that all of the kernels in the subsidence classification with the SVM method were ranked high, meaning its high ability to determine areas vulnerable to subsidence. However, the SVM algorithm has several key parameters which should be set correctly to achieve the best classification results, the exact determination of which may be a timeconsuming process requiring the user's skill and knowledge of the area. The WOI method is considered as appropriate for subsidence zoning due to the ROC rate (90%) and its acceptable accuracy although the results indicate that using the support vector machine method provides higher accuracy compared to the WOI method. In WOI method, the significance of criteria relative to each other can be considered based on the opinion of experts and regional conditions. Assigning weight to each criterion in the overlay process allows the researcher to control the influence of different criteria in the model. The BWM method was utilized to weight the criteria. Selecting the best and worst criteria, reducing the number of comparisons, and increasing the efficiency of the model are among the most significant features of the BWM method compared to other ones, despite its inability to evaluate subjectively and lack of clarity related to mapping the decision makers' perceptions with exact numbers, which requires the use of fuzzy set theory. In addition, the BWM method becomes complicated when there are a large number of criteria to compare,  7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49 Fig. 10 The model performance evaluation based on ROC curves ranged 0.857-0.894 indicating the SVM algorithms had a high validity for land subsidence zoning. Moreover, Tien Bui et al. (2012) zoned the sensitivity of landslide by comparing four models including regression, SVM algorithm with two radial and polynomial functions, Decision Tree (DT), and Bayesian theory. Validation results showed that RBF-SVM and PL-SVM functions had the highest prediction capability. Mohammady et al. (2019b) used the SVM algorithm and the WOE model based on Bayesian theory for assessing the sensitivity of land subsidence and showed the SVM model had a higher validity. Tzampoglou and Loupasakis (2017) examined land subsidence vulnerability and hazard mapping in Amintaio basin, Greece, applying multi-criteria assessment and indicated that the high correlation between produced maps and field mapping data proves the high value of maps and applied techniques in managing and reducing phenomena. In addition, Moharrami et al. (2020) investigated national-scale landslide susceptibility mapping in Austria using fuzzy best-worst multi-criteria (FBWM) decision making and fuzzy analytic hierarchy process (FAHP) and argued that the AUC of the PR for the FAHP and FBWM models equaled to 0.84 and 0.89, respectively. In fact, FBWM is regarded as a vector-based technique compared to the matrix-based FAHP method, which requires fewer pairwise comparisons, resulting in reducing the error and uncertainty in zoning. Further, the RBF function had the highest prediction among the functions of the SVM algorithm. RBF function is usually the most widely used function among other SVM functions and has the highest rate of subsidence prediction in the area due to good learning capability, fewer required parameters, and optimization in nonlinear conditions, which was confirmed by the results of the study of some studies (Xu et al. 2012;Bhavsar and Ganatra 2014;Mhetre and Bapat 2015;Chen et al. 2016;Abdollahi et al. 2019). In sum, the comparison between the studies indicates that each function and method is selected based on the appropriate application and its parameters, accuracy and time required, need and skill of the user, the scale of the studied area, available resources, and the like. In addition, each model can exhibit different results based on such cases and conditions of each region, and selecting the appropriate model requires more research based on different models.

Conclusion
Providing an appropriate model for identifying the areas prone to subsidence can be considered as a suitable solution for protecting human and natural resources and reducing the economic. Therefore, this study aims to model the probability of land subsidence occurrence in Hashtgerd plain based on the Weighted Overlay Index (WOI) model and Support Vector Machine (SVM) algorithm functions. In this study, 19 environmental criteria were used for modeling the hazard of subsidence. Based on the ensemble of four SVMs-WOI models, which has the highest ability to predict the vulnerable areas to subsidence occurrence, 19.3% of Hashtgerd plain is in the zone of high to very high sensitivity. In this plain, extracting from the groundwater resources more than the renewal capacity, construction of illegal agricultural wells without exploitation licenses causes the insufficient recharge of the aquifer and lack of compensation of the excess of extracting the groundwater. Therefore, various environmental, economic, and social problems arose providing the conditions for increasing the occurrence of land subsidence. In this regard, a correct understanding of the rate of groundwater renewability, aquifer recharge, and creating a balance between exploitation and recharge of the groundwater can be considered as one of the appropriate management strategies for reducing and preventing the subsidence phenomenon in the aquifer range of Hashtgerd plain. The authors believe that in future research, the use of other machine learning methods such as Integrated Neural Networks (INN), Back Propagation Neural Network (BPNN), and regression methods including Generalized Additive Model (GAM) and Boosted Regression Trees (BRT) and comparing its results with those of the present research can play an effective role in increasing the accuracy and efficiency of results. Further, comprehensive monitoring and evaluating the subsidence situation of the region by using time series analysis of radar images can help greatly in the land management plans and watershed management plans.
Funding This research did not receive any specific grant from funding agencies in public, commercial, or not-for-profit sectors.
Data Availability Not applicable.
Code availability Not applicable.

Conflicts of interest
The authors have no conflict of interest to declare.
Ethical approval Not applicable.

Consent to participate Not applicable.
Consent for publication Not applicable.