The failure probability of geotechnical structures with spatially varying soil properties is generally computed using Monte Carlo simulation (MCS) methodology. This approach is well known to be very time-consuming when dealing with small failure probabilities. One alternative to MCS is the subset simulation approach. This approach was mainly used in the literature in cases where the uncertain parameters are modelled by random variables. In this paper, it is employed in the case where the uncertain parameters are modelled by random fields, because the spatial variability of the soil properties has proven to greatly affect the behavior of geotechnical structures and to induce a significant change in the variability of their responses. This is illustrated through the probabilistic analysis at the ultimate limit state (ULS) of a strip footing resting on a one- and two-layer purely cohesive soil with a spatially varying cohesion. The soil cohesion parameter was modeled as an anisotropic non-Gaussian (log-normal) random field using a square exponential autocorrelation function. The Expansion Optimal Linear Estimation (EOLE) method was used to discretize this random field. The deterministic model was based on numerical simulations using the finite difference software FLAC3D.