Comment on “Asymptotic dispersion in 2D heterogeneous porous media determined by parallel numerical simulations” by J. R. de Dreuzy et al.

[1] de Dreuzy et al. [2007] present a numerical study aimed at determining the asymptotic dispersion coefficient DLA for flow (uniform in the mean) through porous formations of a two-dimensional structure, for a broad range of heterogeneity levels. The numerical analysis appears to be very accurate, and by and large the results are of definite interest to the stochastic subsurface community. We believe that these are among the most reliable results concerning longitudinal macrodispersivity for flow in strongly heterogeneous formations. [3] For f(κ) lognormal, the log conductivity field underlying the solution (1) still differs from the multi-Gaussian random field modeled by de Dreuzy et al. [2007]. This manifests in the higher-order statistics of Y = ln K, i.e., the multipoint joint pdf, which is different from the multi-Gaussian field, and this may have a significant impact on the derived macrodispersivity for highly heterogeneous formations (σ2 ≫ 1). It is nevertheless instructive to compare (1) with the numerical results of de Dreuzy et al. [2007] in order to assess the applicability of alternative, simple nonlinear approaches to transport in heterogeneous formations. [4] Figure 1 displays DLA as function of σ2 as obtained through formula (1) (thick solid line) and by the numerical simulations of de Dreuzy et al. [2007] (dots); the gray line displays the first-order solution. Along the assumptions of de Dreuzy et al. [2007], a lognormal f(κ) was plugged in (1). The agreement between (1) and the numerical results is quite good for a broad range of σ2, up to the very large value of σ2 = 6.25. The latter characterizes strongly heterogeneous formation. We believe that this result is quite promising, as it compares for the first time a physically simple model with accurate numerical simulations for such large values of the log conductivity variance. [5] The only simulation for which the agreement deteriorates is σ2 = 9, for which the solution (1) overestimates DLA by a factor of two. There are a few possible causes for this disagreement. Thus, one of them was discussed above; that is, the two conductivity structures (numerical and analytical) differ at high order in σ2. Secondly, it was shown in previous work [Janković et al., 2003] that the effective medium approximation works less well in two-dimensional transport than in 3D fields (this may be simply explained by observing that in 3D each element is surrounded by a larger number of neighboring elements than in 2D, rendering the self-consistent argument more plausible). [6] Besides problems concerning the mathematical approach, the derivation of DLA for high σ2 may be strongly affected by a few features of the numerical procedure. We have encountered such problems when determining numerically the travel time distribution of solute particles in highly heterogeneous media [Fiori et al., 2006]. It is found that it is generally very difficult to get stable and reliable DLA for very large σ2 and the problems encountered by de Dreuzy et al. [2007] toward obtaining a stable dispersivity value for σ2 = 9 case confirm our previous findings. [7] Thus, in our previous work [see, e.g., Dagan et al., 2003] we found that presence of zones characterized by very low values of K, i.e., by very low velocity, are primarily determining the value of DLA for large σ2. In fact, the retention of solute particles in those areas and the large residence time, lead to a dramatic increase of the dispersivity. In contrast, preferential, high-velocity channels also cause an increase of DLA, but at a much lesser extent. There are quite a few low-K areas for the lognormal f(κ) and they are isolated. Consequently, the domain and plume areas should be very large in order to adequately sample the K field at the lower extreme end, when σ2 is large. Achieving ergodicity can be quite difficult as convergence of relevant quantities, e.g., high-order moments, may oscillate, with "jumps" and discontinuities [see, e.g., Fiori et al., 2006, Figure 6]. [8] Even if the domain size is extremely large, the low-K zones may not be sampled by the plume for other reasons, e.g., insufficient number of particles, type of average used for the interblock conductivity and errors of the particle tracking algorithm. The combined effect of one or several of the above effects may be lumped for simplicity into a generic "numerical dispersion" effect, which filters out the influence of K values below a certain cutoff. Although extremely small, the latter has a dramatic effect on DLA when σ2 becomes very large, as correctly pointed out by the authors when discussing the impact of local diffusion on dispersivity. [9] In order to mimic the presence of a possible cutoff in K we have introduced in the past such cutoffs in (1) [Dagan et al., 2003; Fiori et al., 2006]. A simple argument brought up by Dagan et al. [2003] relates the cutoff value κ = κc to a local diffusion process, with the approximate relationship Pec ≈ 1/κc, where Pec = uλ/D0 is a cutoff-related Peclet number (D0 is the coefficient of diffusion). The idea is that below κc = Kc/KG the function XD(κ) in the integral (1), which is related to the residence time in the low-K element, does not grow anymore and it levels off at the constant value XD(κc). This reflects the inability, for one or more of the above reasons, of the numerical setup to capture the spatial variability of K below a certain value. [10] For the sake of illustration, we represent in Figure 1 the asymptotic longitudinal dispersivity as function of σ2 on the basis of (1) for a few values of the cutoff Yc = ln Kc = −9.0 ÷ −7.0. Although very small, these cutoff values have a significant impact on DLA for high σ2, as observed in Figure 1. A good agreement between the numerical results of de Dreuzy et al. [2007] and formula (1) is obtained in the entire range of σ2 values for Yc ≈ −9.0. If the above arguments hold true, while other explanations are still possible, this would mean that an "equivalent" Peclet number of the order of Pec ≈ 104 may affect the presumably pure advection numerical simulation of de Dreuzy et al. [2007]. [11] Summarizing, we found that our simple self-consistent solution (1) is in satisfactory agreement with the numerical simulations, all the approximations/limitations notwithstanding. The derivation of a reliable numerical solution of transport in strongly heterogeneous formations by de Dreuzy et al. [2007] constitutes in our view a valuable contribution to the theoretical body of knowledge, for which they are commended. In this comment we wanted to raise a few problems, of both numerical and analytical nature, which are encountered when dealing with aquifers of very high σ2. More work is needed in order to fully understand the transport dynamics in such complex systems and its relation with the underlying conductivity structure. This is even more so for the realistic three-dimensional heterogeneity, as encountered in aquifer applications.