Abstract: The statistics of uncertain rock and soil parameters can be updated with the information from different sources such as in-situ measurements and field observations via probabilistic back analysis, which can be further used for more realistic slope stability assessment. However, the inherent spatial variability of soil properties is almost not incorporated in the current probabilistic back analyses. This paper proposes an efficient approach by integrating multiple response-surface with subset simulation for probabilistic back analysis of slope failure in spatially variable soils. The Congress Street cut in Chicago with two important field observations including slope failure and approximate entry and exit regions of potential slip surfaces is taken as an example, and the posterior statistics of undrained shear strengths in three clay layers are estimated using the proposed approach. The results indicate the proposed approach can effectively back-analyze the posterior statistics of spatially varying soil properties at low-probability levels. The number of samples (N_l) in each intermediate step of subset simulation has an important effect on the posterior statistics of soil parameters, and the common choice of N_l = 500 cannot yield satisfactory results in general. In addition, the spatial variability of soil properties affects the posterior statistics of soil parameters significantly. The updated soil parameters follow non-stationary distributions in the slope profile when the spatial variability of soil properties is considered, which is in good accordance with geotechnical practice, while they still follow stationary distributions if the spatial variability of soil properties is ignored.