Machine Learning-Based Assessment of Watershed Morphometry in Makran

Abstract

This study proposes an artificial intelligence approach to assess watershed morphometry in the Makran subduction zones of South Iran and Pakistan. The approach integrates machine learning algorithms, including artificial neural networks (ANN), support vector regression (SVR), and multivariate linear regression (MLR), on a single platform. The study area was analyzed by extracting watersheds from a Digital Elevation Model (DEM) and calculating eight morphometric indices. The morphometric parameters were normalized using fuzzy membership functions to improve accuracy. The performance of the machine learning algorithms is evaluated by mean squared error (MSE), mean absolute error (MAE), and correlation coefficient (R²) between the output of the method and the actual dataset. The ANN model demonstrated high accuracy with an R² value of 0.974, MSE of 4.14 × 10⁻⁶, and MAE of 0.0015. The results of the machine learning algorithms were compared to the tectonic characteristics of the area, indicating the potential for utilizing the ANN algorithm in similar investigations. This approach offers a novel way to assess watershed morphometry using ML techniques, which may have advantages over other approaches.

1. Introduction

Watershed morphometry is a crucial factor in determining the impact of tectonic processes on the landscape. By analyzing the shape and geometry of watersheds at a regional scale, we can identify the relative significance of tectonic deformation versus erosion in landscape evolution [1,2]. Understanding the impact of these geological forces on the morphology of watersheds and the development of drainage systems is essential, as it can have implications for sediment supply to river reaches and increase the risk of landslides [3]. Furthermore, by quantifying the morphotectonic situation of watersheds, we can gain insight into their evolution and assess the role of regional tectonic control in shaping their development [2,4,5,6]. Recent studies (e.g., [7,8,9,10]) have expanded our understanding of the relationship between morphotectonic factors and the development of drainage systems and landscapes, building upon the foundational work of Horton [11], Strahler [12], and Hack [13]. Through a quantitative assessment of watershed development in relation to regional tectonics, we can better understand the morphotectonic situations and their implications for the broader landscape [14].

The utilization of Geographic Information System (GIS), satellite remote sensing data processing, particularly Global Digital Elevation Models (DEMs), along with the application of Analytic Hierarchy Process (AHP) and Fuzzy Analytic Hierarchy Process (FAHP) techniques, can facilitate the determination of comprehensive geomorphometric analyses (e.g., [15,16,17]).

Due to the inhomogeneity of geological landscapes and vast data, we need complicated models and methods to study geological features. Numerous complex mathematical approaches have been implemented to overcome this issue. To this end, multiple types of research have been conducted using AHP, FAHP, Technique for Order Preference by Similarity to Ideal Solution (TOPSIS), and machine learning algorithms. Some of them are combined (hybrid methods) to take advantage of the best parts of each and obtain the best results [18,19,20].

The use of a Geographic Information System (GIS), a digital database management system specifically designed for handling large-scale and spatially diverse data from various sources, has the potential to reveal the drainage patterns of watersheds [21,22] and is well-suited for advanced zoning applications.

A notable benefit of the Analytic Hierarchy Process (AHP) is its foundation in pairwise comparisons, which effectively facilitate the computation of criteria weights. In addition, the AHP calculates inconsistency indices; ratios of a decision maker’s inconsistency. Nevertheless, decision makers must occasionally perform a huge number of pairwise comparisons, and this condition, particularly with fuzzy AHP, makes using the AHP procedure impracticable [23]. When the researcher is confident in the certainty of the collected data, the classical Analytic Hierarchy Process (AHP) outperforms the Fuzzy Analytic Hierarchy Process (FAHP). However, in cases where the data is uncertain, the FAHP technique is advised [21]. The Fuzzy Analytic Hierarchy Process (FAHP) approach is used in this study to evaluate the weight criteria related to active tectonics based on morphometric parameters.

Artificial intelligence algorithms have attracted increasing attention in recent years as solid computational tools to simulate complicated phenomena in various academic domains [24,25,26,27]. Among the various methods of machine learning, artificial neural networks are considered the backbone of machine learning algorithms. This method’s key advantages are its learning capability based on the training process, which eliminates the requirement for statistical assumptions for the source data, and its ability to cope with nonlinear situations. Researchers in numerous scientific and engineering fields are interested in ANN models since they can correlate huge and complicated multi-parameter datasets without a prior understanding of the relationships between the parameters [28]. Applying machine learning in combination with GIS to analyze morphometric parameters allows us to understand the evolution of watersheds and provides a clear picture of landscape evolution [27].

A few investigations have used geomorphic indices of drainage basins by calculating their arithmetic mean or considering a weight for each index in the AHP to map out a relative tectonic activity (e.g., [15,27,29]), but no study employing artificial intelligence has been published on this subject as of yet. Considering that the morphometry of watersheds can reveal the relative tectonic activity of an area over a long period, this paper attempts to investigate a new method to combine the traditional FAHP technique with artificial intelligence algorithms as an innovation in morphotectonic evaluation. This research employs several geomorphic indices to assess the regional tectonic activity across such a broad area. The main goal of this study is to use machine learning algorithms and FAHP to identify which of the 423 drainage basins studied in the Makran subduction zone are more affected by tectonic features.

2. Materials and Methods

2.1. Study Area

The Makran subduction zone with east-west trending is located in southern Pakistan and Iran. This zone indicates a convergence zone where the Oman oceanic plate subducts beneath the Eurasian continental plate. The subduction probably started during the Late Cretaceous [30]. The western boundary of Makran is often called the Oman Line and Minab-Zendan Fault system, which runs northward and separates a highly seismic region in Zagros [31] from a region of low seismic activity in Makran to the east [32,33]. The eastern boundary of Makran is defined by a transfer zone consisting of three individual faults, namely the Ornach-Nal, Ghazaband, and Chaman faults [34,35].

In contrast to the majority of accretionary complexes in the world, the Makran accretionary wedge lacks an obvious trench [36]. The lack of a trench in this region may be because the angle of the subduction slab at the accretionary front is low, which may be caused by the presence of 7 km thick sediments with low compaction [37]. The seaward 70 km of the forearc comprises semi-consolidated and unconsolidated sediments with high pore fluid pressures and low seismic velocities. These sediments are capable of causing and failing large tsunamigenic slides [38].

The convergence rate decreases east to west along the Makran boundary [39]. In Makran, the average convergence rate is 4 cm/year, decreasing from 4.2 cm/year in the east to 3.65 cm/year in the west. Furthermore, GPS measurements suggest that the highest subduction rate of the Oman oceanic plate beneath the Eurasian plate occurs in the east at about 2.7 cm/year. In comparison, the lowest rate occurs in the west at about 1.95 cm/year [40]. Compared to other subduction zones, such as the Cascadia subduction zone at 35 mm/year [41], the Mexico subduction zone at 41 mm/year [42], the Sumatra subduction zone at 65 mm/year [43], the south Chile subduction zone at 70 mm/year [44], the Japan subduction zone at 80 mm/year [45], and the Tonga subduction zone at 160 mm/year [44], Makran is considered a relatively slow subduction zone.

Zarifi [46] states that the compressional stress direction along the Makran zone is rotating. The stress field in the eastern Makran is influenced by the collision between the Indian and Eurasian plates, while the western Makran stress field is affected by the Arabia-Eurasian collision. Eastern and western parts of the Makran subduction zone exhibit distinctive seismic behavior. The eastern region of Makran displays more seismic activity than the western region [38]. The plate boundary in western Makran lacks well-documented great instrumental and historical earthquakes. In contrast, eastern Makran was ruptured by thrust faulting during the 1945 earthquake and currently experiences earthquakes of varying magnitudes [38].

The main deformation phase from the Late Miocene to the Late Pliocene occurred by imbricating fans at the front of the Makran accretionary prism [47]. The seabed’s current geometry suggests that most of these imbricated fans remain active and that the accretionary prisms continue to be impacted by ongoing deformation from some branches of the imbricated fans. Following the Pliocene, the Makran coast and the mid-slope region experienced normal faulting, uplifting, and ductile flows [47].

The Makran subduction zone, like many other subduction zones around the world, has active mud volcanoes, some of which have formed along anticlinal axes [48]. As the largest mud volcanoes in the world [48], the mud volcanoes in Makran result from mud diapirism, providing evidence for the tectonic expulsion of mud and fluids seaward of the accretionary front. This confirms that tectonic forcing plays a significant role in forming mud volcanoes in Makran [49].

2.2. Extracting the Geomorphic Parameters

In this study, the morphotectonic zoning method was implemented in two steps. The watersheds of the research region were first determined using the Japan Aerospace Exploration Agency’s ALOS World 3D30m images (AW3D30), Ver.3.2 [50], with a resolution of 30 m (Figure 1). Numerous studies have demonstrated that AW3D30 DEM data can be compared to other open-source DEMs, such as SRTM and ASTER, with the advantage of being more accurate in delineating river basins and presenting drainage networks [51,52]. Next, ArcGIS Pro was used to evaluate the morphometric parameters and all physiographic indices to determine the watersheds with higher uplift. The eight indices employed in the morphotectonic process are shown in Figure 2.

2.3. Calculating Criterion Weights by FAHP

The Analytical Hierarchy Process (AHP) provides a systematic approach for evaluating complex decisions by combining qualitative and quantitative factors within a single framework, resulting in a prioritized list of alternatives [53]. Despite its widespread use, there have been concerns raised by researchers regarding certain limitations of the AHP, such as ambiguity in standardizing non-commensurate criteria (i.e., criteria that cannot be compared due to differences in size, type, or scale) and the influence of personal assessments [54]. These limitations can significantly impact the results of the AHP [55,56,57]. To address these issues, the Fuzzy Analytical Hierarchy Process (FAHP) has been developed to overcome the limitations and flaws of the AHP [58,59,60]. In this context, FAHP was employed to assess the weighting criteria related to morphotectonic zoning.

Mohebbi Tafreshi et al. [61] presented a fuzzy modelling method consisting of the following steps. The initial step in fuzzy models involves standardizing parameters using a fuzzy membership function. The experts’ opinions are utilized to assign fuzzy values to raw input values using a transformation function, where values close to 1 are deemed more suitable for the desired outcome, and values close to 0 are considered less suitable. The fuzzy logic extension of the ArcGIS Pro 2.8 software (version 10.8) offers various fuzzy membership functions. The selection of the fuzzification function is based on the nature, significance, and relation of each criterion to the goal. For this preliminary analysis, the Linear and Gaussian functions were chosen from the seven available fuzzy membership functions to standardize the factors.

The fuzzy linear transformation function applies a linear transformation between the minimum and maximum values specified by the user. Any value below the minimum will be assigned a 0 (definitely not a member), and any value above the maximum will be assigned a 1 (definitely a member). On the other hand, the fuzzy Gaussian function transforms the original values into a normal distribution [62].

The process of identifying the critical factors in the morphotectonic analysis was initiated by conducting a literature review, which resulted in identifying the most important indices, as shown in Figure 2. The relative significance of each parameter was then estimated through the use of the Analytic Hierarchy Process (AHP) and a constructed pairwise comparison matrix (8 × 8) based on the input factors, as determined by Saaty’s scale (

Table 1. Saaty’s 1–9 scale of relative importance [63].

Intensity of Importance	Interpretation
1	Equal importance
3	Moderate importance
5	Essential
7	Extreme importance
9	Extreme importance
2, 4, 6, 8	Intermediate values between adjacent scale values

). It is worth noting that the likelihood of inconsistencies in pairwise comparisons increases with the number of comparisons made, as stated in reference [63]. To account for this, AHP incorporates a consistency index (CI) to evaluate the calculated weight matrix. The weight is deemed acceptable if the CI is less than 10% (

Table 2. Pairwise comparison matrix for standardizing factor scores.

Linear and Areal Aspects	Hi	Bs	Cb	Er	Rn	Rs	Cc	Ff	Score
Hypsometric integral (Hi)		0.5	0.5	2	0.5	0.5	2	2	0.136
Basin shape (Bs)			2	2	2	2	2	2	0.037
Circularity basin (Cb)				0.33	0.5	2	2	0.5	0.123
Elongation ratio (Er)					2	2	2	2	0.084
Ruggedness number (Rn)						2	3	2	0.078
River sinuosity (Rs)							0.5	0.5	0.297
CoefficientCompactness (Cc)								0.5	0.050
Form factor (Ff)									0.197
CI	0.03

). Finally, the calculated weights were normalized to a scale between 0 and 1 to facilitate the integration of the weighted map layers.

2.4. Description and Application of the Criterion

Morphometric analysis of watersheds as a prerequisite of hydrological studies would be a valuable method to make informed management choices based on a more comprehensive view of the drainage network’s behavior and the morphology of the watershed.

In this study, eight linear and area-based indices of watersheds are evaluated. These indices provide a broad overview of the drainage basin network with regard to morphology and relief. These indices assess the stream network’s complexity, texture, and distortion due to neotectonic disturbances [64]. The extracted values are valuable in understanding the development of the drainage network concerning lithology and landscape. These indices are explained in detail as follows.

2.4.1. Hypsometric Integral (Hi)

The hypsometric curve represents the elevation distribution concerning the drainage area at various levels, including regional and continental scales [65,66]. The hypsometric integral (Hi) calculates the uneroded volume of a basin by determining the area under the hypsometric curve [67]. This can be computed through Equation (1).

$$\mathrm{Hi}=\frac{H_{avg-}{H}_{min}}{H_{max-}{H}_{min}}$$

(1)

A value of Hi greater than 0.6 signifies a tectonically active region with significant uplift and steep topography [68,69,70]. On the other hand, a mature drainage basin exhibits a moderate-to-low Hi value [71]. In this study, a fuzzy linear membership function was utilized, where a Hi value of 0.65 was assigned the highest weight, and the weight progressively decreased to zero as the Hi value approached 0.04.

2.4.2. Basin Shape (Bs)

The shape of drainage basins in active tectonic zones is generally elongated, but over time it tends to become more circular [67]. The Basin Shape Index (Bs) is calculated using Equation (2) [72,73] to describe this change.

$$\mathrm{Bs}=\frac{\mathrm{Bl}}{\mathrm{Bw}}$$

(2)

This index is based on the length (Bl) and width (Bw) of the basin, measured from the headwater to the mountain ridge and at the broadest point, respectively. These values are calculated using the minimum bounding geometry script. Basins that have lower Bs values are considered more circular in shape and are usually associated with low tectonic activity. On the other hand, steep basins that have high tectonic activity are elongated [29,74]. The most significant weight is assigned to Bs values around 7.9 using the increasing fuzzy linear membership, and as the Bs value decreases to 0.16, the weight decreases until it reaches zero.

2.4.3. Circularity Basin (Cb)

The concept of the circulatory basin (Cb) was introduced by Miller [75] and Strahler [76] and it is defined as the ratio of the area of the basin (A) to the area of a circle with the same perimeter as the basin Equation (3).

$$\mathrm{Cb}=\frac{4\mathsf{\pi }\mathrm{A}}{\mathrm{P}2}$$

(3)

where A represents the area of the basin and P is the basin’s perimeter. The circularity index ranges from 0 for a straight line to 1 for a perfect circle. The more circular the shape of the basin, the higher the value of Cb. Factors such as the length, frequency of streams (Fs), geological structures, land cover, climate, relief, and slope of the basin all play a role in determining the circulatory ratio, making it a key factor in determining the stage of a basin. Low, medium and high Cb values correspond to a basin’s youthful, mature, and ageing periods [77,78]. In this context, the decreasing fuzzy linear membership assigns the lowest weight to a Cb value close to 0.3, and as the Cb value decreases to 0.03, the weight also decreases until it reaches zero.

2.4.4. Elongation Ratio (Er)

The elongation ratio (Er) measures the shape of a drainage basin, representing the ratio of the diameter of a circle with the same area as the basin to the basin’s length [79]. A value of 1 indicates a perfectly circular shape, while a lower value suggests a more elongated and tectonically impacted basin (as described in Equation (4)).

$$\mathrm{Er}=\frac{2\sqrt{\frac{\mathrm{A}}{\mathsf{\pi }}}}{\mathrm{L}}$$

(4)

In this study, the highest and lowest Er values were found to be 0.17 and 0.06, respectively. Using the increasing fuzzy linear membership, the highest weight was assigned to an Er value close to 0.17, steadily decreasing as the Er value decreased to 0.06 and reached zero.

2.4.5. Ruggedness Number (Rn)

The Ruggedness Index, which combines two factors: relief (H) and drainage density (Dd), indicates the steepness and length of slopes. The calculation of the Ruggedness Index is carried out using Equation (5a–c).

$$\mathrm{Rn}=\mathrm{Dd} ·\mathrm{H}$$

(5a)

$$\mathrm{Dd}=\frac{\sum \mathrm{Li} }{\mathrm{A}}$$

(5b)

H = H_max – H_min

(5c)

Where Li is the length of the river, and A is the area of the watershed. H_max and H_min are the highest and lowest elevations of the watershed. The Rn close to 473 is allocated the highest weight utilizing the rising fuzzy linear membership. As the Rn value approaches 1, the weight drops until it reaches zero.

2.4.6. River Sinuosity (Rs)

The morphology of rivers that interact with an active fault zone tend to become uneven due to changes in slope and variations in curvature, leading to increased meandering [80]. Hence, the sinuosity of a river (Rs) can be a valuable indicator of tectonic activity in a drainage basin. A highly sinuous river suggests stability, while a straight river profile suggests ongoing tectonic activity. The formula in Equation (6) can be used to calculate the sinuosity of a river, where c represents the channel length, and v represents the straight length of the valley. Utilizing a decreasing fuzzy linear membership, the highest weight is assigned to the Rs value closest to 1.03, while the weight decreases as the Rs value approaches 2.58.

$$\mathrm{Rs}=\frac{C}{V}$$

(6)

2.4.7. Compactness Coefficient (Cc)

The compactness coefficient, also known as the Gravellus Index, evaluates the basin’s shape irregularity. The more irregular the shape, the higher the value of Cc [81]. This value can be a good indicator of the area’s tectonic activity level Equation (7).

$$\mathrm{Cc}=0.2821\frac{\mathrm{P}}{{\mathrm{A}}^{0.5}}$$

(7)

where P is the perimeter and A is the area of the watershed. With the use of increasing fuzzy linear membership, a Cc value close to 5.03 is given the lowest weight. As the Cc value decreases to 1.83, the weight decreases until it reaches zero.

2.4.8. Form Factor (Ff)

As Horton [11] defined, the form factor represents the relationship between a basin’s area and length, squared. The equation for calculating the form factor is shown in Equation (8).

$$\mathrm{Ff}=\frac{\mathrm{A}}{{\mathrm{L}}_{\mathrm{b}}^2}$$

(8)

Here, A is the basin area in square kilometers, and ${\mathrm{L}}_{\mathrm{b}}^2$ is the square of the basin length. A perfectly circular basin will have a form factor value lower than 0.78. Basins with a lower form factor are considered to be more elongated and potentially influenced by tectonic activity. Using the decreasing fuzzy linear membership, a form factor value close to 0.61 will be given the lowest weight. As the form factor decreases to 0.05, the weight will increase until it reaches one.

2.5. Machine Learning Algorithms

Traditionally, geological research has faced challenges with data sources that are complex and imprecise due to the size and complexity of geological objects [82]. However, with advancements in science and technology, new methods have emerged that enhance the precision of geological data and increase the amount of available information.

Artificial intelligence (AI) is a field within computer science that focuses on developing intelligent systems created by humans. It can be separated into two types: strong AI and weak AI. Weak AI views the creation of reasoning and problem-solving machines as impossible, while strong AI aims to build machines with the ability to think and make decisions. The study of AI encompasses various areas, such as expert systems, machine learning, natural language processing, computer vision, and recommendation systems. Machine learning explicitly explores how computer systems can improve automatically through experience and the fundamental laws that govern all learning systems, including humans, organizations, and computers.

Machine learning algorithms allow us to gain new insights and capabilities. The advent of deep learning has revolutionized the field of AI by addressing the challenges faced in traditional machine learning, such as limited model options, time-consuming training processes, and the complexity of determining model parameters. As a result, deep learning has become a pivotal area of AI advancement.

Representation learning enables machines to automatically identify patterns in raw data and learn the necessary representations for detection or classification. Deep learning methods are a type of representation learning that involve multiple levels of abstraction created by combining simple, nonlinear modules. As these modules are layered, increasingly complex functions can be learned and implemented.

The advancement of big data and AI has opened new opportunities in the field of geology in recent years. Thanks to increased computing power, particularly the progress in GPU technology, the limitations of big data and AI in terms of computing have been significantly reduced. This expands the potential applications and growth prospects of geology. AI is commonly utilized for geological surveys and resource exploration, such as mineral recognition and geochemical anomaly detection. This article delves into the topic of big data and AI in geology.

2.5.1. Artificial Neural Networks (ANNs)

Artificial Neural Networks (ANNs) are computational systems that draw inspiration from functioning biological nervous systems such as the human brain. ANNs comprise many interconnected computational nodes, called neurons, that work together to learn from inputs and optimize the output.

The input, typically a multidimensional vector, is fed into the input layer and distributed to the hidden layers. These layers make decisions based on the previous layer and assess the impact of random changes on the final output, a process referred to as learning. When multiple hidden layers are stacked, it is called deep learning. ANNs have several key features, including the ability to learn and adapt, generalize information, and process and analyze information uniformly and with error tolerance. These features make ANNs powerful tools for pattern recognition, classification, and nonlinear function estimation.

Activation functions are used in ANNs to transform the input into an output, which is then supplied as the input to the next layer. The weights of the connections between neurons in the layers are adjusted during the training phase, allowing the network to learn the patterns between inputs and outputs. ANNs can be divided into two learning models: supervised and unsupervised. In supervised learning, the network is trained with proper outputs for each input pattern. The weights are adjusted to minimize the error between the network output and the actual value. In unsupervised learning, the network discovers the relationship between the data patterns without needing actual responses.

ANNs can also be divided into two categories: recurrent and feed-forward networks. Recurrent networks have a feedback loop, while feed-forward networks do not. The neurons in each layer provide information to both the previous and the subsequent layers.

The ANNs offer these hypotheses:

• Data processing occurs in the units known as neurons. The neurons (or artificial neurons) present a model of brain neurons.
• The exchange of data is facilitated through communication between neurons.
• There is a weight for communicative ways between neurons.
• Every neuron utilizes a nonlinear function to process its inputs (weighted data), producing a specific output [82].

It is possible to identify a neural network via the communicative model between different layers of the network, the number of layers, the number of neurons, the neuron’s operational function, and the learning algorithm. However, no general principle is available regarding the standard size of the network components. It is an innovative approach in most cases where the multilayered networks have a different amount of neurons in each layer, and different activation functions and various learning rates do the training of these networks. Then, it is followed by a selection of the best network. The network training in the learning phase takes place via the adjustment of weights so that outputs can be predicted or classified based on a set of inputs [83].

In this paper, the input parameters of the neural network included Basin Shape, Circularity Basin, Coefficient Compactness, Elongation Basin, Form Factor, Hypsometric Integral, River Sinuosity, Drainage Density, and Number Ruggedness, and the FAHP outputs were considered as the network output parameters. The data on these parameters were divided into training, testing, and data validation. A total of 80% of these data were used for training, 10% of data for validation and the other 10% for testing.

2.5.2. Support Vector Regression (SVR)

Considering a data set with N elements ${\left\{\left(\mathrm{Xi}, \mathrm{yi}\right)\right\}}_{i=1}^N$, where X_i = [x_1,i, …, x_n,i] ∈ ${\mathrm{R}}^{\mathrm{n}}$ and X_i represents the i_th element in a space with n dimensions, y_i (y_i ∈ R) indicates the actual value for X_i, the definition of a nonlinear function is as follows: φ: ${\mathrm{R}}^{\mathrm{n}}$ → R^nh. For mapping the entry data, X_i represents an ${\mathrm{R}}^{\mathrm{nh}}$ space of high dimension known as feature space, which specifies the nonlinear transformation φ. Hence, a linear function $\mathrm{f}$ in a high-dimensional space, and consequently, the entry data, X_i can be related to output y_i. Equation (9) presents the linear function, i.e., SVR.

$$\mathrm{f}\left(\mathrm{X}\right)={\mathrm{W}}^{\mathrm{T}}·\mathsf{\phi }\left(\mathrm{X}\right)+\mathrm{b}$$

(9)

where $\mathrm{b}$ ∈ $\mathrm{R}$ and $\mathrm{W}$ ∈ ${\mathrm{R}}^{\mathrm{n}}$, and $\mathrm{f}\left(\mathrm{X}\right)$ is the foretold value. As indicated in Equation (10), the empirical risk is minimized by the SVR.

$${\mathrm{R}}_{\mathrm{reg}}\left(\mathrm{f}\right)=\mathrm{C}{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{N}}{\mathsf{\Theta }}_{\mathsf{\epsilon }}\left({\mathrm{y}}_{\mathrm{i}}-\mathrm{f}\left({\mathrm{X}}_{\mathrm{i}}\right)\right)+\frac{1}{2}||{\mathrm{W}}^{\mathrm{T}}||$$

(10)

where the cost function is represented by Θε(yi − f(Xi)), regarding the ε-SVR, as shown in Equation (11), a loss function ε-insensitive is utilized:

$${\mathsf{\Theta }}_{\mathsf{\epsilon }}\left(\mathrm{y}-\mathrm{f}\left(\mathrm{X}\right)\right)=\left\{\begin{array}{ll}\left|\mathrm{y}-\mathrm{f}\left(\mathrm{X}\right)\right|-\mathsf{\epsilon }& \mathrm{I}\mathrm{f}\left|\mathrm{y}-\mathrm{f}\left(\mathrm{X}\right)\right|\ge \mathsf{\epsilon } \\ 0& \mathrm{I}\mathrm{n} \mathrm{a}\mathrm{n}\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r} \mathrm{c}\mathrm{a}\mathrm{s}\mathrm{e}\end{array}\right.$$

(11)

The nonlinear function φ is determined using Θ_ε in the R^nh space for finding a function with the ability to fit present training data with a deviation equal to or below ε. Using the mentioned function, the training error is minimized between the data training, and Equation (12) provides the function ε-insensitive [84].

$$\underset{\mathrm{W},\mathrm{b},{\mathsf{\xi }}^{*},\mathsf{\xi }}{\min }{\mathrm{R}}_{\mathrm{reg}}(\mathrm{W},{\mathsf{\xi }}^{*},\mathsf{\xi })=\frac{1}{2}{\mathrm{W}}^{\mathrm{T}}\mathrm{W}+\mathrm{C}{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{N}}({\mathsf{\xi }}_{\mathrm{i}}^{*}+{\mathsf{\xi }}_{\mathrm{i}})$$

(12)

The training errors of f(X) and Y are punished by Equation (12) via the function ε-insensitive. Using the parameter C, the compromise between the points meeting condition |f(X) − y| ≥ ε in Equation (11) and the model complexity (vector W) is determined. With C → ∞, a small model margin is observed that is adjusted to the data. When C → 0, there is a large model margin, which is why it softens. Lastly, ξ_i represents errors more minor than −ε and ξ_i^∗ indicates training errors larger than ε.

For solving the regression problem, the internal product of Equation (9) can be replaced by kernel K() functions. Thus, this operation can be performed in a higher dimension by low-dimensional space data input regardless of knowledge of the transformation φ [85].

2.5.3. Multivariate Linear Regression (MLR)

The regression method is applied to two theories. Firstly, regression analysis is typically utilized for prediction and forecasting, and its application significantly overlaps the machine learning field. Secondly, it is possible to use regression analysis sometimes for determining causal relationships between the dependent and independent variables. Regression alone presents just relationships between a fixed dataset of different variables and a dependent variable.

Based on the regression models, the dependent variables are predicted by the independent variables. The value of the dependent ‘y’ variable is estimated by regression analysis because of the range of independent variable ‘x’ values. This articlediscusses polynomial and linear regression, which fit better into the predictive model. Regression could be multiple regression or a simple linear regression [86].

Simple linear regression is a statistical method used to model the relationship between a dependent variable and a single independent variable. An equation represents it, y = β0 + β1x + ε, where y is the dependent variable, x is the independent variable, β0 and β1 are coefficients, and ε represents the error term. The goal of simple linear regression is to determine the strength and direction of the relationship between the dependent and independent variables and estimate the effect of the independent variable on the dependent variable [86].

Multivariate linear regression is a statistical method for predicting the result of an answer variable, which uses some explanatory variables, as shown in Equation (13). Multiple Linear Regression (MLR) aims to establish a linear association between the dependent variable y and one or more independent variables x, which will then be analyzed.

$$\hat{\mathsf{\beta }}={\left({\mathrm{X}}^{\mathrm{T}}\mathrm{X}\right)}^{-1}{\mathrm{X}}^{\mathrm{T}}\mathrm{y}$$

(13)

2.6. Integrating the FAHP and ML Algorithms

We first used the FAHP model to generate a target database suitable for our ML. This database served as the basis for training and testing the algorithms. Next, we randomly selected 212 out of 423 watersheds across our study area for training and testing the algorithms. By randomly selecting the watersheds, we ensured a representative sample of the study area. Here, raw data from eight different morphometric indices were used to train and test the machine learning algorithms, as shown in Figure 2. These indices were selected based on their relevance and importance in the morphotectonic characterization of watersheds. We then used these raw data as input features for the machine learning algorithms to generate predictions for the target database. Using a representative sample of the study area and including relevant input features, we aimed to create a distinct training network and improve the accuracy and precision of our predictions and our cost-effectiveness. Three distinct ML algorithms, namely, ANNs, SVR, and MLR methods, were employed for each set. It is important to note that the training and testing set includes 212 watershed indices. After optimizing and obtaining the best algorithms using statistical equations, the method was generalized for the whole domain.

To assess the effectiveness of the ANNs, SVR, and MLR techniques, Mean Squared Error (MSE), Mean Absolute Error (MAE), and Correlation Coefficient (R²) were utilized as performance metrics. These metrics are described in Equation (14a–c) [82].

$$\mathrm{MSE}=\frac{1}{\mathrm{N}}{\sum }_{\mathrm{i}=1}^{\mathrm{N}}{\left({\mathrm{y}}_{\mathrm{i}}-{\hat{\mathrm{y}}}_{\mathrm{i}}\right)}^2$$

(14a)

$$\mathrm{MAE}=\frac{1}{\mathrm{N}}{\sum }_{\mathrm{i}=1}^{\mathrm{N}}\left|{\mathrm{y}}_{\mathrm{i}}-{\hat{\mathrm{y}} }_{\mathrm{i}}\right|$$

(14b)

$${\mathrm{R}}^2=1-\frac{{\sum }_{\mathrm{i}}{\left({\mathrm{y}}_{\mathrm{i}}-{ \hat{\mathrm{y}} }_{\mathrm{i}}\right)}^2}{{\sum }_{\mathrm{i}}{\left({\mathrm{y}}_{\mathrm{i}}- \overline{\mathrm{y}} \right)}^2}$$

(14c)

We used ANN for regression here and utilized the supervised method. Also, we use L2 (ridge) regularization to avoid overfitting. To have a faster converging, we need zero-centered data points, and each dimension should be scaled according to its standard deviation. Thus we normalize our data.

Figure 3 shows the performance of the ANN method for the training and validation sets. From epoch 200 and on, the trend was deemed acceptable. After completing the learning process, the ANN model demonstrated high accuracy with an R² value of 0.974, MSE of 4.14 × 10⁻⁶, and MAE of 0.00151.

After developing and training the SVR model with the value of 0.947, MSE and MAE obtained 8.12 × 10⁻⁶ and 1.94 × 10⁻³, respectively.

After developing and training the MLR model with the value of 0.967, MSE and MAE obtained 5.06 × 10⁻⁶ and 1.61 × 10⁻³, respectively.

According to the simulation results obtained, there is a good agreement between the results obtained from ANN, SVR and MLR models but ANN model outperformed the other models. All three algorithms demonstrated similar error patterns in our study, with samples 16 and 96 displaying particularly large prediction errors. We have further investigated the reason for these high prediction errors and found that the location of these two samples may have contributed to the discrepancies. Specifically, sample 16 is located in the northernmost part of the study area, where geological and climatic conditions differ significantly from the rest of the watersheds. Similarly, sample 96 is situated in the southernmost part of the study area, which also has distinct geological and climatic conditions compared to the majority of the watersheds.

These differences in geological and climatic conditions may have played a role in the higher prediction errors observed for samples 16 and 96. Our findings suggest that when analyzing and predicting data in areas with distinct geological and climatic conditions, particular attention should be paid to samples that are situated in the outlier regions of the study area. By doing so, we may better understand and account for the differences in geological and climatic conditions and obtain more accurate predictions.

Overall, our study highlights the importance of considering the spatial distribution of data and the potential impact of varying geological and climatic conditions when making predictions. Further research may be necessary to determine how these findings can be applied in other contexts or how they can inform the development of more accurate predictive models.

3. Results and Discussion

ML techniques use algorithms to learn from data and make predictions or decisions without being explicitly programmed. Meanwhile, the fuzzy technique assigns degrees of truth to statements or rules, making it useful in assessing morphometry in watersheds where data such as rainfall patterns or soil properties can be difficult to quantify accurately. These techniques are useful in engineering problems where the relationship between input and output variables is complex and challenging to model analytically, such as landslides and mass movements [87,88].

However, both techniques have limitations in practical applications. Fuzzy logic relies heavily on expert knowledge to define rules and membership functions, resulting in models that are challenging to interpret and validate. ML techniques are limited by the quality and quantity of available data and the choice of algorithm and parameters used. Overfitting or underfitting data can result in poor predictive performance on new data. Successful application requires careful consideration of their limitations and appropriate use of available data and expertise.

Our study divided the dataset into three groups: a training dataset that accounted for 80% of the collected data, a validation dataset that accounted for 10%, and a testing dataset that comprised the remaining 10%. The performance of the machine learning methods on both the training and testing data is depicted in Figure 4, Figure 5 and Figure 6. These figures compared the target (FAHP output) and output (predicted by algorithms) values and calculated the model error. Numerical results closely matched those obtained from the ANN, SVR, and MLR methods (

Table 3. Performance of the machine learning algorithms.

Methods	MSE	MAE	R2
ANN	4.14 × 10−6	0.00151	0.974
SVR	8.12 × 10−6	0.00194	0.947
MLR	5.06 × 10−6	0.00161	0.967

). However, the ANN model showed superior performance compared to the other methods.

We investigated the potential of coupled artificial intelligence algorithms and FAHP to predict the watershed’s behavior in response to the region’s tectonics. After selecting the best from among the ANN, SVR and MLR methods, the model was applied to whole watersheds in the Makran Subduction zone (Figure 7).

Mapping the results acquired by artificial intelligence algorithms for analog interpretation and comparing them with tectonic features and seismic events indicates the acceptable accuracy of the algorithms used. As shown in Figure 7, the findings measured for each watershed are illustrated in the range of low- to high levels of tectonic activity, which correlates well with the tectonic situation. On the west side of the figure, the location of the Minab fault zone is in good agreement with the area marked on the map as a zone with relatively high tectonic activity. Additionally, this region experiences the effects of two unique geological occurrences: the convergence of the Iranian and Arabian tectonic plates in the west and the subduction of the Oman sea plate beneath the Iranian continent in the east, both ongoing processes. On the eastern side of the region, it is observed that tectonic activity is estimated to be relatively high based on morphometric parameters.

This area aligns with the Chaman fault zone in Pakistan and marks the eastern boundary where the Oman oceanic plate is subducting beneath the Iranian continental plate. Moreover, instrumental epicenters also confirm the higher activity of this zone. The noteworthy point in this figure is the presence of instrumental epicenters in the central region, which are consistent with the results obtained using AI approaches, albeit with a slight shift to the left.

It is affected by the inclination of seismic faults, earthquakes, the geometry of the subducted oceanic plate, and instrumental error. The development of morphometric indicators can be affected by various factors, including geological heterogeneity, land use, and climate conditions. While high erosion rates in certain geological layers may contribute to the lack of clear development of tectonic indicators in some watersheds, it is important to note that this is only one potential explanation among many. Other factors, such as high levels of sedimentation, variations in precipitation patterns, and land use changes, could also contribute to the observed variability in technical indicator development. Additionally, interactions between multiple factors may contribute to the observed patterns.

Furthermore, it should be noted that seismic instrument data belong to only the last hundred years, whereas the morphometric features date back several million years which makes it more reliable. To better understand the underlying causes of the variability, it may be necessary to conduct further analyses that consider these different factors and potential interactions. It is also important to consider the limitations of the data and methods used in the study, as these may affect the accuracy and reliability of the technical indicators.

4. Conclusions

In conclusion, the morphometric analysis of watersheds in a region has the potential to reveal insights into long-term tectonic activity, and the present study aims to develop a novel approach to enhance this analysis. The proposed method leverages advanced machine learning techniques to reduce computational cost and time, especially in the case of large-scale regions. In this study, we examined the morphometry of watersheds in the Makran subduction zone using three artificial intelligence methods: ANN, SVR, and MLR. The results of our analysis indicate that the ANN method is the most accurate of the three, with a value of 0.974. The MSE and MAE values of the ANN method were also found to be lower compared to the other methods. Based on these findings, it can be concluded that the ANN approach can be effectively utilized in the morphometric analysis of watersheds and provides better results than previous techniques.