The objective of this study is to quantify the effect of the spatial variability of soil properties on site response and ground motion coherency. A 2D random field theory is used to simulate soil parameters, including shear wave velocity (Vs) and nonlinear stress–strain curves. Different levels of spatial variability are controlled by coefficient of variation (COV) and correlation length (CL). The wave prorogation is simulated through the 2D finite difference software FLAC2D. Different types of analysis conditions, such as 1D or 2D model, linear or nonlinear analysis, and input motion, are investigated. The increase of Vs variability leads to a decrease in the mean amplification ratio (surface response spectrum to bedrock response spectrum) but an increase in its variation. The effect of the variability of soil nonlinearity on the mean amplification ratio is marginal compared with that induced by Vs variability. However, as the ground motion intensities increase (or more nonlinearity induced), the standard deviation of the simulation normalized by the baseline response increases. Compared with the 2D result, the 1D site response analysis that is commonly used in practice may potentially underestimate the mean amplification ratio but overestimate its variability due to the ignorance of the horizontal soil property variability. In the evaluation of the spatially correlated site response, the ground motion coherency decays significantly, as the COV of Vs and the ground motion intensity increase. The effect of CL on ground motion coherency is minor and limited to a separation distance within two‐times of CL. Therefore, a coherency model, which eliminates the CL term but includes the ground motion intensity (i.e., nonlinear effect) parameter, is proposed in this study. It is revealed to yield a better agreement with the simulation results than other coherency models.