The problem of flow through heterogeneous confined aquifers of variable thickness is analyzed from a stochastic point of view. The analysis is carried out on the basis of the integrated equations of the depth-averaged hydraulic head and integrated specific discharge, which are developed by integrating the continuity equation and equation for the specific discharge over the thickness, respectively. A spectrally based perturbation approach is used to arrive at the general results for the statistics of the flow fields in the Fourier domain, such as the variance of the depth-averaged head, and the mean and variance of integrated discharge. However, the closed-form expressions are obtained under the condition of steady unidirectional mean flow in the horizontal plane. In developing stochastic solutions, the input hydraulic conductivity parameter is viewed as a spatial random field characterized by the theoretical spatial covariance function. The evaluation of the closed-form solutions focuses on the influence of the controlling parameters, namely as a geometrical parameter defining the variation of the aquifer thickness and the correlation scale of log hydraulic conductivity, on the variability of the fluid fields. The application of the present stochastic theory to predict the total specific discharge under uncertainty using the field data is also provided.