  This article presents numerical investigations on accuracy and convergence
properties of several numerical approaches for simulating steady state flows in
heterogeneous aquifers. Finite difference, finite element, discontinuous
Galerkin, spectral, and random walk methods are tested on one- and
two-dimensional benchmark flow problems. Realizations of log-normal hydraulic
conductivity fields are generated by Kraichnan algorithms in closed form as
finite sums of random periodic modes, which allow direct code verification by
comparisons with manufactured reference solutions. The quality of the methods
is assessed for increasing number of random modes and for increasing variance
of the log-hydraulic conductivity fields with Gaussian and exponential
correlation. Experimental orders of convergence are calculated from successive
refinements of the grid. The numerical methods are further validated by
comparisons between statistical inferences obtained from Monte Carlo ensembles
of numerical solutions and theoretical first-order perturbation results. It is
found that while for Gaussian correlation of the log-conductivity field all the
methods perform well, in the exponential case their accuracy deteriorates and,
for large variance and number of modes, the benchmark problems are practically
not tractable with reasonably large computing resources, for all the methods
considered in this study.
