This paper explores the applicability of variability response functions to nonlinear soil–structure interaction problems, focusing on the impacts of spatially variable soil properties on foundation reliability regarding the settlement response. An estimation scheme is proposed to obtain the response functions, which involves a periodic function to approximate the relationship between foundation response and phase angles representing the soil variability patterns. Using a single set of realizations, the response function approach enables the evaluation of foundation reliability under various spatial variability patterns, including different autocovariance distances in three dimensions and the rotated anisotropy features of soil variations. This leads to significant reductions in computational demands compared to previous attempts of random field modelling, which often involved individual Monte Carlo simulations for each combination of spatial variability parameters. The proposed approach is applied to both shallow foundation and piled foundation cases, illustrating its range of applicability. For linear-elastic systems, the approach is shown to be effective for various coefficients of variation in soil variability. For elastic–plastic pile group analyses, the approach leads to efficient evaluation of the statistics of average and differential settlements of the pile group, both of which compare favourably with conventional random field simulation techniques. Since it does not require multiple random field realizations, the approach is particularly efficient in identifying the worst-case scenario of autocovariance distances that corresponds to the largest uncertainty in foundation response. This can become a useful tool for conservative reliability assessments under a project setting, since site-specific estimates of random field parameters are often imprecise due to limited site data.