## 4 Experiments
### Results  
The correlation matrix of features depicted in Figure 11 shows delicate relationships among the geophysical properties that control liquefaction risk. The hvsr_mean and hvsr_std are highly correlated, indicating that variability and the shape of resonance go together; sites with high resonance tend to show increased variability of seismic wave amplification, which may mean perilous vibration. The alpha parameter is less related to other parameters, providing information on harmonic behavior that does not directly relate to site density but is relevant to dynamic response estimation. The gamma parameter is highly correlated with \(K_{g}\), indicating that sites with greater seismic amplification are more vulnerable to liquefaction. The negative correlation with \(f0\) suggests that softer soils (lower fundamental frequency) tend to experience higher amplification, further increasing the liquefaction risk. The \(K_{g}\) (vulnerability index) is strongly negatively correlated with f0, which means that lower fundamental frequency (softer, more liquefiable ground) is  
Figure 10: Resilience and anomalies in seismic vulnerability index (\(K_{g}\)) distribution. Figure 9: Distribution of amplitude (Af) across liquefaction and non-liquefaction sites. associated with a higher \(K_{g}\). In other words, resonance increases dynamic response and increases liquefaction potential. The Z parameter is highly correlated with hvsr_mean and hvsr_std, which means liquefaction susceptibility is related to seismic energy amplification with both horizontal and vertical spectral components, so sites that have a more stable energy distribution or higher risk of amplification. Finally, the high correlation between Af and \(K_{g}\) means high amplification factors, especially with lower f0 values, are at fault for liquefaction susceptibility. Finally, liquefaction susceptibility depends on the frequency, amplification, and energy distribution. The relationships among \(K_{g}\), f0, Af, and Z necessitate a multi-dimensional methodology that accounts for the interplay of amplification and resonance with site dynamics to provide a comprehensive liquefaction hazard assessment. We utilized 60 samples for training, with dimensions of (60, 24), where 24 denotes the number of features. The XGBoost model was configured with the following hyperparameters: n_estimators = 30, max_depth = 5, random_state = 42, use_label_encoder = False, and eval_metric = 'logloss'. These parameters were set to balance model complexity and prevent overfitting, given the size of the dataset. Figure 12 presents the confusion matrices of three classification models--MLP, XGBoost, and the proposed stacked model--for liquefaction and non-liquefaction event prediction. The confusion matrix provides an overview of the model's accuracy. The highest number of true positives and very few false negatives in the results of the best-performing stacked model showcase its ability to correctly identify and minimally misclassify events related to liquefaction. The true-negative and false-positive distributions are nearly balanced, indicating good results in non-liquefaction cases. Figure 11: Correlation analysis of geophysical features for liquefaction susceptibility assessment in the Rif region. XGBoost, which ran second in performance, presents a slightly lower true-positive rate than the stacked model, implying more false negatives. Still, this model performed better than MLP, presenting fewer false positives and better handling of non-liquefaction events. The worst was the MLP model, which had more false negatives and false positives. That is to say that the MLP model is not very efficient in detecting liquefaction events and misclassifies non-liquefaction events, showing its weakness in dealing with an imbalanced dataset. Table 1 highlights the classification accuracies of the MLP, XGBoost, and suggested stacking models.