The problem of depth-averaged groundwater flow in heterogeneous aquifers is looked at from a stochastic point of view. The Galerkin finite element version of the linearized stochastic equation governing the flow is solved analytically in time using an eigenvalue-eigenvector technique. The stochastic solution, which is valid for small variance of aquifer log transmissivity, relates spatial and temporal variabilities of aquifer head to aquifer heterogeneity, stochastic recharge, and random initial head. The computational effort requires the evaluation of integrals of matrices whose elements are linear functions of the nodal mean heads. The solutions are obtained for the exact (i.e., continuous) and quasi-steady approximation of the mean head. Aquifer-head temporal covariances are evaluated using standard matrix operations once the storage matrix is inverted and the associated generalized eigenvalue problem is solved. Implication of the level of spatial discretization on the performance of the solution is examined, and the influence of aquifer heterogeneity on spatial and temporal variances of hydraulic heads is investigated in the presence of a semipervious boundary. ¿ American Geophysical Union 1995