The natural variability of soil properties within geologically distinct and uniform layers has been proven to greatly affect soil behaviour and to induce significant variability in the predicted response. Previous studies concluded that small-scale heterogeneity greatly affects the liquefaction potential of saturated soil deposits, and provided geotechnical design guidelines to account for the effects of various characteristics of spatial variability. Those studies were based on two-dimensional analyses of soil liquefaction (in a vertical plane) assuming plane strain behaviour. Therefore the correlation distance of soil variability in a direction normal to the plane of analysis was implicitly taken as infinite (i.e. no variability in the third direction). In this study, a Monte Carlo simulation approach involving generation of sample functions of non-Gaussian, multivariate, multidimensional random fields and non-linear finite element analyses is used to investigate the effects of soil heterogeneity on the liquefaction potential of a 'stochastically heterogeneous' soil deposit subjected to seismic loading. To assess the 3D effects, Monte Carlo simulation results obtained for a 3D soil deposit are compared with corresponding results from 2D plane strain analyses. The calculations are performed for a range of seismic acceleration intensities, and the results are presented in terms of fragility curves expressing the probability of exceeding various thresholds in the response as a function of earthquake intensity.