Effective stress analyses based on the finite element method are often used as a reliable tool to predict liquefaction occurrence in soil-structure systems during earthquakes. In the analyses, the soil properties are typically specified by using a deterministic model although they intrinsically have spatial variability even in the case of horizontally layered ground. In this study, nonlinear finite element analyses under undrained conditions are performed to investigate the effects of soil heterogeneity on the liquefaction behavior of stochastically heterogeneous soil deposits subjected to seismic loading through a Monte Carlo simulation approach. A series of analyses has revealed that the heterogeneity of the shear wave velocity (or initial shear modulus) has no significant effect on the distribution of the computed excess pore water pressure (EPWP); while the mean value of the maximum EPWP ratio is partially influenced and becomes up to 20% less (in comparison with the deterministic case) by considering the spatial variability in the internal friction angle and the N value under the given seismic loading.