We consider the effect of measuring randomly varying local hydraulic conductivities K(x) on one's ability to predict transient flow within bounded domains, driven by random sources, initial head, and boundary conditions. Our aim is to allow optimum unbiased prediction of local hydraulic heads h(x,t) and Darcy fluxes q(x,t) by means of their ensemble moments, 〈h(x,t)〉c and 〈q(x,t)〉c, conditioned on measurements of K(x). We show that these predictors satisfy a compact deterministic flow equation which contains a space-time integrodifferential "residual flux" term. This term renders 〈q(x,t)〉c nonlocal and non-Darcian so that the concept of effective hydraulic conductivity looses meaning in all but a few special cases. Instead, the residual flux contains kernels that constitute nonlocal parameters in space-time that are additionally conditional on hydraulic conductivity data and thus nonunique. The kernels include symmetric and nonsymmetric second-rank tensors as well as vectors. We also develop nonlocal equations for second conditional moments of head and flux which constitute measures of predictive uncertainty. The nonlocal expressions cannot be evaluated directly without either a closure approximation or high-resolution conditional Monte Carlo simulation. To render our theory workable, we develop recursive closure approximations for the moment equations through expansion in powers of a small parameter which represents the standard estimation error of natural ln K(x). These approximations are valid to arbitrary order for either mildly heterogeneous or well-conditioned strongly heterogeneous media. They allow, in principle, evaluating the conditional moments numerically on relatively coarse grids, without upscaling, by standard methods such as finite elements.