The direct effect of changes in applied total stresses on fluid-saturated geologic materials is reflected in fluctuations in pore pressure. The quantification of these pressure changes has traditionally been based on the classical poroelasticity theory, developed ignoring the natural variability in aquifer properties (parameters). Uncertainty of excess pore pressure prediction is therefore expected to be large when applying the classical model to field aquifer systems which generally have high degree of heterogeneity. Very limited attention to this issue has been paid. The motivation for this work is to quantify the uncertainty associated with the classical model, the spatial variability of the predicted excess pressure head in the field, where the groundwater flow is produced by the changes in total stress. The stochastic methodology used to develop the results of this study is based on the nonstationary spectral approach. It is found that the heterogeneity and correlation length scale of the log hydraulic conductivity process, and the pore compressibility parameter play essential roles in enhancing the variability in excess pressure head. In addition, the pore pressure head distribution predicted using the classical poroelasticity theory is subject to large uncertainty at a large depth in heterogeneous aquifers.