The field and laboratory evidence of nonlinear soil behaviour, even at small strains, emphasizes the importance of employing nonlinear methods in seismic ground response analysis. Additionally, determination of dynamic characteristics of soil layers always includes some degree of uncertainty. Most of previous stochastic studies of ground response analysis have focused only on variability of soil parameters, and the effect of soil sample location has been mostly ignored. This study attempts to couple nonlinear time-domain ground response analysis with variability of soil parameters considering existing boreholes’ location through a geostatistical method using a program written in MATLAB. To evaluate the efficiency of the proposed method, stochastic seismic ground responses at construction location were compared with those of the non-stationary random field method through real site data. The results demonstrate that applying the boreholes’ location significantly affects not only the ground responses but also their coefficient of variation (COV). Furthermore, the mean value of the seismic responses is affected more considerably by the values of soil parameters at the vicinity of the construction location. It is also inferred that considering boreholes’ location may reduce the COV of the seismic responses. Among the surface responses in the studied site, the values of peak ground displacement (PGD) and peak ground acceleration (PGA) reflect the highest and lowest dispersion due to variability of soil properties through both non-stationary random field and geostatistical methods.