In this paper we investigate the combined effect of nonstationarity in log hydraulic conductivity and the presence of boundaries on flow in heterogeneous aquifers. We derive general equations governing the statistical moments of hydraulic head for steady state flow by perturbation expansions. Due to their mathematical complexity, we solve the moment equations by the numerical technique of finite differences. The numerical approach has flexibility in handling (moderately) irregular geometry, different boundary conditions, various trends in the mean log hydraulic conductivity, spatial variabilities in the magnitude and direction of mean flow, and different covariance functions, all of which are important factors to consider for real-world applications. The effect of boundaries on the first two statistical moments involving head is strong and persistent. For example, in the case of stationary log hydraulic conductivity the head variance is always finite in a bounded domain while the head variance may be infinite in an unbounded domain. As in many other stochastic models, the statistical moment equations are derived under the assumption that the variance of log hydraulic conductivity is small. Accounting for a spatially varying, large-scale trend in the log hydraulic conductivity field reduces the variance of log hydraulic conductivity. Although this makes the conductivity field nonstationary and significantly increases the mathematical complexity in the problem, it justifies the small-variance assumption for many aquifers. ¿ 1998 American Geophysical Union