Abstract: Using Monte Carlo simulation (MCS) based on stochastic finite element model to estimate failure probability of gravity dams suffers from a lack of computational efficiency. Surrogate model methods can effectively improve the computational efficiency, however, the conventional collocation methods (e.g., the three-sigma rule, Latin hypercube sampling (LHS) etc.) may suffer from a curse of dimensionality. Thusly, this paper proposes a collocation point selection within subsets of response probability space for surrogate modeling in gravity dam analysis and system reliability analysis method. The proposed method collocates samples in the response space by generalized subset simulation (GSS), then establishes surrogate model of each failure mode of gravity dams, and estimates failure probability by surrogate models based on MCS. This paper exhibits the collocation point selection result based on GSS by the application of a simple two-dimensional “four-boundary system” problem. Subsequently, the paper proves the effectiveness of the proposed method by the application of a gravity dam engineering. For the reliability analysis with a high dimensionality, nonlinear structural performance function, small failure probability, and multiple failure modes, e.g., gravity dam reliability analysis, the proposed method achieves a high accuracy in estimating failure probability. Additionally, it improves the computational efficiency and outperforms the generalized subset simulation method based on finite element model in terms of variance reduction. With the same computational costs, the proposed approach also outperforms the LHS based surrogate model method.