ABSTRACT The fundamental issue of bearing capacity of footings on anisotropic clays holds significant importance in geotechnical engineering. Previous investigations predominantly focused on deterministic analyses, disregarding the spatial variability of soil. A probabilistic analysis of the bearing capacity of footings is conducted in this paper, incorporating the spatial variability of anisotropic clays. To achieve this, Random Adaptive Finite Element Limit Analysis (RAFELA) and Monte Carlo simulations are utilised to capture the full spectrum of potential outcomes under parametric uncertainty. The impact of anisotropic soil strength variability is explored across three input parameters such as the anisotropic strength ratios, coefficients of variation, and dimensionless correlation lengths. In order to establish surrogate models capable of predicting random bearing capacity of anisotropic clays, Artificial Neural Network (ANN) models are developed. The use of the proposed ANN surrogate models presents a more convenient and computationally efficient approach for predicting the ultimate vertical load of footings on spatially random anisotropic clays.