Numerical benchmark study for flow in highly heterogeneous aquifers

Solving the flow problem is the first step in modeling contaminant transport in natural porous media formations. Since typical parameters for aquifers often lead to advection-dominated transport problems, accurate flow solutions are essential for reliable simulations of the effective dispersion of the solute plumes. The numerical feasibility of the flow problem for realistic parameters accounting for the heterogeneity of the aquifer and the spatial scale of the transport problem is addressed in a benchmark study.

The study aims to investigate the accuracy and the 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 two-dimensional benchmark flow problems. The heterogeneity of the aquifer system is described by log-normal hydraulic conductivity fields with Gaussian and exponential correlation structures. For given integral scale both correlation models predict the same effective coefficients, but they pose very different numerical challenges: while the Gaussian correlation ensures the sample-smoothness of the fields, the exponential correlation does not fulfil the theoretical requirements and the numerical representations of the samples are rather noisy.

Realizations of log-normal hydraulic conductivity fields are generated with a Kraichnan algorithm 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 variance of the log-hydraulic conductivity fields, which quantifies the heterogeneity, and for different numbers of random modes, which account for the spatial scale of the simulation. 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 exponential case their accuracy deteriorates and, for large variance and number of modes, the benchmark problems are practically not solvable with reasonably large computing resources, for all the methods considered in this study.