On the accuracy of classic numerical schemes for modeling flow in saturated heterogeneous formations

Quantitative hydrogeology often relies on numerical modeling of flow and transport processes in the earth subsurface. Despite the richness of numerical schemes proposed in the literature most applications are performed by using a few very popular codes based on classical finite volume or finite element tech- niques. An important limitation of these numerical schemes is that they lead to solutions that do not satisfy the refraction law of streamlines at element (or volume) edges. This is not of great concern when the hydraulic conductivity K is spatially homogeneous, or varies smoothly within the computational domain. However, the solution may deteriorate in heterogeneous formations with high contrast between the hydraulic conductivity of adjacent computational cells. We analyze the performance of four widely used classic numerical schemes for solving the flow equation when they are applied to heterogeneous porous media. We first analyze the convergence of the numerical schemes to a known analytical solution in a simple heterogeneous field composed by 4 blocks with contrasting hydraulic conductivities. Then we compare the numerical solutions obtained in both Gaussian and exponential weakly heterogeneous logconductivity fields with existing analytical first- and second-order solutions in the variance of the logconductivity field, r2Y . Our analysis highlights that postprocessing the velocity field to enforce a pos- teriori the refraction law leads to biased results and that the performance of the numerical scheme depends on how mass conservation is discretized on the computational grid. Numerical schemes using inter-block conductivities, based for example on the harmonic mean, modify the spatial structure of the conductivity, with a negative impact on the structure of the velocity field.