<p>Accurate prediction of shallow foundation bearing capacity in layered soils is challenging due to complex interactions among geometric, mechanical, and stratigraphic factors. This study proposes a hybrid modeling framework that integrates Harris Hawks Optimization (HHO) with an artificial neural network (ANN) to predict the normalized bearing capacity (<i>q</i>_u/<i>γ</i>₁<i>B</i>) of strip footings embedded in bi-layered soil systems. A comprehensive dataset generated from high-fidelity isogeometric analysis is used to train the model based on nine dimensionless, physically meaningful input variables describing soil properties, loading conditions, embedment effects, and footing interface characteristics. HHO is employed to systematically optimize ANN hyperparameters, resulting in robust convergence and strong generalization performance. The proposed model is benchmarked against several optimization-based ANN configurations under a unified evaluation protocol and further validated against published numerical limit analysis studies, confirming its physical consistency and reliability. To enhance interpretability, explainable artificial intelligence techniques, including SHAP/DeepSHAP and partial dependence analyses, are applied to identify dominant controlling variables and to elucidate distinct response mechanisms between sand–sand and sand–clay systems, as well as nonlinear interaction effects. Building on the validated surrogate, a web-based prediction tool is developed to support rapid parametric exploration and improve practical accessibility for engineering applications.</p>