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 [19], accurate flow solutions are essential for reliable simulations of the effective dispersion of the solute plumes. The numerical feasibility of the flow problem in realistic conditions has been addressed in a pioneering work by [1], who pointed out that the discretization steps must be much smaller than the heterogeneity scale, which in turn is much smaller than the dimension of the computational domain. Such order relations also served as guidelines for implementing solutions with finite difference methods (FDM) [22, 42, 79] or finite element methods (FEM) [7, 57, 68] used in Monte Carlo simulations. In spectral approaches [32, 46] or probabilistic collocation and chaos expansion approaches [45, 47] the spatial resolution is also related to the number of terms in the series representation of the coefficients and of the solution of the flow problem. Despite the advances in numerical methods, computing accurate flow solutions in case of highly heterogeneous formations faces computational challenges in terms of code efficiency and computational resources [see e.g., 33, 45].

Field experiments show a broad range of heterogeneities of the aquifer systems characterized by variances ${\sigma }^2$ of the log-hydraulic conductivity from less than one, for instance, at Borden site, Ontario [59, 61], up to variances between ${\sigma }^2=3.4$ and ${\sigma }^2=8.7$ at the Macrodispersion Experiment (MADE) site in Columbus, Mississippi [9, 60]. The geostatistical interpretation of the measured conductivity data is often based on a log-normal probability distribution with an assumed exponential correlation model [9, 30, 60]. But, since dispersion coefficients do not depend on the shape of the correlation function [19] and are essentially determined by the correlation scale [59], a Gaussian model can be chosen as well [42, 79]. Moreover, experimental data from Borden, MADE, and other sites can be reinterpreted by using a multifractal analysis and power-law correlations [49, 53, 54]. Random fields with power-law correlations can be generated as linear combinations of random modes with either exponential or Gaussian correlation [25, 77]. Therefore, the influence of the correlation structure on the accuracy of the numerical solutions of the flow problem can be generally analyzed by considering Gaussian and exponential correlations of the log-conductivity field.

Numerical simulations of flow and transport in groundwater were carried out with a FEM method for increasing log-conductivity variance from small to moderate values, e.g., up to ${\sigma }^2=1.6$ by [7]. Similar investigations were performed with FDM methods for variances as high as ${\sigma }^2=9$ in [22] or ${\sigma }^2=16$ in [42]. Though an exponential shape of the log-conductivity correlation is considered in most cases, a Gaussian shape which yields smoother fields could be preferred for numerical reasons [79]. Different approaches were validated through comparisons with linear theoretical approximations of first-order in ${\sigma }^2$ and the numerical methods were further used to explore the limits of the first-order theory. A direct validation of the flow solver by comparisons with analytical manufactured solutions has been done by [42], for both Gaussian and exponential correlations with the same variance ${\sigma }^2=1$ of the log-conductivity field.

Accurate FEM methods were developed to ensure the consistency of the numerical fluxes with Darcy’s law [56, 17], the mass conservation of the numerical scheme [58], and the head continuity at the interface between two materials with contrasting hydraulic properties [12]. With focuss on satisfying the refraction law at the interfaces between blocks with different hydraulic conductivity, [12] used advanced FEM and velocity postprocessing approaches to solve the flow problem for weak and moderately heterogeneous aquifer structures. The convergence of the numerical schemes is assessed by comparisons with analytical solutions for block-constant hydraulic conductivity with sharp contrast. For further validation, the convergence to the same velocity variances of the predictions obtained with different numerical schemes was investigated through a Monte Carlo approach. While almost no discrepancy was observed for the Gaussian correlation model, for the exponential model one obtains large relative differences for the variance of both longitudinal and transverse velocity in case the largest variance of the log-hydraulic conductivity considered in their study, ${\sigma }^2=2$.

Serious challenges are posed by exponential correlation models with higher variances, when conductivity differences may span several orders of magnitude. Using a spectral collocation approach based on expansions of the solution in series of compactly supported and infinitely differentiable “atomic” functions, [33] tested the flow solver for variances up to ${\sigma }^2=8$. The accuracy of the solution was assessed with mass balance errors between control planes, the criterion of small grid Péclet number (defined with advection velocities given by the derivatives of the log-hydraulic conductivity and constant unit diffusion coefficient highlighted by expanding the second order derivatives in the flow equation), as well as by corrections of the solution induced by the increase of the number of collocation points. It was found that increasing the resolution of the log-hydraulic conductivity requires increasing the resolution of the solution to fulfill the accuracy criteria. While in case of Gaussian correlation of log-hydraulic conductivity the desired accuracy is achieved, for ${\sigma }^2$ up to 8, with relatively small number of collocation points per integral scale, in case of exponential correlation, which requires increasing resolution with increasing ${\sigma }^2$, the required resolution of the solution may lead to problem dimensions which exceed the available computational resources [33, Section 5.2].

In line with the investigations discussed above, we present in this article a comparative study of several numerical schemes for the Darcy flow equations, based on different concepts. The schemes will be tested in terms of accuracy and convergence behavior on benchmark problems designed for a broad range of parameters describing the variance and the spectral complexity of the log-hydraulic conductivity fields with both Gaussian and exponential correlation structures used as models for heterogeneous aquifers.

Comparisons of numerical schemes for flow problems have been made for constant hydraulic parameters, the goal being to solve nonlinear reactive transport problems [14, 15], for constant parameters on sub-domains and complex geometry, in order to simulate flow in fractured porous media [23, 28], by considering functional dependence of the hydraulic conductivity on temperature, in case of coupled flow and heat transport [34], or in case of heterogeneous hydraulic conductivity fields with low variability [57]. Such comparative studies use benchmark problems and code intercomparison for the validation of the numerical methods [e.g., 15] or, whenever available, comparisons with analytical solutions for the purpose of code verification [e.g., 34].

In the present study, we adopt a different strategy. The hydraulic conductivity is generated as a realization of a Kraichnan random field [40] in a closed form, consisting of finite sums of cosine random modes, with a robust randomization method [41, 42]. Unlike in turning bands [51] or HYDRO_GEN [8] methods, the Kraichnan approach allows a direct control of the accuracy of the field, by increasing the number of modes, and of its variability, through the variance of the log-hydraulic conductivity random field. Moreover, the explicit analytical form of the generated field allows the computation of the source term occurring when a manufactured analytical solution is forced to verify the flow equation [11, 29, 42, 55, 57, 63, 66]. A realization of the hydraulic conductivity and the source term are computed by summing up cosine modes with parameters given by a fixed realization of the wave numbers and phases produced by the randomization method. These parameters are computed with C codes and stored as data files. Codes written in different programming languages, for instance Matlab, C++, or Python load the same data files to construct precisely the same manufactured solution to be used in code verifications. One obtains in this way a benchmark to evaluate the performance of different numerical methods for a wide range of variances of the log-hydraulic conductivity and of the number of random modes.

Considering, beside the variance of the log-conductivity, the influence of the number of modes is crucial for several reasons. This quantity controls the accuracy of the generated hydraulic conductivity fields, and thus ensures accurate statistical simulations of the solutions of the partial differential equations which govern the Darcy flow. It also affects the ergodicity properties [41, 43] of the simulation methods which infer statistics by performing spatial averages [78]. The number of modes in randomization methods increases with the order of magnitude of the ratio between the maximum and the minimum length scales of the random field structure which has to be captured in numerical simulations. While the maximum length scale should be larger than the correlation length of the random field and is practically given by the scale of the simulation, the minimum scale is a “smoothness microscale” determined by the correlation length of the spatial increments of the field, which is an analogous correlation length for the gradient of the random field [e.g., 41, Sect. 2]. The relation between the number of modes and the scale ratio is determined by the correlation structure of the random field. For instance, in case of power-law correlations it has been shown analytically that this relation is logarithmic or at most linear, depending on the details of the randomization method [41]. In case of an exponential correlation and a small variance of the log-hydraulic field, it has been demonstrated by numerical investigations that, in order to ensure the self-averaging behavior of the simulated transport process, the number of modes in the Kraichnan routine has to be of the order of the number of correlation lengths traveled by the solute plume [27, Fig. 4]. Hence, while the variance characterizes the spatial heterogeneity of the aquifer system, the number of modes used in the Kraichnan randomization method to obtain accurate approximations of the hydraulic conductivity increases with the spatial scale of the flow and transport problem.

We compare the accuracy of the FDM, the FEM, the discontinuous Galerkin method (DGM), the Chebyshev collocation spectral method (CSM), and a newly developed “global random walk” method (GRW) by solving two dimensional problems. A Galerkin spectral method (GSM) is tested only in the one-dimensional case. Code verification is done by evaluating the errors with respect to manufactured analytical hydraulic head solutions for combinations of parameters of the Kraichnan routine, i.e. variances ${\sigma }^2$ of the log-hydraulic conductivity between 0.1 and 6 and two different number of modes, $N={10}^2$ and $N={10}^3$. Estimations of the convergence order are obtained numerically by successive refinements of the discretization. The codes verified in this way are further used to solve the homogeneous equation for the hydraulic head and to compute the Darcy velocity. Ensembles of 100 realizations of the flow solution are computed and used in Monte Carlo inferences of the statistics of the hydraulic head and of the components of the Darcy velocity. The numerical methods are validated by comparisons with first-order theoretical results [6, 19, 52].

In order to test the numerical methods in conditions of hydraulic conductivity with increasing variability, comparisons are done for log-hydraulic conductivities with exponential and Gaussian correlations. Testing the numerical schemes for these correlation models is highly relevant from the point of view the numerical difficulty of solving the flow problem. Indeed, while a Gaussian correlation ensures the smoothness of the generated hydraulic conductivity samples, in case of exponential correlation the smoothness significantly deteriorates with the increase of the number of modes, the generated samples approaching non-differentiable functions in the limit of an infinite number of modes [72, 74].

The remainder of this paper is organized as follows. Section 2 presents the benchmark problems, the Kraichnan randomization method, and some considerations of the sample-smoothness of the generated hydraulic conductivity fields with Gaussian and exponential correlation. The numerical schemes considered in this study are described in Section 3. The results of the code verification by manufactured solutions are presented in Section 4.1 and the convergence results in Section 4.2. Section 5 presents Monte Carlo inferences and comparisons with first-order theory. Finally, some conclusions of the present study are drawn in Section 6. A presents the computation of the manufactured solutions and B contains tables with estimated order of convergence. Supplementary materials S1–S4 associated with this article are available in the online version.

We consider the two-dimensional steady-state velocity field in saturated porous media with constant porosity determined by the gradient of the hydraulic head $h(x,y)$ according to Darcy’s law [19]

where $K(x,y)$ is an isotropic hydraulic conductivity consisting of a fixed realization of a spatial random function. The velocity solves a conservation equation $\partial V_x/{\partial }_x+\partial V_y/{\partial }_y=0$, which after the substitution of the velocity components (1) gives the following equation for the hydraulic head,

We are looking for numerical solutions of the equations (1)-(2) in a rectangular domain, $(x,y)\in \mathrm{\Omega }=[0,L_x]\times [0,L_y]$, subject to boundary conditions given by

where $H$ is a constant head.

The hydraulic conductivity fields $K$ considered in this study are isotropic, log-normally distributed, with exponential and Gaussian correlations. The random log-hydraulic conductivity field $Y=\ln (K)$ is specified by the ensemble mean $⟨Y⟩$, the variance ${\sigma }^2$, and the correlation length $\lambda$. The Kraichnan algorithm [40] used to generate samples of fluctuating part of the log-hydraulic conductivity field $Y^{\prime }=Y-⟨Y⟩$ is a randomized spectral representation of statistically homogeneous Gaussian fields [41]. A general randomization formula, valid for both exponential and Gaussian correlation types, is given by a finite sum of $N$ cosine functions,

The probability density of the random vector $(k_{i,1},k_{i,2})$ is given by the Fourier transform of the correlation function of the statistically homogeneous random field $Y^{\prime }$ divided by the variance ${\sigma }^2$ and the phases ${\varphi }_i$ are random numbers uniformly distributed in $[0,2\pi ]$ [42, Eq. (22)].

The smoothness of the $K$ field is essentially determined by the correlation of the random field $Y^{\prime }$. Sufficient conditions for sample smoothness are given by the existence of higher order derivatives of the correlation function $C(𝐫)=⟨Y^{\prime }(𝐱)Y^{\prime }(𝐱+𝐫)⟩$. For instance, in the one-dimensional case, if $C(r)$ has a second order derivative at $r=0$, then the one-dimensional field $Y^{\prime }(x)$ is equivalent to a field with samples that are continuous with probability one in every finite interval. Moreover, if the fourth derivative $C^{(4)}(0)$ exists, $Y^{\prime }(x)$ is equivalent to a field the sample functions of which have, with probability one, continuous derivatives in every finite interval. These sufficient conditions can hardly be relaxed and one expects that they are very close to the necessary conditions [18, 83].

The smoothness conditions of [18] are fulfilled for instance if the correlation of $Y^{\prime }$ has a Gaussian shape $C(r)\sim e^{-r^2}$ and is infinitely differentiable, which results in a random field (5) with continuous sample derivatives. A counterexample is the field $Y^{\prime }$ with an exponential correlation $C(r)\sim e^{-r}$, which is not differentiable at the origin. Since, rigorously speaking, the differentiability of $C(r)$ at the origin is only a sufficient condition, we proceed with numerical investigations on the sample-smoothness of the $K$ field in case of Gaussian and exponential correlations of the log-hydraulic conductivity $Y$ with $\lambda =1$ and ${\sigma }^2=0.1$. The estimations of the spatial derivative $dK/dx$ of a given realization of the numerically generated one-dimensional $K$ field with $N={10}^2$ random modes in (5), computed at $x=1+i\mathrm{\Delta }x,i=1\mathrm{\dots }100$, are presented in Figs. 2 and 2. In case of Gaussian correlation the estimations approach the exact value of $dK/dx$ at $x=1$ obtained by differentiating the expression (22) with $y=1$ fixed. In case of exponential correlation of $Y$, the realization of the $K$ field apparently approaches a non-differentiable function for $\mathrm{\Delta }x\le {10}^2$, then, by further decreasing $\mathrm{\Delta }x$ one observes a smooth behavior and the approach to the true value of $dK/dx$ at $x=1$.

The impact of increasing the number of modes is illustrated by the convergence behavior of the numerical estimates of $dK/dx$ presented in Figs. 4 and 4. One remarks that while for Gaussian correlation the number of modes has no impact, the convergence towards the exact value of the derivative being practically achieved for $\mathrm{\Delta }x={10}^{-2}$, in case of exponential correlation the resolution has to be increased four orders of magnitude ($\mathrm{\Delta }x={10}^{-6}$) to obtain an accurate approximation of the derivative when $N$ increases one order of magnitude.

In order to investigate the robustness and the limits of applicability of the numerical schemes tested in this study, we consider highly oscillating coefficients $K$ generated by (5) for isotropic correlations of the random field $Y^{\prime }$, $C(r)=⟨Y^{\prime }(x,y)Y^{\prime }(x+r_x,y+r_y)⟩$, where $r={(r_x^2+r_y^2)}^{1/2}$. As correlation functions, we consider the two extreme cases presented above, i.e., the Gaussian correlation,

which yields analytical samples of $Y^{\prime }$, and, respectively, the exponential correlation,

which, in the limit of infinite number of modes $N$, produces non-smooth samples.

To compute the hydraulic conductivity $K=\exp (Y)=\exp (⟨Y⟩)\exp (Y^{\prime })$, the geometric mean $K_g=\exp (⟨Y⟩)$ has to be specified either from experimental data or through its relation with other parameters of the flow problem, such as the mean slope of the hydraulic head, $-H/L_x=-J$, and the mean velocity $U$. In this study we consider the second approach. According to the boundary conditions (3-4), the mean velocity $(U,0)$ is aligned to the longitudinal axis $x$. As an input in the field generator, we use a mean hydraulic conductivity $⟨K⟩$ related according to Darcy’s law (1) to the mean velocity and the mean slope by $U=⟨K⟩J$. Further, as far as we do not compare velocity fields for different ${\sigma }^2$, we generate random conductivity fields by using the relation between the arithmetic and the geometric mean of the log-normal distribution, $⟨K⟩=K_g\exp ({\sigma }^2/2)$. The arithmetic mean approximates the effective hydraulic conductivity $⟨Kh⟩/J$ if ${\sigma }^2\ll 1$ [30, 82]. Higher order corrections were also derived as asymptotic approximations for small ${\sigma }^2$ [24]. But our numerical tests have shown that neither of these corrections is appropriate for statistical inferences, because they lead to underestimations of the mean velocity. Therefore, to compare results obtained with different values of the parameter ${\sigma }^2$ presented in Section 5, we normalize the velocity components by $K_{g}J$ [see also 42].

The mean hydraulic conductivity is fixed to $⟨K⟩=15$ m/day, a value representative for gravel or coarse sand aquifers [19]. The correlation length $\lambda$ in (6) and (7) is set to $\lambda =1$ m and the dimensions of the spatial domain are given in $\lambda$ units. In order to compare the results obtained by different schemes for a given realization of the field $K$, the sets of wavenumbers $k_{i,1}$, $k_{i,2}$, and phases ${\varphi }_i$ are computed with the algorithms for Gaussian and exponential correlation functions described in [42, Sects. 5.2 and 5.3] implemented as C-codes. These coefficients, computed for $N_{max}={10}^4$ modes $i$, are stored as data files. Samples of the $K$ field are further constructed for the desired numbers of modes $N\le N_{max}$ and values of the variance ${\sigma }^2$, according to (5). The same numerical parameters can be used in either Matlab, C++, or Python codes, to ensure identical values of the $K$ fields, in the limit of double precision, irrespective of the programming language. For the purpose of statistical inferences, ensemble of solutions of the flow equations are computed with the same algorithms implemented as functions in each particular scheme.

The longitudinal dimension of the domain, $L_x$, has to be tens of times larger than $\lambda$ (for instance, [1] recommend $L_x/\lambda =25$). Constrained by considerations of computational tractability, we fix the dimensions of the domain to $L_x=20$ and $L_y=10$. The constant head at the inflow boundary is fixed to $H=1$. To enable comparisons of the two dimensional solutions obtained with the FDM, FEM, DGM, and GRW numerical schemes tested in the present benchmark, the number of discretization points is set to $n_x={10}^3+1$, in the longitudinal direction, and to $n_y=5\cdot {10}^2+1$ in the transverse direction. Consequently, the space steps are fixed at $\mathrm{\Delta }x=\mathrm{\Delta }y=2\cdot {10}^{-2}$, which, for unit correlation length considered in this study, gives $\lambda /\mathrm{\Delta }x=50$. This value is larger than the resolution of the log-conductivity field usually considered in similar finite difference/element approaches reported in the literature [see e.g., 8, 12, 22, 42]. The fixed resolution considered in this study is appropriate for flow solutions in case of Gaussian correlation of the $\ln (K)$ field but could be to coarse in case of exponential correlation (see Figs.2 – 4).

In case of spectral approaches, the number of collocation points, for both CSM and GSM schemes, is chosen after preliminary tests as large as necessary, possible optimal, to ensure the desired accuracy (see Section 3.4).

The codes based on the numerical schemes considered in this study are first verified by comparisons with manufactured solutions. The approach consists of choosing analytical functions of the dependent variables, $h(x,y)$, on which one applies the differential operator from the left hand side of Eq. (2), the result being a source term $f$. Hence, $h(x,y)$ is an exact solutions of Eq. (2) with a source term $f$ added in the right hand side and modified boundary functions according to the chosen manufactured solution [e.g., 63].

An example of choosing the appropriate manufactured solution for the flow solver within a numerical setting similar to that of the present study can be found in [42, Appendix B]. In order to test the numerical approximations of the differential operators in the flow equation (2), the manufactured solutions used in this study are simple, but not trivial analytical functions with nonzero derivatives up to the highest order of differentiation. The construction of the manufactured solution and of the source term $f$ is presented in A. For the purpose of code verification we chose 10 parameter pairs formed with ${\sigma }^2=0.1,\mathrm{\hspace{0.33em}1},\mathrm{\hspace{0.33em}2},\mathrm{\hspace{0.33em}4},\mathrm{\hspace{0.33em}6}$, and the two numbers of modes considered in Figs. 4 and 4, $N={10}^2$ and $N={10}^3$. A larger set of 21 parameter pairs, also including ${\sigma }^2=8$, ${\sigma }^2=10$, and $N={10}^4$, has been considered in the extended version of this study [3].

Note that the two values of $N$ considered here are relevant for practical applications of the numerical codes. A number of modes of the order ${10}^2$ is required to obtain accurate simulations of the transport process over tens of correlation scales even for a small variance ${\sigma }^2=0.1$ and a linearized solution of the flow equation [27]. The larger number of modes $N={10}^3$ is needed, in the same conditions, for simulation lengths of hundreds of correlation scales. In the latter case, the small spatial domain $\mathrm{\Omega }$ considered above can be thought of as a subdomain in a domain decomposition approach. On the other hand, the problem is linear and the correlation of the log-hydraulic conductivity only depends on spatial lags scaled by the correlation length. Thus, the results presented in this benchmark study also can represent solutions to problems with spatial dimensions scaled by the correlation length.

The code verification is completed by convergence tests aiming at comparing the computational order of convergence with the theoretical one. In case of FDM, FEM, DGM, and GRW approaches, computational orders of convergence are estimated by computing error norms of the solutions obtained by successively halving the spatial discretization step with respect to a reference solution obtained with the finest discretization [see e.g. 65, 50]. The approach is different in case of spectral methods, where the convergence is indicated by the decrease towards the round-off plateau of the coefficients of the spectral representation [5]. Convergence tests are performed for $N={10}^2,{\mathrm{\hspace{0.33em}10}}^3$ and a larger range of variances, ${\sigma }^2=0.1,\mathrm{\hspace{0.33em}4},\mathrm{\hspace{0.33em}10}$.

Validations tests are finally conducted by comparing mean values and variances of the velocity components inferred from ensembles of 100 realizations of the solution of the homogeneous boundary value problem (1-4) for ${\sigma }^2=0.1,\mathrm{\hspace{0.33em}0.5},\mathrm{\hspace{0.33em}1},\mathrm{\hspace{0.33em}1.5},\mathrm{\hspace{0.33em}2}$ and fixed $N={10}^2$ with theoretical and numerical results given in the literature.

The formulation of the benchmark problems, data files, numerical codes and functions are given in a Git repository [4] which can help the interested readers to test their numerical methods for flows in highly heterogeneous aquifers.

In this section, we describe the methods on which the codes participating in the present benchmark study are based. For the FDM, FEM, and DGM approaches, based on well established finite difference/element methods, only a few aspects relevant for the construction of the numerical scheme will be summarized in the following. The spectral approaches CSM and GSM adapted to this benchmark, as well as the GRW algorithm using computational particles proposed in this article will be presented in more detail.

Ensembles of numerical solutions of the flow problem are often used for the purpose of validation of the numerical schemes through comparisons with first-order approximation results (e.g., [7, 22, 68, 71]). The latter are obtained by perturbation expansions truncated at the first order in ${\sigma }^2$, equivalent to a linearization of the flow problem (1-4) (e.g., [19, 31]). Here, we consider first-order approximations of the velocity field numerically generated with a Kraichnan procedure, based on proportionality of the spectral densities of the random velocity and that of the $\ln (K)$ field [78, Appendix A]. The number of random periodic modes used in the Kraichnan generators of the hydraulic conductivity and of the velocity field (in the case of the linear approximation) is set to $N=100$. Validation tests of the five schemes are performed for $\ln (K)$ fields with Gaussian and exponential correlations and a fixed small variance ${\sigma }^2=0.1$. The deviation from the predictions of the linear approximation is further investigated with the FDM, FEM, DGM, and GRW schemes for increasing ${\sigma }^2$ of 0.5, 1, 1.5, and 2. The number of realizations, in all cases investigated here, is fixed to $R=100$. Numerical investigations indicate that the chosen values of the parameters $N$ and $R$ ensure unbiased statistical estimations of the dispersion of a passive scalar transported over tens of correlation lengths in heterogeneous aquifers characterized by a small variance ${\sigma }^2=0.1$ of the $\ln (K)$ field [27].

The homogenous problem (1-4) is solved in the same two-dimensional domain, with $L_x=20$ and $L_y=10$, as that used for code verification tests in the previous section. Graphical illustrations and analysis of numerical solutions are given in Supplementary materials S1 to S4. The ensemble of realizations were computed for space steps $\mathrm{\Delta }x=\mathrm{\Delta }y=0.02$, with the FDM scheme and, to reduce the computational effort, for larger steps $\mathrm{\Delta }x=\mathrm{\Delta }y=0.05$, with the FEM and DGM schemes, and $\mathrm{\Delta }x=\mathrm{\Delta }y=0.1$, with the GRW scheme. Constrained by the maximum Matlab array size allowed on the computer used to obtain CSM solutions, we use a limited number of 160 Chebyshev collocation points in both $x$- and $y$-directions, even though a larger number of points in $x$-direction is required to obtain accurate solutions [see 3, Fig. 9]. Sample solutions obtained with FDM are illustrated in Supplementary material S1. FEM, DGM, and GRW solutions, not shown in S1, have a similar appearance. Instead, among the realizations of the CSM solution for exponential correlation of the $\ln (K)$ field there also are unacceptable solutions with extremely large and irregular variations. Therefore, we selected an ensemble of 100 valid realizations by disregarding the “bad” solutions after the visual inspection of the larger initial ensemble of solutions, as illustrated in Supplementary material S2. The issue of “bad realizations” is common in situations where, due to computational constraints, numerical simulations have to be conducted on relatively coarse grids (we already encountered this issue during the preparation of the transport simulations presented in [57, 76]).

Monte Carlo inferences of the mean and variance are obtained by computing averages over the ensembles of $R$ realizations followed by spatial averages [75, Appendix B3]. Error bounds of the ensemble-space averages are given by standard deviations of the ensemble averages at grid points, estimated by spatial averages. The spatial averaging domain is an inner region of the computational domain unaffected by the boundary conditions approximately delimited by visual inspection of the variances estimated on longitudinal and transverse sections through the computation domain presented in Supplementary materials S3 and S4 [see also 7, 42].

Theoretical investigations indicate that the variance of the hydraulic head is much smaller than that of the hydraulic conductivity. Spectral perturbation approaches for flow in two- and three-dimensions predict variances of the order ${\sigma }_h^2\sim {(\sigma \lambda J)}^2$ [1, 6, 52]. The estimated variances presented in Table 6 show an acceptable agreement with the spectral theory, which for our numerical parameters predict ${\sigma }_h^2\sim {10}^{-4}$.

The first-order approximate relations between the variance of the velocity components, ${\sigma }_{V_x}^2$ and ${\sigma }_{V_y}^2$, and the variance ${\sigma }^2$ of the $\ln (K)$ field are given by ${\sigma }_{V_x}^2=\frac{3}{8}{\sigma }^{2}U$ and ${\sigma }_{V_y}^2=\frac{1}{8}{\sigma }^{2}U$, respectively [19, Eq. (3.7.27)]. Note that these relations do not depend on the shape of the correlation function of the $\ln (K)$ field. Since, as mentioned in Section 2, for the purpose of statistical inferences we use dimensionless velocities, $U=1$. Then, with the smallest variance of the $\ln (K)$ field considered in this study, ${\sigma }^2=0.1$, the velocity variances predicted by the linear theory are ${\sigma }_{V_x}^2=3.75\cdot {10}^{-2}$ and ${\sigma }_{V_y}^2=1.25\cdot {10}^{-2}$, respectively.

The estimated mean values and variances of the velocity components are presented in Tables 7 and 8. FDM, FEM, DGM, and GRW results are generally close to theoretical first-order predictions within a range of less than 10% in case of FDM, FEM, DGM, and GRW approaches, excepting FEM estimated variance of the transverse velocity in the Gaussian case which deviates with about 30%. CSM results show deviations of about 30%, for the variance of the transverse velocity in the Gaussian case, and deviations even larger than 100% for both velocity variances, in the exponential case. These differences can be again attributed to the limited number of collocation points used in the CSM approach.

The results obtained with the FDM, FEM, DGM, and GRW approaches for larger ${\sigma }^2$ are shown in Figs. 11 and 11. We remark that the velocity variances progressively increase above the trend predicted by the linear theory, similarly to the Monte Carlo results presented in the literature (e.g., [7, 42, 68]). The differences between the results obtained with the three methods remain within the estimated error bounds determined by the present numerical setup (number of realizations, resolution, and dimension of the computational domain). For ${\sigma }^2\le 0.5$ the estimated variances practically coincide with the first-order predictions. Within this range, velocity fields can be quite well approximated by the linearized solution of the flow problem and simulations of solute transport in groundwater can be carried out with Kraichnan generated velocity fields at lower computational costs [75, 73, 78].

We conclude this section with a remark about the GRW results. The ensembles of GRW simulations were computed for ${10}^6$ iterations in all cases, excepting the Gaussian case with ${\sigma }^2=2$, where the number of iterations was set to $5\cdot {10}^6$. As indicated by Figs. 9 and 9, as well as by Supplementary materials S3 and S4, these numbers of iterations are not large enough to ensure the convergence of the solution for homogeneous problems. Surprisingly, the variance estimates from Tables 7-8 deviate with maximum 10% from the predictions of the first order theory for small ${\sigma }^2$. For larger ${\sigma }^2$, the GRW results are also close to those obtained with the FDM, FEM, and DGM schemes for ${\sigma }^2\le 2$ in Gaussian case (Fig. 11) and for ${\sigma }^2\le 0.5$ in exponential case (Fig. 11).

The numerical approaches tested in the present benchmark study passed, at least partially, the validation test consisting of comparisons of the numerical solutions obtained for ${\sigma }^2=0.1$ and $N=100$ with theoretical results provided by first-order perturbation approaches. For larger variances, up to ${\sigma }^2=2$, the statistical inferences performed by averaging over Monte Carlo ensembles of solutions given by the FDM, FEM, DGM, and GRW schemes were close to those obtained by simulations with similar numerical setup reported in the literature. However, these encouraging results are not a guarantee for the accuracy of the numerical method in case of higher values of ${\sigma }^2$ and $N$, needed to solve the flow problem for highly heterogeneous aquifers and large spatial scales. This is indicated, for instance, by the situation of the GRW ensembles of solutions for Gaussian correlation of the $\ln (K)$ field which provide satisfactory statistical estimates even though they are not yet convergent (see Section 5) and cannot be used for single-realization approaches in subsurface hydrology.

For Gaussian correlation of the $\ln (K)$ field the code verification tests were accurately solved, with the exception of the FDM solution for $N={10}^3$ and ${\sigma }^2=6$ and that of the GRW solutions for $N={10}^3$ and ${\sigma }^2\ge 2$. Instead, in case of exponential correlation, for all the schemes, excepting the CSM scheme, there are ranges of parameters for which accurate solutions are not realizable: $(N\ge {10}^2,{\sigma }^2\ge 1)$ for FDM, $(N\ge {10}^2,{\sigma }^2\ge 2)$ for FEM, $(N={10}^2,{\sigma }^2\ge 4)$ and $(N={10}^3,{\sigma }^2\ge 2)$ for DGM, $(N={10}^2,{\sigma }^2\ge 2)$ and $(N={10}^3,{\sigma }^2\ge 1)$ for GRW. The limited solvability of the flow problem in case of exponential correlation found in our benchmark study, for a resolution of the numerical solution which is one order of magnitude higher than usually recommended [see e.g., 1, 22], completes similar results presented in the literature. For instance, [12] report difficulties to obtain the convergence of their numerical schemes for exponential correlation of the $\ln (K)$ field with ${\sigma }^2=2$ and [33] show that for high resolution of the exponentially correlated $\ln (K)$ field with ${\sigma }^2>2$ the required discretization of the hydraulic head solution exceeds the limit of the available computing resources.

The observed large errors produced by the three related discretization schemes, FDM, FEM, and DGM, in the test scenarios with exponential correlation and large values of the parameters $N$ and ${\sigma }^2$ (cf. Tables 3 and 4) may be explained as follows: The listed errors consist of the discretization error of the respective numerical scheme and the error of the linear solver. In practice, the latter one cannot be computed, however, the observed high condition numbers and the computed residuals indicate that the error of the linear solver grows when using the above-mentioned parameters (in fact, the product of condition number and residual is an upper bound for the arithmetic error [37]). More precisely, we found that the condition number increases with ${\sigma }^2$ in a comparable way for both Gaussian and exponential correlation models and for both numbers of modes. The residuals remain almost unchanged in the Gaussian case, but in exponential case, the increase of ${\sigma }^2$ and $N$ cause the increase of the residuals and, consequently, the increase of the error bound [3, Sect. 9]. Regarding the linear solver: for the FDM and DGM schemes we used Matlab’s direct solver mldivide, which for our system matrices, calls UMFPACK’s Unsymmetric MultiFrontal method [20], and for the FEniCS implementation of the FEM scheme we used the default solver in FEniCS, which calls UMFPACK as well. We also tested various other solvers (e.g., UMFPACK’s SPQR [21] and ILUPACK’s AMG_solver [10]), however, a smaller residual was not producible.

In exponential cases of large $N$ and ${\sigma }^2$, the coefficients $K$ contain large numbers of periodic modes with amplitudes spanning over several orders of magnitude. The high resolution necessary to prevent the occurrence of ill-conditioned matrices and large residuals leads to extremely large degrees of freedom of the linear system of equations which render impractical, from the point of view of computational resources, the use of direct flow solvers such as FDM, FEM, and DGM considered in this study. A possible remedy could be provided by numerical upscaling techniques developed for problems with rough coefficients, which capture the effect of small scales on large scales without resolving the small scale details [13, 36, 39].

CSM schemes verify the manufactured analytical solutions with extremely small errors in both one- and two-dimensional cases. This high precision is rendered possible by the use of the analytical expressions of the derivatives of the $K$ field given by the Kraichnan algorithm. The convergence is demonstrated by the decrease of Chebyshev coefficients towards the roundoff plateau. The non-homogeneous flow problem can be solved with relatively small numbers of collocation points (less than 100 points on both spatial directions), which are enough for the spectral representation of the simple, smooth manufactured solutions chosen for the purpose of code verification. But in case of highly oscillating solutions of the homogeneous problem for exponential correlation of the $\ln (K)$ field, a much larger number of points is required, which increases the dimension of the problem beyond the limit of the available memory. Due to this constraint, the maximum number of collocation points in each direction has been fixed to 160, which did not ensure reliable solutions for all the realizations of $\ln (K)$.

The GSM scheme has been proposed to avoid the use of the analytical derivative of the hydraulic conductivity $K$ and has been used only in solving the one-dimensional non-homogeneous problem for comparisons with manufactured solutions. The errors are again very small in the Gaussian case, but extremely large in the exponential case. Some tests made in the attempt to improve the precision indicate that exceedingly large computing resources would be necessary to render the problem computationally feasible (see Table 5).

The GRW scheme is a transitory scheme using computational particles, which is accurate and unconditionally stable. The drawback is the large number of time iterations needed to reach the stationary state corresponding to the steady-state solution of the flow problem. In the Gaussian case, approximations with satisfactory accuracy of the manufactured solution can be obtained with a coarse discretization in reasonable computing times from minutes to hours. In the exponential case, a much finer discretization is needed and the computing time increases from hours to days.

The numerical schemes tested so far on benchmark problems have their own strengths which make them useful in specific applications, depending on the objective pursued. FDM schemes are appropriate for fast computations of Monte Carlo ensembles in case of moderate aquifer heterogeneity, FEM and DGM schemes are appropriate to cope with higher heterogeneity in either single-realization or Monte Carlo approaches, GRW is rather appropriate for single-realization flow problems with heterogeneous coefficients, and CSM/GSM approaches could be useful in solving problems with heterogeneous coefficients and source terms given by smooth functions. We have seen that all the schemes tested in this study face the issue of computational feasibility in case of an exponential structure of the log-conductivity field. It is therefore recommendable that, unless well documented by experimental data for a specific problem, the exponential correlation model should be avoided. Instead, Gaussian models can be used to account for effective dispersion associated with an integral scale [19, 79], as well as in flow and transport simulations for fractal aquifer structures [e.g., 25].

The authors are grateful to Dr. Maria Crăciun, for her help in testing the random field generators, and to Dr. Emil Cătinaş for carefully reading the manuscript. Andreas Rupp acknowledges the financial support of the Deutsche Forschungsgemeinschaft (DFG) – Germany under Grants RU 2179 “MAD Soil - Microaggregates: Formation and turnover of the structural building blocks of soils” and EXC 2181 “STRUCTURES: A unifying approach to emergent phenomena in the physical world, mathematics, and complex data”. The work of Nicolae Suciu was supported by the DFG Grant SU 415/4-1 “Integrated global random walk model for reactive transport in groundwater adapted to measurement spatio-temporal scales”.