Mapping Risk to Land Subsidence: Developing a Two-Level Modeling Strategy by Combining Multi-Criteria Decision-Making and Artiﬁcial Intelligence Techniques

: Groundwater over-abstraction may cause land subsidence (LS), and the LS mapping suffers the subjectivity associated with expert judgment. The paper seeks to reduce the subjectivity associated with the hazard, vulnerability, and risk mapping by formulating an inclusive multiple modeling (IMM), which combines two common approaches of multi-criteria decision-making (MCDM) at Level 1 and artiﬁcial intelligence (AI) at Level 2. Fuzzy catastrophe scheme (FCS) is used as MCDM, and support vector machine (SVM) is employed as AI. The developed methodology is applied in Iran’s Tasuj plain, which has experienced groundwater depletion. The result highlights hotspots within the study area in terms of hazard, vulnerability, and risk. According to the receiver operating characteristic and the area under curve (AUC), signiﬁcant signals are identiﬁed at both levels; however, IMM increases the modeling performance from Level 1 to Level 2, as a result of its multiple modeling capabilities. In addition, the AUC values indicate that LS in the study area is caused by intrinsic vulnerability rather than man-made hazards. Still, the hazard plays the triggering role in the risk realization.


Introduction
Land subsidence (LS) due to groundwater over-abstraction is a man-made problem, which threatens water availability, the environment, and structures. Risk assessment, as the first step of risk management, can play a pivotal role in proactively managing or mitigating LS problems. Risk assessment can be carried out by risk indexing or mapping when historical records are unavailable, and, consequently, frequency analysis is not feasible. Different techniques are available in the ongoing research activities to formulate risk indexing or mapping by incorporating a set of data layers, as discussed in due course. However, they suffer from subjectivity associated with the data layers due to expert judgment. The paper aims to delineate risk maps to LS and reduce the subjectivity by incorporating multicriteria decision-making (MCDM) and artificial intelligence (AI) techniques, two widely used approaches in hydrology and environmental studies.
The literature review highlights two different approaches of MCDM and AI in risk mapping of groundwater issues. MCDM formulates the issue as a decision-making problem Both the MCDM and AI techniques require the measure LS values as validation or target datasets. Generally, LS measurement methods are classified into direct and indirect methods. In the direct method, global positioning system (GPS) and accurate alignment are used to detect LS. These methods are accurate but expensive due to installation and maintenance costs and covering only restricted areas. The indirect method does not have these limitations and employs interferometric synthetic aperture radar (InSAR) as the remote sensing technology. Different techniques for the InSAR analysis are available in the literate including, differential InSAR (D-InSAR) [23], ALOS InSAR [24,25], small baseline (SB) [26], and persistent scatterer InSAR (PSI) [27,28].
The paper aspires to decrease the subjectivity in LS mapping by combining both MCDM and AI techniques, which are discussed in the literature review. The paper's formulation calculates the vulnerability index based on the MCDM techniques at the first level. It then incorporates the MCDM results as a target dataset in an AI-based formulation at the second level. Notably, FCS is used as an MCDM technique at the first level, and SVM is used as the AI technique at the second level. The paper also transfers vulnerability indexing to risk indexing by multiplying the hazard and vulnerability data layers.

Methodology
In this paper, a multiple modeling strategy, called inclusive multiple modeling (IMM), is formulated to reduce the subjectivity in rates and weights of the ALPRIFT data layers, and thereby improve modeling performance. Previous studies have outlined a mathematical basis that demonstrates the reduced error rate of multi-modeling compared to singlemodeling [29,30]. On this basis, IMM formulates the mapping of hazard, vulnerability, and Water 2021, 13, 2622 4 of 15 risk to LS. At Level 1, FCS is employed to map the hazard, vulnerability, and risk, using a set of incorporated data layers, referred to as ALPRIFT data layers suggested by [12]. At Level 2, SVM is used to map the supervised hazard, vulnerability, and risk. SVM is trained using ALPRIFT data layers as input data and a data-fused of LS map with Level 1 results as the target data. The LS map is prepared using interferometric synthetic-aperture radar (InSAR) processing and data fusion performed by catastrophe theory. The FCS advantages include the ability to rely on the MCDM technique and SVM takes advantage of artificial intelligence (AI) models to employ supervised learning capabilities. By combining both capabilities through IMM, the paper aims to improve modeling performance.

Basic ALPRIF Framework
The ALPRIFT framework was developed by [12] for mapping the vulnerability of LS by including seven data layers, such as aquifer media (A), land use (L), pumping of groundwater (P), recharge (R), aquifer thickness impact (I), fault distance (F), and decline of water table (T). These data layers are rated and weighted as per prescribed values to calculate subsidence vulnerability index (SVI) as follows, where the subscripts w and r represent weight and rate, respectively.
Transforming Vulnerability to Risk The risk concept in a system covers both concepts of hazard and vulnerability, in which hazard refers to actuating factor that triggers a risk; and vulnerability refers to the system's resistance. In the LS problem, P and T data layers are time-variant that triggers the risk, and A, L, R, I, F data layers are time-invariant and refer to the intrinsic vulnerability of an aquifer. Table 2 presents the required input datasets and GIS-processing steps for preparing hazard and vulnerability data layers. According to the above formulation of risk, the subsidence risk index (SRI) is calculated by the product of hazard and vulnerability as follows [32]:

InSAR Processing
Based on the difference in the phase between two Sentinel-1 SAR observations, InSAR processes obtain information about the earth's surface. InSAR processing leverages the amount of phase change between two consecutive and complete sine wave cycles. In particular, the phase follows the topography of the terrain. From the images and phase information in SAR systems, we can determine the strength of radar recognition from the amplitude information. The study used interferometric wide (IW) swath products with a spatial resolution of 5 m × 20 m and a swath of 250 km. A total of three sub-swaths are available here, based on terrain observation with progressive scans SAR (TOPSAR). A uniform signal-to-noise ratio (SNR) is provided along with a distributed target ambiguity ratio (DTAR) to produce homogeneous image quality throughout the swath. Several procedures are necessary to process InSAR data, including co-registration, interferogram formation, coherence estimation, phase removal from topography, phase filtering, phase unwrapping, terrain correction, and converting displacement along line-of-sight to vertical displacement. Further information about this procedure is available by [33].

Modeling Strategy at Level 1 by Fuzzy Catastrophe Scheme
FCS is used in the paper for mapping hazard, vulnerability, and risk at Level 1, where FCS reduces subjectivity in rates and weights of APLRIFT data layers. In MCDM, FCS is a technique developed by [34], bringing fuzzy membership analysis and catastrophe theory together. In fuzzy logic [35], membership functions, fuzzy sets, and fuzzy inference engines are used. According to [36], a catastrophe is the result of a set of dependent variables (also called state variables) and independent variables (also called control parameters). Recently, FCS was utilized in different fields of water resources, including ecological environment sustainability by [37], predicting harmful algae blooms by [38], aquifer vulnerability to LS by [6], and aquifer vulnerability to saltwater intrusion by [5].
As part of the catastrophe theory, a ranked list of functions can be assigned to the data layers, which include fold, cusp, swallowtail, butterfly, and wigwam. A state variable and control parameters ranging from 1 to 5 define these functions. Table 3 represents different types of catastrophe functions and related control parameters. Table 3. Catastrophe functions [6,34]. The formulated FCS includes: (i) hazard with two control parameters (PT data layers) then uses cusp, and (ii) vulnerability with five control parameters (ALRIF data layers) then uses wigwam. Therefore, FCS specifies the following functions for both hazard, vulnerability, and risk:

Name State Variable Control Parameter Catastrophe Fuzzy Membership Functions
where Equations (4) and (5) contains two and five data layers, respectively. Data layers with higher weight gain the higher power, according to the recommendation of [12]. The mean operator calculates the system status based on the complementary principle; the alternative approach may be the minimum operator, which does not arise in the problem of LS. Notably, data layers in Equation (3) are normalized between 0 and 1 by linear membership functions as per the following equations: where i counts pixels; X max and X min are maximum and minimum values, respectively; and X n i normalizes particular values at the ith pixels. Notably, Equation (7) normalizes data layers that are directly proportional to LS (A, L, P, I, T data layers); and Equation (8) normalizes data layers with inverse proportion (R, F data layers).

Modeling Strategy at Level 2 by SVM
SVM at Level 2 incorporates the Level 1 results in its structure. The input and target datasets of SVM is formulated as follows: (i) ALPRIFT data layers are input dataset; (ii) the data-fused of Level 1 result and InSAR subsidence map by catastrophe theory is considered as target datasets. Three SVMs are trained and tested for mapping supervised maps of hazard, vulnerability, and risk by SVM1, SVM2, and SVM3, respectively. Figure 1 summarizes the details in the architecture of SVMs. The target datasets for SVM1-SVM3 are calculated through Equations (9)-(11) based on the cusp catastrophe equation (see Table 3) as follows:

Performance Metrics
The performance of hazard, vulnerability, and risk indices concerning the InSAR result is evaluated by the receiver operating characteristic (ROC) curve and the area under curve (AUC). These criteria were developed by [31] to evaluate the precision of a diagnostic system, where the diagnosis refers to Earth's displacement in the LS problem. The events related to diagnosis are divided into four groups of true negative (TN), true positive (TP), false positive (FP), and false negative (FN). Concerning the threshold of various settings, the ROC curve plots FP proportion against the TP proportion. An upper left corner deviation of the ROC curve corresponds to an undesirable performance. AUC measures the ratio of the area under the ROC curve compared to the total area ranging from 0.5 to 1. The AUC values equal to 1 and 0.5 correspond to the perfect and random performance, respectively.
The determination coefficient (R 2 ) and root mean square error (RMSE) metrics evaluate the performance of SVM models. A value of R 2 close to 1 represents the perfect performance, and a value close to zero shows unsatisfactory performance. Additionally, RMSE values close to 0 represent the best models, and it has no upper limit for unsatisfactory performance. The determination coefficient (R 2 ) and root mean square error (RMSE) metrics evaluate the performance of SVM models. A value of R 2 close to 1 represents the perfect performance, and a value close to zero shows unsatisfactory performance. Additionally, RMSE values close to 0 represent the best models, and it has no upper limit for unsatisfactory performance.

Study Area
This study is being conducted in the Tasuj plain, which lies with the province of East Azerbaijan in northwest Iran (see Figure 2). Additionally, this plain lies on the northern shores of Lake Urmia, where water levels have declined by about 15 m since 2000. At Tasuj climatological station, the average annual precipitation is 290 mm from 2008-2018, characterized by [39] as an arid and cold climate. The maximum and minimum temperatures are 33°C and −11°C during the same period. The area of the plain is about 300 km 2 and generally is composed of alluvial sediments. In addition to the sedimentary deposit, igneous and metamorphosed formations are observable from Precambrian to the recent era. Figure 2 illustrates the lithological formations and the spatial location of the fault lines.
The sedimentary deposits form an unconfined aquifer as the primary source of water supply for agricultural and domestic purposes in the plain. The groundwater and surface water flow direction is towards Lake Urmia (see the flow direction in Figure 2b). The surface water is seasonal and is not a reliable source of water supply. The groundwater is discharged by 144 tube wells, 40 springs, and 70 qanats, which is equivalent to 16 × 10 6 m 3 volume of water annually [40]. The groundwater level (GWL) in the aquifer is monitored monthly by 28 observation wells. The locations of the observation wells and abstraction wells are given in Figure 2b,c. The evaluation of the GWL time series in these wells shows that the GWL declined is 8.6 cm per month, which can trigger environmental problems such as LS. The paper presents the spatial distribution of subsidence captured by InSAR from 2017 to 2018 in due course.

Study Area
This study is being conducted in the Tasuj plain, which lies with the province of East Azerbaijan in northwest Iran (see Figure 2). Additionally, this plain lies on the northern shores of Lake Urmia, where water levels have declined by about 15 m since 2000. At Tasuj climatological station, the average annual precipitation is 290 mm from 2008-2018, characterized by [39] as an arid and cold climate. The maximum and minimum temperatures are 33 • C and −11 • C during the same period. The area of the plain is about 300 km 2 and generally is composed of alluvial sediments. In addition to the sedimentary deposit, igneous and metamorphosed formations are observable from Precambrian to the recent era. Figure 2 illustrates the lithological formations and the spatial location of the fault lines.
The sedimentary deposits form an unconfined aquifer as the primary source of water supply for agricultural and domestic purposes in the plain. The groundwater and surface water flow direction is towards Lake Urmia (see the flow direction in Figure 2b). The surface water is seasonal and is not a reliable source of water supply. The groundwater is discharged by 144 tube wells, 40 springs, and 70 qanats, which is equivalent to 16 × 10 6 m 3 volume of water annually [40]. The groundwater level (GWL) in the aquifer is monitored monthly by 28 observation wells. The locations of the observation wells and abstraction wells are given in Figure 2b,c. The evaluation of the GWL time series in these wells shows that the GWL declined is 8.6 cm per month, which can trigger environmental problems such as LS. The paper presents the spatial distribution of subsidence captured by InSAR from 2017 to 2018 in due course.

Preparation of ALPRIFT Data Layer
The data layers related to hazard (PT) and vulnerability (ALRIF) are prepared after pre-processing and GIS-processing as per Table 2. Figure 3 illustrates the spatial distribution of these data layers. A plot of the P-and T-data layers are depicted in Figure 3a,b based on pumping volume in abstraction wells and decline trend in observation wells, interpolated by the inverse distance weighting (IDW) technique.
In Figure 3c, aquifer media are calculated in geological logs as per aquifer media rates recommended by [12] and interpolated by IDW. The land use data layer is prepared by image processing through ENVI software. The required procedures are correcting the geometric and atmospheric, identifying different land using normalized difference vegetation index (NDVI), and interpreting images by the supervised classification approach to

Preparation of ALPRIFT Data Layer
The data layers related to hazard (PT) and vulnerability (ALRIF) are prepared after preprocessing and GIS-processing as per Table 2. Figure 3 illustrates the spatial distribution of these data layers. A plot of the P-and T-data layers are depicted in Figure 3a,b based on pumping volume in abstraction wells and decline trend in observation wells, interpolated by the inverse distance weighting (IDW) technique.
classify different land uses based on the maximum likelihood method. Figure 3d represents the land use data layer within the study area. The recharge data layer is calculated as per [41] using precipitation, infiltration, and slope data, and the result is presented in Figure 3e. The thickness of the aquifer in Figure 3f is calculated by geological logs and interpolated by IDW. In Figure 3g, a fault distance is calculated using the Euclidian distance toolbox available in the ArcGIS software.  (1) and the prescribed rates and weights by [12]. The comparison between Figure 4a,b indicates that the result of basic ALPRIFT is not entirely defensible, and there is room for improvement. Figure 4c,d represent, respectively, the hazard and vulnerability indices as per Equations (4) and (5). There are some similarities and differences between the spatial pattern of hazard and vulnerability indices, which is expected and stems from the incorporated data layers. The risk index pools together these similarities and differences as per Equation (6). Notably, hazard triggers a risk in areas with higher vulnerability, and this issue is illustrated in Figure 4e. A visual comparison between the InSAR result ( Figure 4a) and the risk index (Figure 4e) indicates higher agreement than basic ALPRIFT, but there is still room for improvement. In Figure 3c, aquifer media are calculated in geological logs as per aquifer media rates recommended by [12] and interpolated by IDW. The land use data layer is prepared by image processing through ENVI software. The required procedures are correcting the geometric and atmospheric, identifying different land using normalized difference vegetation index (NDVI), and interpreting images by the supervised classification approach to classify different land uses based on the maximum likelihood method. Figure 3d represents the land use data layer within the study area. The recharge data layer is calculated as per [41] using precipitation, infiltration, and slope data, and the result is presented in Figure 3e. The thickness of the aquifer in Figure 3f is calculated by geological logs and interpolated by IDW. In Figure 3g, a fault distance is calculated using the Euclidian distance toolbox available in the ArcGIS software. Figure 4a illustrates the result of LS by InSAR processing from 2017 to 2018. The figure indicates that LS is mainly concentrated in the south-central part of the study area despite the scattered LS records elsewhere. The result of basic ALPRIFT is given in Figure 4b as per Equation (1) and the prescribed rates and weights by [12]. The comparison between Figure 4a,b indicates that the result of basic ALPRIFT is not entirely defensible, and there is room for improvement. Figure 4c,d represent, respectively, the hazard and vulnerability indices as per Equations (4) and (5). There are some similarities and differences between the spatial pattern of hazard and vulnerability indices, which is expected and stems from the incorporated data layers. The risk index pools together these similarities and differences as per Equation (6). Notably, hazard triggers a risk in areas with higher vulnerability, and this issue is illustrated in Figure 4e

Hazard, Vulnerability, and Risk Indices at Level 2
The hazard, vulnerability, and risk indices are calculated as per supervised learning by SVM1, SVM2, and SVM3, described in detail in Figure 1. Table 4 presents the performance criteria and SVM parameters for the hazard, vulnerability, and risk indices. The table evaluates the results regarding R 2 and RMSE metrics for training, testing, and total datasets. According to the table, the models have similar performance in terms of both metrics in the training and testing phases. The spatial distribution of the results is given in Figure 5. Although there are similarities between this figure and the corresponding results in Figure 4, the results at Level 2 provide evidence that the Level 2 results are more compatible with the result of InSAR. The results at Level 2 identify the southcentral part of the study area as the hotspot area in terms of hazard, vulnerability, and risk.
The risk map in Figure 5c represents the hotspots with higher priorities for modifying management policies. However, this figure provides a relatively conservative result compared to the results of the hazard and vulnerability because the areas swept by Band 1 show a marked increase in risk potential. Additionally, some information in the hazard and vulnerability indices is not reflected in the risk index. For example, there are areas with a higher hazard index but a relatively lower vulnerability index in the north part of the plain. Thus, although human activities threaten these areas, they tend to be relatively less vulnerable, leading to a lower risk index.

Hazard, Vulnerability, and Risk Indices at Level 2
The hazard, vulnerability, and risk indices are calculated as per supervised learning by SVM1, SVM2, and SVM3, described in detail in Figure 1. Table 4 presents the performance criteria and SVM parameters for the hazard, vulnerability, and risk indices. The table evaluates the results regarding R 2 and RMSE metrics for training, testing, and total datasets. According to the table, the models have similar performance in terms of both metrics in the training and testing phases. The spatial distribution of the results is given in Figure 5. Although there are similarities between this figure and the corresponding results in Figure 4, the results at Level 2 provide evidence that the Level 2 results are more compatible with the result of InSAR. The results at Level 2 identify the southcentral part of the study area as the hotspot area in terms of hazard, vulnerability, and risk.    Figure 6 evaluates the indices of basic ALPRIFT, hazard, vulnerability, and risk at Level 1 ( Figure 6a) and Level 2 (Figure 6b) by the ROC curve and the AUC value. The figure indicates that the basic ALPRIFT provides noisy signals with the closest distance to the random classifier (the diagonal line). The AUC value for basic ALPRIFT is 0.56, which denotes that there is room for improvement. The ROC curve for hazard indices at Level 1 is slightly away from the random classifier line with an AUC value of 0.63. A similar improvement also is observed for the vulnerability and risk indices at Level 1, respectively, with AUC values of 0.72 and 0.71. This improvement refers to decreasing the subjectivities with rates and weights by FCS.

Evaluation of Results in Terms of ROC and AUC
The results at Level 2 provide considerable improvements compared with the results at Level 1 as follows: (i) the AUC values for the hazard increases from 0.63 to 0.67; (ii) the AUC values for the vulnerability increases from 0.72 to 0.80; and (iii) the AUC values for the risk increases from 0.71 to 0.82. This improvement is expected due to using a learning technique by SVM. The lower AUC values for the hazard indices compared to vulnerability at both levels implies that the LS occurrence is more affected by intrinsic vulnerability than a man-made hazard. However, the hazard has a triggering role to risk since the ROC curve for the risk is somewhat higher than the vulnerability, and the AUC value for the risk index is slightly higher than the vulnerability index. The risk map in Figure 5c represents the hotspots with higher priorities for modifying management policies. However, this figure provides a relatively conservative result compared to the results of the hazard and vulnerability because the areas swept by Band 1 show a marked increase in risk potential. Additionally, some information in the hazard and vulnerability indices is not reflected in the risk index. For example, there are areas with a higher hazard index but a relatively lower vulnerability index in the north part of the plain. Thus, although human activities threaten these areas, they tend to be relatively less vulnerable, leading to a lower risk index.

Discussion
ALPRIFT and similar frameworks are data-driven, and their results rely on data and the characteristics of study areas. These frameworks are in their infancy, and further investigations, such as using different modeling strategies, are required for proofing the concept. As discussed in Section 1, both MCDM and AI techniques were implemented on ALPRIFT in the previous studies. However, the paper combined two common approaches, and, in fact, it extended the MCDM-based study by [6] through considering a new level, which incorporates the Level 1 result as the target dataset in the formulation of Level 2.
The paper utilized a two-level multiple modeling strategy, referred to as inclusive multiple modeling (IMM), and the previous studies (e.g., [11]) indicated that IMM could significantly improve the modeling performance. As discussed earlier, the paper selected FCS at Level 1 and SVM at Level 2. There is no theoretical basis for model selection at both levels, and, generally, they are selected by a trial-and-error procedure. The idea for employing an MCDM technique at Level 1 and an AI technique at Level 2 stems from different formulations in the literature, which uses only AI approaches at both levels. However, FCS and SVM can be substituted with other techniques of MCDM and AI in future studies.
There is a significant difference between the AUC values related to the hazard and vulnerability maps. The paper identified that LS occurrence is more affected by intrinsic vulnerability than man-made hazard as per ROC/AUC. However, it should be noted that the man-made hazard triggers the LS, and it does not occur in the absence of over-abstraction, even with system vulnerabilities. Sadeghfam et al. [6] identified the man-made data layers as more effective than intrinsic data layers. This is expected because this issue is attributable to the characteristics of study areas. However, quantifying the role of these components provides more insight into risk management. Further investigation by statistical techniques is highly recommended to be undertaken to precisely identify the impact percentage of all components by future studies.
The investigation revealed that the hotspots identified by FCS at Level 1 occupy more areas compared to the hotspots by SVM. This is traceable to two issues: (i) The extent of areas with LS records is limited in the study area; and (ii) AI techniques are more powerful techniques in making a relationship between the input data layers and LS captured by InSAR even though the target dataset is fused by the Level 1 results. There seems to be room for investigating the role of others data fusion techniques in preparing target datasets in future studies. The results at Level 2 provide considerable improvements compared with the results at Level 1 as follows: (i) the AUC values for the hazard increases from 0.63 to 0.67; (ii) the AUC values for the vulnerability increases from 0.72 to 0.80; and (iii) the AUC values for the risk increases from 0.71 to 0.82. This improvement is expected due to using a learning technique by SVM. The lower AUC values for the hazard indices compared to vulnerability at both levels implies that the LS occurrence is more affected by intrinsic vulnerability than a man-made hazard. However, the hazard has a triggering role to risk since the ROC curve for the risk is somewhat higher than the vulnerability, and the AUC value for the risk index is slightly higher than the vulnerability index.

Discussion
ALPRIFT and similar frameworks are data-driven, and their results rely on data and the characteristics of study areas. These frameworks are in their infancy, and further investigations, such as using different modeling strategies, are required for proofing the concept. As discussed in Section 1, both MCDM and AI techniques were implemented on ALPRIFT in the previous studies. However, the paper combined two common approaches, and, in fact, it extended the MCDM-based study by [6] through considering a new level, which incorporates the Level 1 result as the target dataset in the formulation of Level 2.
The paper utilized a two-level multiple modeling strategy, referred to as inclusive multiple modeling (IMM), and the previous studies (e.g., [11]) indicated that IMM could significantly improve the modeling performance. As discussed earlier, the paper selected FCS at Level 1 and SVM at Level 2. There is no theoretical basis for model selection at both levels, and, generally, they are selected by a trial-and-error procedure. The idea for employing an MCDM technique at Level 1 and an AI technique at Level 2 stems from different formulations in the literature, which uses only AI approaches at both levels. However, FCS and SVM can be substituted with other techniques of MCDM and AI in future studies.
There is a significant difference between the AUC values related to the hazard and vulnerability maps. The paper identified that LS occurrence is more affected by intrinsic vulnerability than man-made hazard as per ROC/AUC. However, it should be noted that the man-made hazard triggers the LS, and it does not occur in the absence of over-abstraction, even with system vulnerabilities. Sadeghfam et al. [6] identified the man-made data layers as more effective than intrinsic data layers. This is expected because this issue is attributable to the characteristics of study areas. However, quantifying the role of these components provides more insight into risk management. Further investigation by statistical techniques is highly recommended to be undertaken to precisely identify the impact percentage of all components by future studies.
The investigation revealed that the hotspots identified by FCS at Level 1 occupy more areas compared to the hotspots by SVM. This is traceable to two issues: (i) The extent of areas with LS records is limited in the study area; and (ii) AI techniques are more powerful techniques in making a relationship between the input data layers and LS captured by InSAR even though the target dataset is fused by the Level 1 results. There seems to be room for investigating the role of others data fusion techniques in preparing target datasets in future studies.

Conclusions
ALPRIFT as a standard framework quantifies aquifers' vulnerability to land subsidence (LS), but it suffers subjectivity associated with the incorporated data layers. The paper formulated a methodology based on inclusive multiple modeling (IMM) at two levels to decrease the subjectivity and increase the modeling performance. IMM also map the hazard, vulnerability, and risk indices to LS, in which a fuzzy catastrophe scheme as a multi-criteria decision making (MCDM) technique was used at Level 1; and the obtained results were feed as the target dataset to the support vector machine as an artificial intelligence (AI) technique at Level 2. Therefore, the formulated methodology combines two capabilities of MCDM and AI. The formulation was implemented in Tasuj plain, located in northwest Iran, which suffers from over-abstraction. Results at both levels identify significant signals as per the receiver operating characteristic and the area under curve (AUC) performance metrics for the hazard, vulnerability, and risk indices. Additionally, further improvements were achieved from Levels 1 to 2 for the vulnerability and risk indices, which is expected due to the capability of IMM. However, there is room for further improvements in the performance of results by future studies. The higher AUC values for the vulnerability index on recorded LS compared to the hazard index indicated that the LS occurrence in the study area is more affected by the intrinsic vulnerability of the system, but anthropogenic activities play an actuation role in the risk realization.