Based on the FLAC3D finite difference software, the dynamic reliability of seismic liquefaction response of sandy soil foundation covered with clay layer was studied, and a plane-strain model of clay-sandy soil foundation was established. By simulating the spatial variability of soil shear modulus, the area of liquefaction zone, excess pore water pressure and surface displacement were evaluated using the Monte Carlo simulation method and the generated non-Gaussian random field. The influence of different coefficient of variation (COV) on the reliability of seismic dynamic liquefaction of foundation was studied. The results show that the larger the COV of shear modulus, the more difficult it is to reduce the liquefaction range and the more dispersed the distribution range of peak liquefaction range; the horizontal displacement and maximum differential settlement of the surface increase with the increase of COV of shear modulus.