We consider the effect of measuring randomly varying hydraulic conductivities K(x) on one's ability to predict numerically, without resorting to either Monte Carlo simulation or upscaling, steady state flow in bounded domains driven by random source and boundary terms. Our aim is to allow optimum unbiased prediction of hydraulic heads h(x) and fluxes q(x) by means of their ensemble moments, ⟨h(x)⟩c and ⟨q(x)⟩c, respectively, conditioned on measurements of K(x). These predictors have been shown by Neuman and Orr <1993a> to satisfy exactly an integrodifferential conditional mean flow equation in which ⟨q(x)⟩c is nonlocal and non-Darcian. Here we develop complementary integrodifferential equations for second conditional moments of head and flux which serve as measures of predictive uncertainty; obtain recursive closure approximations for both the first and second conditional moment equations through expansion in powers of a small parameter &sgr;Y which represents the standard estimation error of ln K(x); and show how to solve these equations to first order in &sgr;Y2 by finite elements on a rectangular grid in two dimensions. In the special case where one treats K(x) as if it was locally homogeneous and mean flow as if it was locally uniform, one obtains a localized Darcian approximation ⟨q(x)⟩c≈-Kc(x)∇⟨h(x)⟩c in which Kc(x) is a space-dependent conditional hydraulic conductivity tensor. This leads to the traditional deterministic, Darcian steady state flow equation which, however, acquires a nontraditional meaning in that its parameters and state variables are data dependent and therefore inherently nonunique. It further explains why parameter estimates obtained by traditional inverse methods tend to vary as one modifies the database. Localized equations yield no information about predictive uncertainty. Our stochastic derivation of these otherwise standard deterministic flow equations makes clear that uncertainty measures associated with estimates of head and flux, obtained by traditional inverse methods, are generally smaller (often considerably so) than measures of corresponding predictive uncertainty, whichcan be assessed only by means of stochastic models such as ours. We present a detailed comparison between finite element solutions of nonlocal and localized moment equations and Monte Carlo simulations under superimposed mean-uniform and convergent flow regimes in two dimensions. Paper 1 presents the theory and computational approach, and paper 2 <Guadagnini and Neuman, this issue> describes unconditional and conditional computational results. ¿ 1999 American Geophysical Union