We consider the numerical prediction of transient flow in bounded, randomly heterogeneous porous media driven by random sources, initial heads, and boundary conditions without resorting to Monte Carlo simulation. After applying the Laplace transform to the governing stochastic flow equations, we derive exact nonlocal (integrodifferential) equations for the mean and variance-covariance of transformed head and flux, conditioned on measured values of log conductivity Y = ln K. Approximating these conditional moment equations recursively to second order in the standard deviation σ of Y, we solve them by finite elements for superimposed mean-uniform and convergent flows in a two-dimensional domain. An alternative conditional mean solution is obtained through localization of the exact moment expressions. The nonlocal and localized solutions are obtained using a highly efficient parallel algorithm and inverted numerically back into the time domain. A comparison with Monte Carlo simulations demonstrates that the moment solutions are remarkably accurate for strongly heterogeneous media with σ as large as 4. The nonlocal solution is only slightly more accurate than the much simpler localized solution, but the latter does not yield information about predictive uncertainty. The accuracy of each solution improves markedly with conditioning. A preliminary comparison of computational efficiency suggests that both the nonlocal and localized solutions for mean head and its variance require significantly less computer time than is required for Monte Carlo statistics to stabilize when the same direct matrix solver is used for all three (we do not presently know how using iterative solvers would have affected this conclusion). This is true whether the Laplace inversion and Monte Carlo simulations are conducted sequentially or in parallel on multiple processors and regardless of problem size. The underlying exact and recursive moment equations, as well as the proposed computational algorithm, are valid in both two and three dimensions; only the numerical implementation of our algorithm is two-dimensional.