We consider the effect of measuring randomly varying hydraulic conductivities K(x) on numerical prediction of transient flow in bounded domains driven by random source and boundary terms. Tartakovsky and Neuman(1,2) derived exact integro-differential equations for optimum unbiased predictors of hydraulic quantities and associated predictive uncertainty. They also explored ways to localize the equations so as to render them differential. The authors developed recursive closure approximations for the nonlocal equations through expansion in powers of a, the standard estimation error of In K(x). Here we solve these recursive equations numerically to first order in sigma(Y)(2) using a parallel finite element scheme in the complex Laplace domain and numerical inversion. Results for superimposed meanuniform and convergent flows in two dimensions compare well with standard Monte Carlo simulations for sigma(Y)(2) considerably higher than 1. Localization is shown to work well except near the point sink; it however does not provide information about predictive uncertainty. The degree to which parallelization enhances computational efficiency is explored.