Geotechnical engineering applications are characterized by various sources of uncertainties, most of them attributed to the stochastic nature of soil parameters and their properties. In particular, soil's inherent random heterogeneity, inexact measurements and insufficient data necessitate numerical methods that incorporate the stochastic soil properties for a realistic representation of the soil behavior. In this paper, the process of consolidation of saturated soils is examined on the basis of the coupled u-p finite element formulation. A generalized Newmark implicit time integration scheme is implemented to treat the time integration of the coupled consolidation equations. A benchmark geotechnical engineering problem of a strip footing resting on a saturated soil layer is analyzed. The soil permeability coefficient k, as well as the elastic modulus E, are treated as lognormal random fields in two dimensions. The investigation of the effect of the spatial variability of the soil properties on the response of a footing-soil system is examined by means of the direct Monte Carlo simulation. The influence of the coefficient of variation and correlation length of the stochastic fields is quantified in terms of footing settlements, as well as excess soil water pore pressure. The effects of spatial variability of the permeability coefficient k and the elastic modulus E on the system response are demonstrated. It is shown that the footing differential settlement, along with generated excess pore pressures, is highly affected by the variation of the soil properties considered, as well as the correlation length of the underlying random fields.