Shile et al., [1] propose a radial basis function (RBF) collocation method to simulate steady-state flows in heterogeneous aquifers. The RBF approach is tested on benchmark problems proposed in [2], using the Matlab codes published in [3]. The hydraulic-head solution is implemented with the codes posted at [4]. Note that these codes are not properly cited, the reference [15] cited by Shile et al., [1], containing a sample Matlab code, does not provide a link to the codes and functions used to solve the RBF interpolation problems presented in their paper. Also, the reference [3] is not cited as a distinct publication, a link to the Git repository being instead added to reference [2], so that the reader cannot know where different elements of the benchmark were taken from.

The benchmark study of Alecsa et al., 2020a [2] is two-fold. First, code verification tests, with manufactured solutions, were conducted for a finite difference method (FDM), a global random walk algorithm (GRW), a finite element method (FEM), and a discontinuous Galerkin method (DGM), as well as for a collocation spectral method (CSM). The realizations of the hydraulic conductivity were provided by a randomization method as analytic space functions. The manufactured solutions were constructed for a fixed realization of the hydraulic conductivity as smooth space functions. These two ingredients favor the CSM, which, using analytical differentiation of the hydraulic conductivity and relatively small numbers of collocations points, yields solutions with an accuracy that is orders of magnitude higher than that of the classical discretization approaches FDM, FEM, GRW, and DGM.

The second part of the study presented in [2] consists of statistical inferences on ensembles of Monte Carlo simulations of transport computed for increasing variance ${\sigma }^2$ of the log-hydraulic conductivity. This setup, representative for realistic situations, favors the classical discretization methods. Instead, the feasibility of the CSM is limited to small values ${\sigma }^2$, due to the need to increase the number of collocation points as the heterogeneity of the hydraulic conductivity increases.

The RBF approach proposed by Shile et al., [1] passes the code verification tests, with errors with respect to the manufactured solutions smaller than those of the classical discretization methods and much larger than the errors of the CSM method presented in [2]. Instead, the statistical inferences results show striking discrepancies with results obtained within the same numerical setup. As shown in the following, the results from [1] cannot be explained by numerical artifacts and there is a strong evidence that the Monte Carlo results are fabricated and the comparisons with previous studies are falsified.

The numerical procedures used in Monte Carlo simulations of transport in groundwater are validated through comparisons with linear theory and similar numerical approaches (see [5, 6, 7, 8, 2]). One of the mostly used validation test in two-dimensional simulations consists of comparisons of the variances of the longitudinal velocity $V_x$ and of the transverse velocity $V_y$, defined by

where $⟨\cdot ⟩$ stands for the average over the ensemble of realizations.

The studies cited above indicate a similar, monotonic increase of the velocity variances with the increase of the variance ${\sigma }^2$ of the log-hydraulic conductivity above the variances predicted by the linear theory, which, for unit mean longitudinal velocity considered in [2, 1], are given by

The departure from the linear approximation can be illustrated by the semi-analytical second-order approximation, calibrated by Monte Carlo simulations, obtained by Gotovac et al., [8]:

The percentage deviation from the linear approximation (2) will be evaluated in the following by

In Table 1 the results of Shile et al., [1, Table 9] for ${\sigma }^2=1.5$ and exponential correlation of the log-hydraulic conductivity are compared via (4) to the DGM results of Alecsa et al., 2020a [2], to the results from [5, 6, 7], which are close at less than 1.5% to each other [7], and to the second-order approximation (3). The deviation from linear theory of the RBF results sown in Table 1 is one order of magnitude smaller than in earlier published studies.

The same striking discrepancy is shown in Table 2 for the results of Shile et al., [1] for the whole range ${\sigma }^2$ and for both Gaussian and exponential correlations. Tables 3 and 4 show deviations from the linear theory of the DGM results [2] and of the second-order approximation (3), respectively, which are again one order of magnitude larger that those of Shile et al., [1]. And, finally, Table 5 shows that the RBF variances results [1] are systematically smaller than those given by the second-order approximation (3) and the discrepancy increases with ${\sigma }^2$.

A visual comparison of the results of Shile et al., [1] and earlier results from Alecsa et al., 2020a [2] and Gotovac et al., [8], computed with Eq. (3), is presented in Figs. 2 - 4. Here again one can see that the results presented by Shile et al., [1, Table 9] largely underestimate the departure of the velocity variances from the prediction of the linear theory.

The second striking discrepancy is the non-monotonous behavior of the standard deviations in Shile et al., [1, Table 9], indicated by error bars in Fig. 2 and Fig. 4, as compared with monotonically increasing standard deviations in Fig. 2 and Fig. 4 which present results from [2], that are similar to those published by Bellin et al., [5, Figs. 6a and 6b].

The standard deviations of the ensemble average velocity variances (1) are computed by space averages according to

where ${⟨\cdot ⟩}_s$ stands for the spatial average over the inner core region of the computational domain not affected by the boundary conditions [2]. If the random velocity is statistically homogeneous ${\sigma }_{V_x}^2$ and ${\sigma }_{V_y}^2$ do not depend on the spatial variable and their standard deviations defined by (5) vanish. However, since the numerically generated fields are not strictly stationary, the ensemble averages vary from point to point, resulting in a consistent increase of the standard deviations (5) with the increase of ${\sigma }^2$, reported in many studies published since the pioneering article of Bellin et al., [5].

Shile et al., [1] used the codes from [3] to generate the random hydraulic conductivity and claim [1, Sect. 4.4] that they used the same numerical setup as Alecsa et al., 2020a [2] and the definition (5) to compute standard deviations of spatially variable velocity variances. Moreover, the code verification tests with two-dimensional manufactured solutions [1, Tables 5 and 6] show that, for ${\sigma }^2\le 2$, their RBF solver is at least as accurate as the DGM code used in [2, Table 1 and 3]. In these conditions the results obtained with the two approaches should be, in principle, almost identical.

On the other side, the manufactured solutions used in code verification tests are very smooth, with slow oscillations in space (see [1, Fig. 12]), and can be approximated by relatively small numbers of collocation points used in the RBF code, as well as in the CSM code tested in [2]. In these circumstances, the two collocation methods mentioned above benefit from the explicit analytical expression of the hydraulic conductivity used in the benchmark problems [2, Eq. (A3)] to compute its derivatives. Therefore, they outperform the classical approaches which discretize the derivatives (see discussion in [2, Sect. 6]). Instead, the solutions of the homogeneous flow equation used in Monte Carlo simulations are highly heterogeneous and in a first-order approximation follow the shape of the hydraulic conductivity shown in [1, Fig. 2]. Accurate representations of such solutions require more collocation points than in case of the smooth manufactured solutions.

We have thus two possible situations. First, if the same number of collocation points as in code verification tests allows accurate representations of the solutions, as the authors claim in [1, Sect. 4.4], then we are in the case described above, when the RBF and DGM results should be almost identical. This is not the case. The second possibility is that the number of collocation points is not large enough. Then, the representation of the solution is inaccurate and there are realizations with large deviations from the expected almost-laminar behavior (see Supplementary material 2. in [2] which shows “bad” realizations of the CSM collocation method). In this case, larger spatial variability yields larger standard deviations of the velocity variances than the DGM code. This also is not the case (compare Figs. 2 and 2, for Gaussian correlation, and Figs. 4 and 4, for exponential correlation of the $\ln K$ field). A first conclusion of this analysis is that the standard deviations presented by Shile et al., [1, Table 9] cannot be produced by their RBF code.

To elucidate the issue of non-monotonous standard deviations shown in Fig. 2 and Fig. 4, I computed ensembles of 100 realizations of the $\ln K$ field, for Gaussian and exponential correlations, by using the same parameters and generated with the same codes as in [1, 2], posted in the new folder “Comment on [Shile et al., 2024]” of the Git repository [3]. The estimated mean and standard deviation of variance ${\sigma }_{\ln K}^2$ are plotted in Fig. 6. As expected, the standard deviation $SD({\sigma }_{\ln K}^2)$ increases monotonically with the nominal parameter ${\sigma }^2$. With these, the nominal variance ${\sigma }^2$ used in the Monte Carlo simulations is approximated with a $95\%$ confidence level by

From (6) and (2) we get the following linear approximation of the velocity variances,

If number of random modes used in the field generator and the number of realizations go to infinity, one approaches a statistical homogeneous $\ln K$ field, $SD({\sigma }_{\ln K}^2)\to 0$ and (7) reduces to (2).

The quantities $\frac{3}{8}SD({\sigma }_{\ln K}^2)$ and $\frac{1}{8}SD({\sigma }_{\ln K}^2)$ in (7) are confidence half-intervals of the mean variances $\frac{3}{8}\overline{{\sigma }_{\ln K}^2}$ and $\frac{1}{8}\overline{{\sigma }_{\ln K}^2}$ induced by the statistical inhomogeneity of the $\ln K$ field. The same repository also contains ensembles of linear approximations of the corresponding velocity components generated by the Kraichnan method described in [9, Appendix C.3], without the need to solve the flow equation. The comparison presented in Table 6 shows that the confidence intervals from (7) estimate the standard deviations of the velocity variances in a first-order approximation. Since numerical methods to solve the flow equation are prone to errors which increase the statistical inhomogeneity, the estimates from Table 6, computed without solving the flow equation, represent lower bounds of the standard deviations.

The detailed comparison of standard deviations for exponential correlation from Fig. 6 shows a total mismatch between the RBF and DGM results. Moreover, the RBF data points show non-monotonous dependence of the standard deviations on ${\sigma }^2$ and values smaller than the lower bounds $\frac{3}{8}SD({\sigma }_{\ln K}^2)$ and $\frac{1}{8}SD({\sigma }_{\ln K}^2)$. Since such results cannot be obtained with any numerical method, which was also validated by comparisons with manufactured solutions, the only possible explanation for these flagrant discrepancies of the results from Shile et al., [1, Table 9] is that they are made up.

Summarizing, The results presented by Shile et al., [1] mimic the general feature that the distance between the numerically estimated variance of the velocity components and their linear approximation increases with the increase of ${\sigma }^2$, while neglecting three other aspects:

a. the magnitude of the distance from the linear approximation, which is almost the same in previously published studies;

b. the monotonous increase of the standard deviations of the velocity variances;

c. the lower bound of these standard deviations given by the standard deviation of the variance of the $\ln K$ field multiplied by 3/8 and 1/8, for longitudinal and transverse velocity, respectively.

Based on a comparison plotted in their Fig. 13, Shile et al., [1] conclude the presentation of the statistical inferences by stating that “there is good agreement between all the stochastic solutions, supporting the verification of RBFPUM-QR with previously published stochastic results”. However, as we’ll see in the following, this statement is not supported by the numerical results presented by the authors in [1, Table 9].

The analysis of the results shown in Figs. 8 and 8 demonstrates that the claim “Fig. 13 confirms that for higher values of ${\sigma }^2$, RBFPUM-QR produces results similar to those obtained using the Eqs. (42) and (43)” is false.

First, I compare in Fig. 8 the data points given by Eq. (3) and the RBF data from [1, Table 9] for ${\sigma }^2\le 2$ to the data points from [1, Fig. 13] for ${\sigma }^2=1$ and ${\sigma }^2=2$. The distance of 0.123 between the points corresponding to RBF and to Eq. (3) on the lower branch for ${\sigma }^2=2$ corresponds to about half the height of the diamond marker in [1, Fig. 13]. As such, this difference should have been visible in [1, Fig. 13] and the marker should have been placed just below the line representing Eq. (43).

The mismatch between the results presented in Shile et al., [1, Table 9 and Fig. 13] is even more obvious in Fig. 8, where the Eq. (3) and the RBF data points from Fig. 8 are completed with the points computed by (3) on which I superposed the RBF points, as in [1, Fig. 13]. The double data points on the upper an lower branches are far away from the trend of RBF points from Fig. 8. This demonstrates that the RBF points cannot be perfectly aligned with the points given by the second-order approximation (3). Hence [1, Fig. 13] falsifies the comparison between the RBF and the previously published results.

The article of Shile et al., [1] reproduces, with unnecessary details, the description of the benchmark for flow in saturated aquifers proposed by Alecsa et al., 2020a [2] and the principle of the RBF interpolation approach. The codes from [3] and [4] are further used to conduct code verification tests of the RBF solution. Such tests constitute a necessary preparatory step for applications of the numerical methods to more complex and realistic problems.

The RBF method performs better than classical discretization methods and could be a promising alternative approach in numerical simulations of groundwater flow based on Kraichnan generators of hydraulic conductivity which provide differentiable coefficients of the flow equation. The challenge is to prove that the RBF code can manage realistic spatial heterogeneity conditions which, unlike the smooth manufactured solutions used in code verification tests, require significantly larger numbers of collocation points. Unfortunately, the results of the Monte Carlo simulations presented by Shile et al., [1] are fabricated and the comparison with similar studies, which is supposed to validate the RBF approach, is falsified.