Comparison of Liquefaction Potential Prediction Results Using Machine Learning Methods Based on CPT Data

1. Introduction

Soil liquefaction is an earthquake phenomenon in which saturated granular soils lose their strength and stiffness under seismic loading and behave more like a viscous fluid than a solid medium. This transformation, recognized by the 1964 Niigata and Great Alaska events, has been observed and reported repeatedly in many earthquakes for about half a century. After the awareness of liquefaction, a significant portion of earthquake damage has been attributed to soil liquefaction, and many studies and procedures have been included in the literature to estimate liquefaction potential and prevent liquefaction hazards [1,2,3,4].

At the core of liquefaction lies the mechanism of cyclic shearing under seismic forces. As strong ground shaking occurs, cyclic loading causes a progressive buildup of pore water pressure in the saturated soils. When pore pressure nearly equals the initial confining stress, the effective stress-and hence the soil’s shear strength–drops dramatically. While loose sands often experience complete loss of strength and significant deformations, medium-density sands, silty sands, and sandy silts may not fail entirely but still undergo considerable softening and cyclic strain [5,6,7].

Traditional liquefaction potential assessments are semi-empirical methods based on both field tests and past liquefaction observations due to the difficulty of taking undisturbed samples from frequently liquefaction-occurred loose sandy soils. These methods rely on key parameters such as soil properties (e.g., relative density and gradation), environmental conditions (e.g., groundwater depth), and earthquake (e.g., peak ground acceleration) to classify a site as susceptible or resistant to liquefaction. One of the seminal approaches in this area is based on the standard penetration test (SPT) or cone penetration test (CPT) combined with the computation of the cyclic stress ratio (CSR), which utilizes parameters like peak ground acceleration to quantify the potential for liquefaction [6,8,9].

Despite their widespread use and proven track record, traditional methods face inherent challenges. In particular, the influence of various in situ state parameters-such as fines content, plasticity index, cementation, and overconsolidation-on $q_c$ and $V_s$ is complex and not fully captured by deterministic approaches. For example, while $q_c$ may decrease with increasing fines content and increase with cementation and overconsolidation, $V_s$ tends to behave inversely when fines content exceeds a critical threshold. Such complexities underscore the limitations of conventional methods in accurately modeling the multifaceted nature of soil behavior under seismic loading [9].

Although these empirical techniques have formed the backbone of liquefaction hazard assessment for decades, in recent years, machine learning (ML) methods have emerged as transformative tools in civil and geotechnical engineering, offering innovative solutions to complex problems that traditional techniques often struggle to address. One of the key advantages of machine learning techniques over conventional empirical methods is their ability to model nonlinear and multidimensional interactions without the need for a predetermined mathematical form between input and output variables. Traditional approaches typically rely on simplified assumptions and empirical correlations that may overlook the intricate interplay among various soil properties, environmental conditions, and earthquake characteristics. In contrast, ML algorithms can flexibly adapt to the complex behavior of soils under seismic loading, capturing subtle variations that influence liquefaction potential and thereby enhancing prediction accuracy [10] ML methods can automatically learn from large datasets-such as those derived from standard penetration tests (SPT), cone penetration tests (CPT), and shear wave velocity ($V_s$) measurements-by adjusting internal parameters to detect underlying patterns and relationships in the data [5].

Moreover, the adaptive learning process inherent in ML techniques offers a dynamic framework that evolves with emerging research and changing site conditions. This continuous updating capacity makes ML-based approaches particularly well-suited for capturing the transient and variable nature of soil behavior under seismic loading, ultimately leading to more resilient and cost-effective engineering designs. Consequently, machine learning is paving the way for a new era of innovation in geotechnical engineering by overcoming traditional limitations through advanced data analytics and computational intelligence [10,11].

2. Materials and Methods

2.1. Dataset

Within the scope of the study, the dataset is prepared for loose sandy soils that is known as susceptible to liquefaction, soils such as gravel, clay, unsaturated soils are excluded. Dataset used in the analyses consists of both liquefied and non-liquefied cases due to earthquakes [12,13,14]. The selected sources, which reported observed liquefaction and non-liquefaction cases for the 22 earthquake events, are listed in Table 1. The number of cases reported as liquefied and non-liquefied for each earthquake event is given in Figure 1. The total number of data points is 480, and proportions of observed liquefaction and non-liquefaction cases are 66% (316 cases) and 34% (164 cases), respectively, as shown in Figure 2.

Each data point in the created dataset contains nine quantities, which are earthquake moment magnitude ($M_w$), peak ground acceleration ($a_{max}$), measured depth of $q_c$ and $f_s$ ($d$), cone tip resistance ($q_c$), cone friction resistance ($f_s$), cyclic stress ratio (CSR), soil behavior type index ($I_c$), total stress (${\sigma }_v$) and effective stress (${\sigma }_v^{\prime }$),. $M_w$, $a_{max}$, $d$, $q_c$, $f_s$, ${\sigma }_v$, ${\sigma }_v^{\prime }$ quantities are determined by measurement, and CSR, $I_c$ quantities are determined depending on the measured quantities.

Researchers conduct studies in cases where liquefaction is observed as a result of earthquakes, and they provide valuable contributions to the literature by reporting their results. Some of these studies added to the literature may include laboratory study results such as dynamic triaxial, monotonic triaxial, consistency limit testing, and sieve analyses. Even simple tests such as sieve analysis and consistency limits require laboratory work. Studies based on field data results can be completed faster and more economically. Therefore, some studies in the literature may not contain laboratory results. Thus, approximately 40% of the studies evaluated within the scope of this study do not have fine material data [6,9,10,11,14,15,16]. Field tests are performed as a standard practice in studies about cases of liquefaction after an earthquake. Therefore, almost all of the case studies evaluated within the scope of this study have field test data. The number of data points is one of the major factors that governs the accuracy and reliability of the machine learning methods [17]. For this reason, the high number of input data in machine learning methods was prioritized, in this study. For this purpose, it has been preferred to use the soil type ($I_c$) parameter instead of fine material content (FC) parameter that may be missing in some studies. $I_c$ is a parameter that indicates the soil type (gravel, sand, silt, silty sand, etc.), and it can be determined based on CPT test results for all data points [3].

The following parameters earthquake moment magnitude ($M_w$), peak ground acceleration ($a_{max}$), depth ($d$), cone tip resistance ($q_c$), cone friction resistance ($f_s$), total stress (${\sigma }_v$), effective stress (${\sigma }_v^{\prime }$), cyclic stress ratio (CSR), soil behavior type index ($I_c$) were evaluated as input parameters of the prediction models. The liquefaction status in the database, which is evaluated as the output parameter, is considered as a binary classification problem to separating patterns between samples where liquefaction occurs (positive) and those where liquefaction does not occur (negative).The histogram and scatter plots of the input parameters used in the study are given in Figure 3.

In machine learning models, variables with different scales can negatively impact model performance. Specifically, some algorithms may be biased toward variables with larger values, causing the model to become unbalanced. To prevent this issue, a normalization process is applied during the data preprocessing stage. In this study, the min-max normalization method given in Equation 1 was used to scale the variables to a certain range. Min-max normalization ensures that all values fall between 0 and 1, allowing variables of different magnitudes to contribute equally.

$$X_{normalized}={\displaystyle \frac{X-X_{min}}{X_{max}-X_{min}}}$$

(1)

Here $X$ represents the original value, $X_{min}$ represents the lowest value in the dataset, and $X_{max}$ represents the highest value [18,19,20,21].

As part of this study, min-max normalization was applied to all independent variables used for liquefaction prediction. This process helped balance the effect of variables with different scales, increased the learning efficiency of the model, and minimized errors due to differences in scale between variables [7,17,22]. Traditional liquefaction potential prediction method [3].

Liquefaction potential has been computed for several decades by evaluating field tests such as SPT, CPT and shear wave velocity tests [23]. Cyclic stress ratio (CSR) is calculated using Equation 2, was proposed by Seed and Idriss [24], and it is a worldwide accepted and frequently preferred index for liquefaction computations. The cyclic resistance ratio (CRR) can be computed is calculated using the CPT data [3].

$$CSR=0.65{\displaystyle \frac{a_{max}}{g}}\left({\displaystyle \frac{{\sigma }_{v0}}{{\sigma }_{v0}^{\prime }}}\right)r_d$$

(2)

Here $a_{max}$ is maximum horizontal surface ground acceleration, $g$ is ground acceleration, ${\sigma }_{v0}$ is total vertical stress, ${\sigma }_{v0}^{\prime }$ is effective vertical stress and $r_d$ is stress reduction coefficient.

The Stress Reduction Coefficient $r_d$, which is included in Equation 2 and allows the flexibility of the soil profile to be considered, is calculated by Equation 3. The z value in the equation indicates the depth (m).

$$r_d={\displaystyle \frac{(1-0.4113z^{0.5}+0.04052z+0.001753z^{1.5})}{(1-0.4177z^{0.5}+0.05729z-0.006205z^{1.5}+0.001210z^2)}}$$

(3)

The curve separating the areas where liquefaction is observed/not observed, shown in Figure 4, which was presented by combining the studies in the literature to be used for estimating CRR using CPT tip resistance, is approximately given by Equation 10. This graph is known to provide more valid results for clean sands [3], so the values used in the calculations need to be converted to values suitable for clean sands.

The soil behavior table given in Figure 5 was developed using the appropriate value of $n=1.0$ for clayey soil types. However, a value of $n=0.5$ would be appropriate for clean sands, and a value between 0.5 and 1.0 for silts and sandy silts [3]. The soil behavior type boundaries shown in Figure 5 were determined according to the $I_c$ value ranges given in Table 2. The soil behavior type index $I_c$ is not valid for zones 1, 8 or 9. Throughout the normally consolidated zone, the soil behavior type index increases in direct proportion to the fines content and plasticity. This relationship is also shown in Figure 6.

The Cyclic Stress Ratio (CRR) is calculated by following the process steps explained below and summarized in the flow chart given in Figure 8 [3].

Firstly, the soil types defined as clay are distinguished from sand and silt. This distinction is calculated by assuming $n=1.0$ in the equation given by Equation 4 and calculating the Q value.

$$Q=\left[\left(q_c-{\sigma }_{vo}\right)/P_a\right]\left[{\left(P_a/{\sigma }_{vo}^{\prime }\right)}^n\right]$$

(4)

Here, $Q$ is the dimensionless CPT end resistance, $q_c$ is the CPT end resistance, ${\sigma }_{vo}$ is the vertical total stress, $P_a$ is the air pressure (100 kPa = 1 atm), ${\sigma }_{vo}^{\prime }$ is the vertical effective stress, and $n$ is a value ranging from 0.5 to 1.0 representing the soil type.

$$F=\left[f_s/\left(q_c-{\sigma }_{vo}\right)\right]\times 100\mathrm{\%}$$

(5)

The soil behavior type index $I_c$ given by Equation 6 is calculated using the $Q$ value calculated by Equation 4 and the $F$ value calculated by Equation 5 . Here, $F$ is the dimensionless friction ratio and $f_s$ is the CPT friction resistance.

$$I_c={\left[{\left(3.47-\log Q\right)}^2+{\left(1.22+\log F\right)}^2\right]}^{0.5}$$

(6)

If the calculated $I_c$ value is greater than 2.6, the soil is classified as clayey, and the analysis is completed by assuming that it is too rich in clay to liquefy. However, soil samples should be taken and tested with laboratory tests to verify the soil type and liquefaction resistance. If the calculated $I_c$ value is less than 2.6, the soil is classified as granular and therefore, $Q$ and $I_c$ values are recalculated by accepting $n=0.5$ in the equation given by Equation 4.

If the recalculated $I_c$ value is less than 2.6, the soil is classified as non-plastic and granular. This $I_c$ value is used to estimate the liquefaction resistance in the following steps. If the recalculated $I_c$ value is greater than 2.6, the soil is very silty and probably plastic. In this case, the normalized CPT tip resistance $q_{c1N}$ value for silty sands given by Equation 7 should be recalculated by taking $n=0.7$.

$$q_{c1N}=\left({\displaystyle \frac{q_c}{P_a}}\right){\left({\displaystyle \frac{P_a}{{\sigma }_{vo}^{\prime }}}\right)}^n$$

(7)

The $I_c$ value is then recalculated using the recalculated $q_{c1N}$ value instead of the $Q$ value in Equation 6. This $I_c$ value is then used to calculate the liquefaction resistance.

For silty sands, the normalized CPT tip resistance $q_{c1N}$ is corrected to the equivalent clean sand value ${\left(q_{c1N}\right)}_{cs}$ using Equation 8.

$${\left(q_{c1N}\right)}_{cs}=K_c\times q_{c1N}$$

(8)

Here $K_c$ is the correction factor for the grain characteristics calculated by Equation 9.

$$\begin{gathered} if{I}_c\le 1.64,K_c=1.0 \\ if{I}_c>1.64,K_c=-0.403I_c^4+5.581I_c^3-21.630I_c^2+33.750I_c-17.880 \end{gathered}$$

(9)

The $K_c$ curve given by Equation 9 is shown in Figure 7. For $I_c>2.6$, the curve shown as a dashed line indicates that the soils in this $I_c$ range must have a high clay content or be plastic to liquefy.

The obtained ${\left(q_{c1N}\right)}_{cs}$ value is used in the relationship given by Equation 10 and as a result of calculating the ${CRR}_{7.5}$ value, the cyclic resistance ratio (CRR) required for the liquefaction of the soil and the cyclic stress ratio (CSR) created by the earthquake are compared with each other and the liquefaction potential can be estimated with the ratio obtained. As a result of these calculations, it is concluded that there is a risk of liquefaction in the layer if the FS safety factor is less than 1.0. Here ${\left(q_{c1N}\right)}_{cs}$ is the clean sand CPT tip resistance normalized to 100 kPa (1 atm).

$$\begin{gathered} {if\left(q_{c1N}\right)}_{cs}<50,{CRR}_{7.5}=0.833\left[{\left(q_{c1N}\right)}_{cs}/1000\right]+0.05 \\ if50\le {\left(q_{c1N}\right)}_{cs}<160,{CRR}_{7.5}=93{\left[{\left(q_{c1N}\right)}_{cs}/1000\right]}^3+0.08 \end{gathered}$$

(10)

2.2. Machine Learning Methods

Liquefaction, defined as the loss of soil strength during earthquakes, poses a significant risk to engineering structures. While traditional empirical and experimental methods can be used for liquefaction prediction, machine learning methods have recently emerged as a powerful alternative in this field.

Machine learning enables the development of data-driven models that can learn patterns from large datasets to make predictions. Various machine learning algorithms, such as logistic regression, decision trees, support vector machines (SVM), random forest, and artificial neural networks (ANN), have been utilized for liquefaction prediction.

Machine learning-based approaches can predict liquefaction potential by utilizing different variables, including soil properties (e.g., grain size, void ratio), groundwater level, and seismic parameters (e.g., acceleration, magnitude). The advantages of these methods include:

More flexible and faster than traditional methods.

Capable of making high-accuracy predictions with large datasets.

Easily updatable with new data.

In this study, different models were created using Random Forest (RF) and Artificial Neural Networks (ANN) methods for liquefaction prediction and their performances were evaluated and compared with the traditional method. The data set used in the modeling includes field observations with and without liquefaction, and the algorithms are aimed to learn these patterns and make accurate predictions. It is thought that machine learning-based predictions can be a useful tool in engineering decision-making processes.

2.2.1. Random Forest (RF)

The Random Forest method was developed by Leo Breiman in 2001 [25]. While developing the Random Forest method, Breiman was inspired by several previous studies and combined the ideas in these studies to create a new approach. As Breiman stated in his article [25], this method is based on the combination of two main concepts: Bagging (Bootstrap Aggregating) [26] and Random Subspace Method [27].

Random Forest is a machine learning algorithm used for both classification and regression problems that works by combining multiple decision trees. It is based on the basic principle that individual trees may be weak, but when they work together, they can form a strong predictive model.

One of the significant differences between decision trees and random forest methods is the management and reduction of overfitting problem. Overfitting problem is the situation where the model makes predictions that are very close to the learning data, in other words, it memorizes the problem. In this case, the model can make very high accuracy predictions with the learning data, but it can make very low accuracy predictions with the test data. While this problem is evident in the decision tree method, it is controlled in the Random Forest method with the techniques of “Bagging” [26] and “Random Subspace” [27]. These techniques reduce the variance of random forests, increase the generalization ability of the model, and reduce the risk of overfitting. In this way, the random forest method generally provides the creation of high-performance and reliable models.

2.2.1.1. Bagging (Bootstrap Aggregating) [26]

Bagging is an ensemble learning technique used particularly in statistical learning. (Figure 9)

The main purpose of Bagging is to increase generalization power and reduce overfitting by combining many independent models. In this context, Bagging consists of Bootstrap Sampling, Training of Independent Models and Aggregating Model Results steps.

Bootstrap Sampling is the process of creating a new training data set by taking random samples from the original training data. The feature of this sampling is that each sample is selected by replacing the place in the training set. Therefore, the same data point can be selected more than once.

Training of Independent Models is the process of creating a model that is trained independently for each Bootstrap Sampling and each learns on a different data set.

Combining Model Results is the process of combining the predictions of each model after the training phase is completed. For classification problems, the class that receives the most votes from the class predictions made by each model is determined, and for regression problems, the final prediction is made by calculating the arithmetic average of the predictions of all models.

2.2.1.2. Random Subspace Method [27]

The Random Subspace Method is also used to reduce over-learning and increase model generalization power, like the Begging Technique. This method ensures that each model is trained by selecting a random subset from all available features (input parameters). The main differences between the two methods are shown in Table 3.

2.2.1.3. Structure of Random Forest Method

Instead of branching each node using the best branch among all variables, Random Forest branches each node using the best branch among randomly selected variables at each node (Figure 10). Each dataset is generated from the original dataset with replacement. Then, trees are developed using random feature (input parameter) selection. In simple terms, Random Forest creates multiple decision trees and combines them to obtain a more accurate and stable prediction.

Decision trees consist of three main structures: root, decision node and leaf (Figure 11).

Root is the top node of the decision tree and is where the analysis of the entire data set begins. Each branch of the decision tree is assumed to start from this root.

Decision Node is the node where a decision rule is applied and branches. Decision nodes divide the data set by means of a feature and a threshold value determined for this feature. (For example, a decision node can check whether a feature value is greater or less than a certain threshold value. This divides the data set into two subsets.)

Leafs are the last nodes of the decision tree and contain the results. A class label (for classification problems) or a prediction value (for regression problems) is found in the leaf nodes. Leaves are where a branch ends, and the results are obtained. In other words, a leaf node is reached at the end of each branch.

2.2.1.4. Splitting Criteria

Splitting criteria are a crucial step in the construction of decision trees. These crises decide which feature and which equality persists at a node. The goal is to maximize the information gain at each node by splitting the dataset into purer (homogeneous) subgroups. Different categories can apply different splitting criteria. The most used splitting criteria are:

Gini impurity is a metric that measures the mixing of classes in a node and is calculated by Equation 11. Here $C$ is the number of classes and $p_i$ is the proportion of observations belonging to class $i$. Gini impurity is widely used by decision trees (e.g. CART - Classification and Regression Trees). When a node is pure (i.e. contains a single class), its Gini value is 0.

$$Gini=1-\sum _{i=1}^C{p}_i^2$$

(11)

Information gain measures the difference in entropy before and after the split. Entropy is used to measure the confusion or uncertainty in the data set (Equation 12 and Equation 13). Here $D$ is the dataset, $D_k$ are the subsets after the split. Information gain measures the difference in entropy before and after the split. Entropy is used to measure the confusion or uncertainty in the dataset.

$$Entropy=-\sum _{i=1}^C{p}_i·{\log }_2\left(p_i\right)$$

(12)

$$Information Gain=Entropy\left(D\right)-\sum _k{\displaystyle \frac{\left|D_k\right|}{\left|D\right|}}·Entropy\left(D_k\right)$$

(13)

The information gain can be biased for features (input parameters) with many classes. For example, if a feature has a different value for each observation, the information gain can be maximum. The information gain ratio normalizes the information gain by dividing it by the partition entropy to correct for such biases (Equation 14). This method is widely used in the C4.5 algorithm.

$$Information Gain Ratio={\displaystyle \frac{Information Gain}{Entropy\left(split\right)}}$$

(14)

In this study, since the Random Forest method, which is a decision tree method, was preferred, the Gini Impurity method was used as the division criterion applied at the nodes.

3. Results

3.1. Hyperparameter Optimization

Hyperparameters are predefined settings that directly affect the learning process and overall accuracy of the model. Wrongly selected hyperparameters can lead to problems such as overfitting or underfitting, causing the model to perform poorly on real-world data. Therefore, determining the optimal hyperparameters is a critical step to increase the generalization ability of the model and maximize its accuracy.

In this study, hyperparameter optimization was applied to improve the performance of RF and ANN models. The optimized hyperparameters used in the RF model are shown in Table 4, and the optimized hyperparameters used in the ANN model are shown in Table 5. This optimization process increased the generalization ability of the models and provided more reliable results in liquefaction prediction.

3.2. Correct and Incorrect Estimated Values

True Positive (TP), False Positive (FP), True Negative (TN) and False Negative (FN) concepts are the basic parameters used when evaluating the prediction results of a model in classification problems. These concepts, summarized in Table 6. , are calculated because of comparing the results predicted by the model with the actual results.

Classification models have four key outcomes. A True Positive (TP) occurs when the model correctly identifies a positive class as positive, accurately recognizing the presence of a condition. Conversely, a False Positive (FP) happens when the model incorrectly predicts a negative class as positive, mistakenly indicating the presence of a condition when it is absent. Similarly, a True Negative (TN) is observed when the model accurately classifies a negative class as negative, correctly acknowledging the absence of a condition. Lastly, a False Negative (FN) takes place when the model incorrectly predicts a positive class as negative, thereby failing to detect a condition that is present.

3.3. Performance Metrics

The accuracy and effectiveness of predict results need to be measured with various metrics.

3.3.1. Recall

Recall expresses the rate at which the model correctly predicts positive classes (Equation 15). It is also known as sensitivity.

$$Recall={\displaystyle \frac{TP}{TP+FN}}$$

(15)

3.3.2. Precision

Precision measures how accurate the positive predictions of the model are. It is expressed as the ratio of true positive predictions to all positive predictions (Equation 16).

$$Precision={\displaystyle \frac{TP}{TP+FP}}$$

(16)

3.3.3. Specificity

Specificity refers to the rate at which the model correctly predicts negative classes as negative (Equation 17).

$$Specificity={\displaystyle \frac{TN}{TN+FP}}$$

(17)

3.3.4. F-Measure

It represents the harmonic mean of precision and recall (Equation 18). It is preferred especially in cases where the positive class is more important.

$$Fmeasure=2·{\displaystyle \frac{Precision·Recall}{Precision+Recall}}$$

(18)

3.3.5. Accuracy

The proportion of classes that the model predicts correctly (Equation 19) is defined as accuracy.

$$Accuracy={\displaystyle \frac{TP+TN}{TP+TN+FP+FN}}$$

(19)

3.3.6. Cohen's Kappa

Cohen's Kappa compares the classification performance of the model to the probability of making a random guess. It eliminates the effect of random guess. It is calculated by Equation 20. Here $P_o$ is defined as the observed accuracy rate ($(TP+TN)/Sum$), $P_e$ is defined as the expected accuracy rate (probability of random guess).

$$\kappa ={\displaystyle \frac{P_o-P_e}{1-P_e}}$$

(20)

3.4. Evaluation of Performance Metrics

The number of correct and incorrect predictions and performance metrics obtained as a result of the study are shown in Table 7 for RF, Table 8 for ANN, and Table 9 for the traditional method; all results are summarized in Table 10.

The accuracy values obtained for RF, ANN and the traditional method are given in Figure 12. The highest accuracy value was reached with the RF model with 89%. Accuracy rates of 87% were obtained for the traditional method and 86% for ANN.

3.5. ROC Curve and AUC

ROC curve (Receiver Operating Characteristic Curve) is a graph used to analyze how a classification model performs at different threshold values. It is drawn by calculating the True Positive Rate (TPR) and False Positive Rate (FPR) values. TPR (Sensitivity) is the rate at which the model correctly predicts true positives (Equation 21). FPR (False Positive Rate - Fall-out) is the rate at which the model falsely predicts true negatives as positives (Equation 22).

$$TPR={\displaystyle \frac{TP}{TP+FN}}$$

(21)

$$FPR={\displaystyle \frac{FP}{FP+TN}}$$

(22)

If the ROC curve is close to the upper left corner, it indicates that the model is performing well. If the curve is close to the 45-degree diagonal line, it means that the model is making random predictions. The ROC curve is given in Figure 13 for RF and in Figure 14 for ANN.

AUC (Area Under the Curve) refers to the area under the ROC curve and measures the overall classification performance of the model. If the AUC value is close to 1, the model makes excellent discrimination, if AUC = 0.5, the model makes random predictions, and if AUC < 0.5, the model makes poor classification. The AUC value is one of the most important metrics that show how well the model distinguishes positive and negative classes. A value between 0.7-0.8 is considered good, 0.8-0.9 is considered very good, and greater than 0.9 is considered an excellent model. The ROC curve and AUC are important metrics for comparing different models, evaluating the model in datasets with class imbalance, and determining the optimal threshold value. The ROC curves given in Figure 13 and Figure 14 were evaluated and the area under the graph was calculated as 0.932 for RF and 0.872 for ANN.

3.6. Feature Importance

The split number indicates how many times a feature is used for the split operation in the nodes in the trees. In each decision tree, the nodes divide the data into two subsets. During this splitting operation, one of the features is selected and splits the data according to a certain threshold value. The split number indicates how often the model uses that feature in the decision processes. A high split number indicates that the feature plays an important role in the classification or prediction process of the model. If the split number of a feature is very low (or even zero), this may indicate that the feature is generally unimportant to the model or is used less than other features.

The number of candidates indicates how many times a feature is considered a candidate for splitting. At each node, the model makes the best split decision from a randomly selected subset of features (candidate features). The number of candidates indicates how many times a feature is included in this subset. This number usually depends on the total number of nodes and the size of the randomly selected subset. A high number of candidates indicates that the feature is frequently included in randomly selected subsets and therefore frequently considered for splitting. However, being a candidate does not necessarily mean that a feature will be used. The model selects the feature that best splits.

Candidate > Split: Not every feature may be used for splitting in every case it is selected as a candidate. In other words, a feature may be frequently nominated in nodes, but may be used less because it does not perform better than other features for splitting.

Split ≈ Candidate: If a feature is both frequently nominated and frequently used for splitting, this means that the feature is very important for the model.

Split ≪ Candidate: If the candidate number of a feature is high but the number of splits is low, this may indicate that the feature generally falls short of the splitting criterion.

The split and candidate numbers in the first three levels of Random Forest trees are shown in Table 11. Here, the ratio of the split numbers in the first three levels to the total number of splits for each feature (input parameter) is significant in terms of feature importance. When the split number ratios given in Figure 15 are evaluated; It is seen that the parameters $q_c$, $a_{max}$, $f_s$, CSR and $I_c$ are more important than other parameters in the Random Forest model, respectively. Since feature importance is not directly calculated in ANN models compared to the RF module, only feature importance for the RF model was evaluated within the scope of the study.

4. Discussion

In this study, the potential of machine learning (ML) techniques for predicting seismic soil liquefaction was evaluated using Random Forest (RF) and Artificial Neural Network (ANN) models, and their performances were compared with the traditional Cone Penetration Test (CPT) based method. The findings clearly demonstrate the capability of ML algorithms in analyzing complex soil behaviors under seismic loading and highlight their advantages over conventional empirical approaches.

The analyses conducted to compare the ML models and the traditional method reveal significant differences in classification success:

• Accuracy Rates: According to the performance metrics, the RF model achieved the highest accuracy rate at 89%. This was followed by the traditional CPT-based method with 87% and the ANN model with 86%.
• ROC-AUC Analysis: When evaluating the Area Under the Curve (AUC) values, the RF model yielded an AUC of 0.932, whereas the ANN model scored 0.872. An AUC value greater than 0.9 indicates that the RF model possesses an "excellent" classification capability in distinguishing between liquefied and non-liquefied cases. Although the ANN model is considered a "very good" model with an AUC between 0.8 and 0.9, it fell behind the RF model.
• Overfitting Control: The overfitting problem, which is frequently encountered in decision tree methodologies, was successfully managed in the RF model through the use of "Bagging" and "Random Subspace" techniques. This structural advantage effectively reduces the variance of the model and increases its generalization capability, confirming why RF provided more robust predictions compared to the ANN model.

The fact that the traditional method achieved a relatively high accuracy of 87% proves its long-standing empirical validity. However, traditional deterministic approaches often struggle to fully capture the complex, nonlinear interactions between soil properties and seismic parameters.

Conventional methods typically rely on simplified assumptions and deterministic correlations to account for multifaceted parameters such as fines content or plasticity.

In contrast, machine learning algorithms can model nonlinear and multidimensional interactions flexibly, without the need for a predetermined mathematical form between input and output variables.

Furthermore, the dynamic framework of ML techniques allows them to continuously update and adapt as new data emerges, offering a distinct advantage over static traditional models in capturing the transient and variable nature of soil behavior.

Identifying which parameters play a more significant role in prediction models directly impacts engineering decision-making processes. The feature importance analysis applied to the RF model revealed the parameters that the model relied on most heavily:

The analysis demonstrated that the most influential variables used in the splitting operations of the trees were, respectively, cone tip resistance ($q_c$), peak ground acceleration ($a_{max}$), cone friction resistance ($f_s$), cyclic stress ratio (CSR), and soil behavior type index ($I_c$).

Considering that the fundamental mechanism underlying liquefaction is the cyclic shearing induced by seismic forces, it is theoretically consistent that $a_{max}$ and CSR are among the most critical parameters for the model.

Similarly, the high importance of direct field measurement metrics such as $q_c$ and $f_s$, which characterize the in-situ physical state of the soil, provides a strong physical validation that reinforces the model's accuracy.