With advances in computational mechanics and computer science, the assessment of variability in engineering system responses has become feasible. This study evaluates the limit state of doubly eccentric loaded footing settlements in a 3D finite element model of a two-layer stratified soil domain. The Drucker–Prager constitutive model and Modified Cam Clay model are implemented within a high-fidelity finite element framework. The spatial variability of key soil parameters, including the reload path slope \(\kappa \) , the hydraulic permeability k (governing Darcy’s law), and the critical state line inclination c, is investigated using crude Monte Carlo simulations, with Latin hypercube sampling employed as the sampling strategy. Despite the strong nonlinearity of the system response, the Gaussian stochastic nature of the input variables is preserved in the output variables. Soil configurations with cohesive clay layers in the upper strata tend to produce higher statistical moments of the limit loads. In contrast, sandy soils, in several cases, yield higher mean values of maximum displacement and footing rotation. Similar trends are observed in the corresponding measures of output variability. Furthermore, the probabilities associated with the initiation point of the Meyerhof curve are estimated. The proposed computational framework can therefore support footing settlement design and decision-making processes.
