Abstract. Much of the stochastic analysis of flow field variability in heterogeneous aquifers in the literature assumes that the parameters in the associated stochastic flow equation are weakly (second order) stationary. On this basis, the spectral representation approach can then be used to quantify the variability of the flow fields given known covariance functions of the input parameters. However, the condition of second-order stationarity is rarely encountered in nature and is difficult to verify using the limited experimental data available. The purpose (or novelty) of this work, therefore, is to develop a new framework for modeling the variability of the flow fields that generalizes the stochastic theory that applies to stationary second-order random input parameters to intrinsic (nonstationary) random input parameters. In this work, the log hydraulic conductivity and log aquifer thickness are assumed to be intrinsic random functions for flow through heterogeneous confined aquifers of variable thickness. On this basis, semivariograms of depth-averaged hydraulic head and integrated specific discharge fields are developed to characterize the variability of flow fields. The application of the proposed stochastic theory to the case where the variability of a random input parameter can be characterized by a linear semivariogram model is provided.