Prediction And Potential Detection of Land Subsidence By Integrating AHP And Programming Methods: Case Study of South Eastern Iran

One of the biggest risks that need to be considered before any construction operation such as road construction or construction is to determine the exact location of subsidence and predict the places with the potential to cause these complications.There are various methods for identifying and predicting the occurrence of these phenomena. In this article, hierarchical analysis methods (AHP and FAHP) and programming in MATLAB environment with the help of Evals function to study and map these hazards in the southwest of Iran has acted as a case study.Studies of lithological layers, slope, distance from fault and distance from waterway in both methods express very good and signicant results in predicting this type of hazards.According to these results, the most hazards occur in places that have evaporitic lithology and have the shortest distance from waterways and faults with a slope between 20 to 70 degrees.The output of these studies is to provide a map of subsidence (sinkhole) hazards for the rst time by integrating the stated methods with operational accuracy above 90% and introducing a very useful software in MATLAB program for experts. A programming method is introduced in MATLAB environment to nd subsidence points. Offers the integration of hierarchical methods and programming a new method in the study of subsidence. A model for subsidence hazard mapping was introduced.


Introduction
Subsidence refers to the phenomenon of landslide caused by various factors such as uncontrolled abstraction of groundwater, dissolution of various layers of gypsum or carbonate (Banks et al, 2020 [9], Brick and Alexander, 2021) [12]. In fact, excessive consumption of groundwater and excessive digging of wells, permitted and unauthorized water, reduces and decreases the level of groundwater resources, followed by subsidence of plains and sinkholes, Bhattacharya, 2021 [10].Subsidence has always been considered as a risk factor and recognizing areas with the potential to cause this natural phenomenon prevents the occurrence of many hazards,Mahm et al, 2017 [24], Gauci et al, 2019[18].sinkhole is formed when parts of the earth that have an evaporative, carbonate, or saline composition are dissolved by groundwater circulation, Waltham et al., 2016 [41]; Biswas et al., 2020 [11].These sinkholes can be several meters to several hundred meters deep, Uri et al, 202139]. This phenomenon can be observed in various forms such as bowl, horseshoe, and walls (Soumyajit, 2019).Khuzestan province in southwestern Iran as one of the agricultural hubs of Iran, due to uncontrolled abstraction of groundwater and the presence of evaporative formations as well as Gachsaran Formation and carbonate formations such as Asmari Formation as one of the high-risk areas for sinkholes Is considered, Pourmorad and Jahan, 2020[28].In this study, the aim is to study and identify areas with potential for subsidence in the construction parts of this province with the help of remote sensing studies using AHP and programming methods.Since the study of subsidence for the parts of the plain with the altitudes is different in terms of the existence of factors and layers of studies, in these studies, due to the higher importance of highlands, has focused on studying parts of the province that at different altitudes and in the vicinity of different formations, such as evaporites.

Geological Location Of The Study Area
The study area is located in southwestern Iran and in Khuzestan province (Fig. 1).The study area is part of The Zagros fold -thrust belt (Zagros FTB) in southwestern Iran. The Zagros fold-thrust belt in Iran forms the external part of the Zagros active orogenic wedge, Alavi, 2004 [6].
The Zagros Basin is characterized by a thick sequence of 7 to 14 km of sedimentary sediments in areas of extreme length and width along the northnortheastern edge of the Arabian plate, Bagheri et al, 2020 [8]. The formations that exist in the study area can be divided into 4 categories of clastic formations (including Bakhtiari and Aghajari formations), evaporitic (Gachsaran formation), carbonate (for example Asmari formation) and marlestone (Mishan formation).

Study Method
Studies have been conducted to determine areas with subsidence potential in Khuzestan province in two different ways to achieve high accuracy for applied studies, each of which is examined separately.
3.1. AHP method: In general, the steps of this method can be expressed in two steps: A-Collecting the required data and preparing information layers: At this stage, while reviewing previous studies, the required data and information were collected from various sources and effective information layers in subsidence, including the geological map of Khuzestan province, digital elevation model map of Khuzestan province (DEM), information Hydrometric stations and waterways of Khuzestan province and the role of faults in Khuzestan province were collected with the help of the Geological Survey of Iran.

B-Analysis and interpretation of data
At this stage, the collected information is analyzed using GIS and the decision tree is plotted as criteria and sub-criteria (Fig. 2). The effective values for each criterion were speci ed in the sinkhole potential and standardized by fuzzy membership and rasterization functions.
Then, using two methods of Analytic Hierarchy Process (AHP) and Fuzzy Analytic Hierarchy Process (FAHP), which are multi-criteria decision methods, the importance of the criteria is estimated and by combining all of them with the index overlap method, the potential zoning map Subsidence has been prepared.

Fuzzy inference system programming
These studies, similar to AHP studies conducted in the GIS environment, included the three main sections of library studies (gathering various information, including articles and related roles), operational studies in the MATLAB environment, and nally eld validation studies. In this study model, the aim is to determine the amount of subsidence susceptibility in different regions of Khuzestan province with the help of lithology layers and distance from waterway and distance from fault and slope of layers so that it can be updated.

Discussion
In this section, the two methods are discussed separately.

AHP method:
This method in 6 stages of studies includes determining the effective factor, selecting criteria, Accreditation, data processing, fuzzy and overlap.

Selection of criteria
In selecting the effective criteria in evaluating the general rule, it is assumed that these criteria should be determined in relation to the problem situation,Munier and Honoria, .
[26] 2021 The process of generating benchmark maps is based on GIS functions, which include entering geographic data (acquisition, reformatting, land referencing, collecting and documenting relevant data), data storage (non-spatial data), It is the processing and analysis of data (to obtain information) and the output of data,Carter et al, 2021 [13]. In order to compare different measurement scales (based on different indicators), scaling or standardization should be used, which is measured through the elements of dimensionless converted indicators,Fullér et al, 2019 [16].Based on the type of information available to create maps, benchmark maps can be classi ed into de nite, probabilistic, and fuzzy states,Hsu et al, 2021 [20]. In these studies, the fuzzy method has been used.

Accreditation
The importance of benchmark maps in achieving output is not the same. Therefore, it is necessary to score benchmark maps, or in other words, to be validated,Wang, 2016 [42]. The validity of each criterion indicates its importance and value compared to other criteria in location or routing operations,Tang and Hsu, 2018 [38].
Accreditation is based on expert knowledge and based on the opinion of experts, taking into account various factors such as study area, location or routing parameters, the impact of each parameter, etc. Shen et al, 2019 [35]. Accreditation is done in several different ways, in which the Analytic Hierarchy Process (AHP) model is used.
In fact, the process of hierarchical analysis is one of the decision-making systems for multiple criteria that is based on expert knowledge, Yamagishi andBhandary, 2017 [44]. In hierarchical analysis, it is possible to formulate the problem and consider different quantitative and qualitative criteria,Abdel Rahman et al, 2021 [1].After forming the pairwise comparison matrix and completing it with the mentioned methods, the weight of each criterion is determined. To do this, use the program written in the MATLAB software environment and Expert Choice software, and by entering the data of the pairwise comparison table, which are fuzzy, the validity of each criterion is determined (Tables 1 and 2).

Slope layer processing
Basin slope is one of the most important factors that control the surface ow time and water concentration in the river,Ali et al, 2019 [7].The slope of a basin affects the amount of in ltration, surface ow rate, soil moisture and the amount of groundwater in the basin,Nayak,2021 [27]. Basically, steeply sloping basins often have ood ows Abuzaid and Fadl, 2018 [4].In areas with high slopes, the velocity of water ow increases and as a result, with the dissolution of surface layers, especially in evaporative layers, the probability of subsidence increases Abuzaid et al, 2020 [3].
In mountainous areas, increasing the current intensity along with permeability can be considered as an effective factor in creating subsidence Abdel Rahman et al, 2018. [2] In these studies, after preparing the digital elevation model (DEM), the slope map of Khuzestan province was prepared using the Slope function (Figs. 3 and 4).

Waterway compaction layer processing
The network of basin waterways shows how runoff is discharged from the basin surface,Saidi et al, 2021 [32]. Therefore, the more the waterways of a basin evolve and develop, the more subsidence of that basin indicates,Meng et al, 2017 [25].
To prepare this layer, rst a digital le (Shap le) of waterways in Khuzestan province was prepared and combined with the waterways from the DEM layer in ArcGIS program.
The output layer in Google Earth software was then validated with real streams, and new streams were added and some removed (Fig. 5). After calculating the length of the nal lines, the density layer of the streams was calculated using the Density function in the ArcGIS environment (Fig. 6).

Lithology layer processing
The most important factor in causing subsidence is the resistance of the constituents of the basin surface to erosion,Lange et al, 2019 [22].The strength of formations depends on the type, type and composition of the constituent minerals, the way of crystallization, the amount of pores, strati cation, the way of erosion, interlayer, etc. Semeraro et al, 2016 [33].To prepare this layer and classify the different lithologies of the basin in terms of subsidence potential, the geological map of Khuzestan province prepared by the Geological Survey of Iran (center of Ahvaz) has been used (Fig. 7).In this layer, due to the different importance of rock layers in the rate of subsidence, rst, according to the type of formations, 4 types of lithology with evaporitic names (gypsum and anhydrite), detrital (conglomerate, sandstone, silt and clay), carbonate (dolomite) And lime) and marl were divided.In the next step, with the help of the Reclassify function, the maximum validity and weight are given to the evaporating layer, then to the carbonate layers, and nally to the marls and detritus.

Fault density layer processing
Another factor in creating subsidence is the distance of layers from faults in the region and the impact of these faults on them Gasperini et al, 2021. [17]By applying stress to different rock units, the tectonic agent causes faults and fractures in them, and these fractures create secondary porosity in the rock, which increases the permeability and ultimately increases the risk of subsidence. Hoc Yakupoğlu et al, 2019 [43]. The higher the density of these fractures, the greater the permeability, Semeraro et al, 2016 [33].Given the size of the study area, it is virtually impossible to locate and measure fractures. However, since the presence of fractures depends on the fault, on a large scale, by locating the faults and calculating their density, the highest fracture density can be logically achieved. To prepare this layer, the faults drawn in the geological maps of Khuzestan province were prepared with the faults drawn using integration satellite images and the distribution map of the faults of the province (Fig. 8). Finally, using the fault layer and applying the Density function in ArcGIS environment, a fault density map of Khuzestan province was prepared (Fig. 9).

Layer standardization or fuzzy
In the eld of GIS, to determine subsidence potential, standardize GIS-based raster criteria by assigning a value between [1 − 0] to each pixel using fuzzy membership functions such as Sigmoidal, J shape, Linear, or many other forms. The complexity is non-uniform,Carter et al, 2021 [13].In this study, standardization of data according to their quantitative and qualitative nature was done in two ways. Quantitative data were standardized using fuzzy membership functions and qualitative data were standardized using the rasterization method and assigning a value between [0-1].

Standardization of quantitative layers
As mentioned in the section above, quantitative layer standardization was performed using fuzzy membership functions. Quantitative layers include slope, fault density, and waterway density (Figs. 10 and 11).

Standardization of quality layers
Standardization of qualitative layers (geological layer) has been done using rasteration. In this method, values between zero and one [1 − 0] were given to the geological map and nally converted to raster using these values.

Overlap
Overlap is a spatial function that can combine layers of spatial data obtained from separate sources for location applications using hybrid models. The new layer (output) is a function of two or more input layers.
Due to the fact that effective criteria and sub-criteria have different validities and all of them must participate in the overlap, the index overlap method has been used for this purpose. In this method, the standardized layer obtained from each criterion (x_i) is multiplied by its validity (its weight) by the criterion (w_j). This is done for all criteria and sub-criteria and new layers are obtained that are overlapped using various functions such as AND, OR, SUM, Product and GAMA.
A i = ∑ w j x ij , ∑ w j = 1 After overlapping the main criteria with different functions, nally the output of a function is considered as the best option, in accordance with eld and statistical evidence, for zoning of ood potential. At this stage, the output of AHP method is considered as the nal output. Is (Fig. 12).

Fuzzy inference system programming in MATLAB environment:
The fuzzy inference system is an intelligent system whose task is to learn expertise from experts and to be able to make decisions based on the expertise it has learned Munier and Hontoria, 2021 [26]. In this system, the aim is to de ne a system that determines the risk level of the project based on different information layers,Singh and Lone, 2021 [34].
In these studies, based on the value and validity of each information layer, a speci c weight for that layer in MATLAB program is introduced and nally, based on the de ned weight, the level of risk in the study area is determined, the most important of which include the following:

Introduction and initial validation of layers:
Lithology layer: As mentioned before, this layer is the most important layer affecting the results, so with high accuracy, the whole of Khuzestan province was divided into four types of evaporitic lithology, carbonate, marl and shale and nally detrital, and the effect of each lithology with a code 1 to 4 (lithology, respectively) were introduced.
Layers of distance from the fault and distance from the waterway: The closer the distance from the fault and the distance from the waterway, the greater the risk, so from two scales of distance less than 300 meters (high risk) and more than 300 meters (low risk) for risk It has been used in these studies.
Slope layer: The degree of risk based on the slope in these studies were divided into three categories: 0 to 20 degrees (low risk), 20 to 70 degrees (steep and high risk), 70 to 90 degrees (low slope and low risk).

Fuzzy input variables (fasi cation):
The attribution function is used to fuzzy the 4 input layers because it becomes a fuzzy function when its attribution function is known and its linguistic variable is knownHsu et al, 2021 [20].Therefore, at this stage, after the initial segmentation and validation of the layers, with the help of the fuzzy part of the MATLAB program, it has been used to de ne the different rules of these layers.In these studies, Mamdani FIS Type is used, which is the most widely used fuzzy inference system used in public works.
Mamdani fuzzy inference was rst introduced as a method to create a control system by synthesizing a set of linguistic control rules obtained from experienced human operators [1]. In a Mamdani system, the output of each rule is a fuzzy set, Mamdani and Assilian, 1975[15].Since Mamdani systems have more intuitive and easier to understand rule bases, they are well-suited to expert system applications where the rules are created from human expert knowledge, such as Geological diagnostics Darbor et al, 2020. [14]The inputs along with the output were divided into 5 sections using a triangular function (Fig. 13).
The intensity of changes in these intervals is indicated by the letters ( VL, L, M, H, VH)- (Fig. 14). For lithology, 4 types of change intensities including detrital (V-LOW), marl (LOW), carbonate (MED) and evaporative (High) were used. For the distance from the fault and the distance from the waterway, high change intensity (0 to 300 m) and low change intensity (300 to 1000 m) were used .For slope, the intensity of change was very low (0 to 20 degrees), high (20 to 70) and low (70 to 90).

De nition of rules:
Duction of correct and practical rules plays a very important role in the nal decision of the care system, Harahap et al, 2021 [19]. These rules must be made very carefully and with the help of the opinions of experienced experts, Yanar et al, 2020 [45]. In these studies, 49 laws were introduced for this system according to the opinions of different experts (Fig. 15).

Import inference system and inputs in MATLAB application using Eval s function:
With the help of this function, general inputs can be introduced to the fuzzy system and the output of all of them can be received, Akgun et al, 2021 [5].In fact eval s ( s,input) evaluates the fuzzy inference system s for the input values in input and returns the resulting output values in output,Yamagishiand Bhandary, 2017 [44].
According to Fig. 16, the lithology entrance is marked with an L, the distance from the waterway is marked with an R, the distance from the fault is marked with an F, and the slope is entered with a mark D. As can be seen, by entering the exact values of the study area, the degree of risk of the study area in eld operations can be determined. According to Fig. 16, if the type of Avari lithology (1) and the distance from the river is 600 meters and the slope of the layers is 80 degrees and the distance from the fault is 400 meters, the output is 0.6, which according to the de nitions for this number system This is the study section ( Fig. 16).Also, according to Fig. 17, if the lithology of the area is carbonate (number 3) and the distance from the river is 100 meters and the distance from the fault is 80 meters and the slope is 60 degrees, the risk level of this area is 2.5, i.e., medium to weak risk. The results seem quite reasonable considering that the amount of hazards in carbonate layers is more than detrital layers. It should be noted that the numerical scales of lithology include 1: detrital, 2: marl and shale, 3: carbonate and 4: evaporitic, respectively.This method will be very useful in preparing geological maps with scales of 1: 5000 to 1: 25000 and also in preparing maps of high-risk geological areas that are prepared with the aim of carefully studying a small part of the earth.

Validation of results with eld and statistical evidence
In order to interpret, analyze and validate the prepared zoning, 2000 points from different regions of the province that have statistics and evidence of subsidence were evaluated on the output of the prepared subsidence map.
The results of this study show that the outputs of this research are more than 90% close to reality. These studies have been carried out in different cities of Khuzestan province that have low to high subsidence potentials. In the lower part of Gotvand city in the north of Khuzestan province is presented as an example of eld studies and telemetry.
Gotvand city is located in the southwest of Iran and in the north of Khuzestan province and on the edge of Khuzestan plain. In the northwest of this city, there are heights of Bakhtiari Formations (with conglomerate lithology), Aghajari (with sandstone lithology), Gachsaran (with evaporitic lithology) and Mishan (with Marne lithology). They show a very high agreement with the results of the study of the two proposed models (Fig. 18). Field studies conducted in this area showed that in some parts of the Gachsaran Formation, which are close to seasonal and permanent waterways, several sinkholes are being formed (Fig. 19).
In this region, with the proximity to the Avari and Marni formations (Aghajari and Mishan), the amount of subsidence has been reduced so that in these formations, despite the proximity to the waterway and different slopes, no subsidence effects are observed (Fig. 20).
As stated, the results of both models presented to investigate the potential for subsidence in this area show a very good agreement with eld observations so that in parts of the Gachsaran Formation that are less distant than waterways and faults According to these models, they are less dangerous than the areas that are closer to these complications.
According to the models introduced for this region, Bakhtiari, Aghajari and Mishan formations are divided as very low risk areas and Gachsaran formation is divided into two high-risk and very high-risk areas.

Conclusion
The results of the study of areas with subsidence potential show very good results for the study area. According to these results, the highest risk of subsidence is observed in the northern and northwestern parts of Khuzestan province. In these areas, due to the presence of evaporitic formations such as Gachsaran Formation, the presence of many waterways, slope changes and the presence of many faults, the greatest potential for subsidence hazards is observed.These studies for the rst time introduced a program in MATLAB program that various experts in earth sciences, especially engineering geology with lithological information, distance from waterway, distance from fault and slope, the ability to access the risky area before any Civil operations or other geological studies.
This program is designed to be upgradeable and updated over time.The study of land subsidence in the study areas using the AHP method had similar results to the method introduced in the MATLAB program, which corresponds to the two methods with eld observations indicate the accuracy of over 90% of these methods.Comparison of the two methods shows that the programming method for areas with a very small area has a much higher accuracy that this method can be used in important civil studies. Finally, for the rst time, the subsidence map with the methods introduced for the study area of southwestern Iran was introduced as a model for this type of study.

Figure 1
Location of the study area in southwestern Iran (Khuzestan province).   Slope map of Khuzestan province  Density map of faults in Khuzestan province Figure 10 Standard map of the slope layer of Khuzestan province Figure 11 Standard map of the fault density layer of Khuzestan province  General shape of the inputs and outputs in the fuzzy logic method with Mamdani algorithm in the study area Figure 14 The intensity of the changes made for each layer and their output Figure 15 A picture of how to introduce different rules in the MATLAB program Figure 16 An example of studies performed to validate the programming method Figure 17 An example of studies performed to validate the programming method Figure 18 Location and name of existing formations in the veri ed area north of Khuzestan province Figure 19 Examples of subsidence of Gachsaran Formation in the studied areas A) in the east of Gotvand city, B) in the northeast of Gotvand city in the southwest of Iran.