[1] Eliassi and Glass [2001] analyze the continuum-scale modeling of gravity-driven fingers in two-dimensional (2-D) unsaturated flows in porous media. They report that the Richards equation does not contain the essential physics required to model gravity-driven fingers, at least not when using standard monotonic constitutive equations and hysteretic relations. Eliassi and Glass [2001] have conducted an examination of the modeling of the gravity-driven finger phenomenon and have verified their results by testing using several numerical schemes and packages. The purpose of this comment is to point out that (a) the basis of their assumptions on monotonicity are not fully met, (b) the temporal truncation errors lead naturally to oscillation in the solutions, and (c) these oscillations pose difficulties in the choice of switch value in the hysteresis relations. [2] Eliassi and Glass [2001] cite the papers by Rubin and Steinhardt [1963], Rubin et al. [1964], and Youngs [1995] to support the proposition that the Richards equation should admit only monotonic solutions. Youngs [1995] makes no comment on the monotonicity properties of the solutions of the Richards equation. Rubin and Steinhardt [1963] do address the monotonicity property and draw heavily on the theoretical results obtained by Philip [1957a, 1957b]. These theoretical results are based on assumptions of 1-D vertical flow and nonhysteretic behavior. Philip [1957b] discusses at some length the hysteretic properties of a propagating water front and gives conditions such that the hysteresis is not activated. These general conditions are satisfied in the numerical experiments of Eliassi and Glass [2001], at least on the centerline or line of symmetry of the finger (see their Figure 4a). [3] However, Eliassi and Glass's [2001] numerical solutions are calculated in a properly 2-D region, are not to the same level of theoretical proof as the results of Philip [1957a, 1957b], and hence are open to question. They could consider the pseudo-1-D case where the applied flux at the top of the region [see Eliassi and Glass, 2001, Figure 2] is uniform from η = 0 to η = L*. If their centerline results are repeated across the full range of η, then they have supported their suggestions about monotonicity. [4] In two dimensions, there does not seem to be any specific literature relating to monotonic solutions of the Richards equation or its linearized form. The concept of monotonicity is normally applied to the study of functions of one independent variable. Philip [1957b] clearly states that “…as presented here, the theorems are established only for one-dimensional regions…” (p. 437). In two or more dimensions a maximum principle can be obtained bounding the solutions of linear parabolic partial differential equations, where Dirichlet boundary conditions are applied [Egorov and Shubin, 1992, pp. 163–166]. Indeed, it can be shown that the solution attains its maximum (and minimum) value on the bounding region. In Figure 1 an arbitrary region is shown for a linear parabolic equation subject to Dirichlet conditions. A maximum (and minimum) of the condition imposed on the boundary is shown, along with schematic solution contours. [5] In two dimensions a direction or curve needs to be specified so that dimensional concepts of monotonicity can be applied along the direction or curve. The solution may be monotonic on sections of the bounding curve. [6] Eliassi and Glass [2001] do seem to have attempted to apply a result, which holds on one-space dimension, to a two-space dimension problem. Such problems may give solutions which are monotonic in time; that is, a direction is specified. The question of 2-D spatial monotonicity is yet to be properly formulated. [8] Eliassi and Glass [2001] use a small initial time step of 10−10 but allow the time step to alter as the solution progresses. The maximum possible time step, τmax, is given by Δτmax = 0.1(Δ/vmax), where vmax is the dimensionless maximum Darcy velocity, and Δ = Δη = Δξ is the spatial grid spacing. The initial time step is certainly much smaller than the spatial time steps used in the numerical solutions. However, there is no detail reported in the paper on the history of the time steps used in the solution or of the magnitude of the Darcy velocity. Given the length of the computation, the reader is left to conjecture that the time and space steps may eventually be of the same order, with the time step near its maximum constraint. The possibility that Eliassi and Glass's results [2001, Figures 1 and 3] for a reduced spatial step size may be due to a reduction in the temporal step size is not yet ruled out. [13] Eliassi and Glass [2001] provide graphic output from their computer model, showing the development of a finger under a variety of conditions (see their Figure 3; also note additional oscillations behind the front and up through the center of the plume in their Figures 9a and 9b). These oscillations are located on the region inside the finger, and these oscillations are probably traveling or standing wave solutions arising from the temporal truncation error leading to the telegraph equation. Along the center of the plume the water content does not change markedly, and the assumption of approximately constant Γ and constant κ will apply. [14] Figure 9 of Eliassi and Glass [2001] shows the wave-like structure of the local truncation error (LTE) in the tip of the plume. Eliassi and Glass [2001, Figure 9a] show oscillations in LTE in the η direction, while their Figure 9b shows a more complex set of oscillations in both the η and ξ spatial directions, this being where hysteresis was not incorporated. These oscillations are localized in the tip of the finger and bounded by the steep solution gradients near the edges. Certainly the horizontal (η coordinate) oscillations are standing or traveling waves. [18] Eliassi and Glass [2001] use essentially the same switching mechanism to handle the hysteretic change of phase. They set δ to a small number slightly larger than the computer's round off error and suggest that δ ≈ 10−15, which is definitely down in the numerical noise region, given the accuracy and convergence criteria used. The results in Eliassi and Glass's [2001] Figure 9a (with hysteresis) and Figure 9b (without hysteresis) suggest that hysteresis is damping some of the oscillation in the calculation. Eliassi and Glass [2001, section 3.5] also suggest that the plume in their Figure 3c shows numerical oscillations, these presumably arising from the error near the front. However, Eliassi and Glass's [2001] Figure 3c also shows some oscillations near the top boundary, and these are not negligible compared to the size of the hysteretic switch. Eliassi and Glass's [2001] figures show water content values, but the level of resolution is not at the machine rounding level of 10−15. The results of Eliassi and Glass [2001] do show differences between the solutions for cases where hysteresis is (or is not) included. This leads to the question as to the size of the switch mechanism to be used in such models and whether it is to be based on a numerical artifact or on some underlying physical solution. It is of some interest as to whether the solution then settles on a genuine physically driven hysteretic flow. [19] The switching mechanism of Kool and Parker [1987] is a spatially local device for the (i, j) node, driven by the properties of the numerical solution, and pays no attention to what is happening to its neighboring nodes. As such, it is based only on the temporal scales inherent in the solution. A more sophisticated switching mechanism would need to consider the spatial scales and what is happening at neighboring nodes as well as look at the timescales for the change of phase. [20] This raises the question as to the appropriate size of δ for use in such a switch of phase. The progressive wave (equation (7)) and standing wave (equation (8)) solutions are obtained from a linearized form of the computational equation: the telegraph equation. These elemental solutions need to be incorporated into Fourier series, where the amplitude spectrum is determined from the initial conditions and forcing terms arising from higher-order truncation errors. These amplitudes would be extremely difficult to estimate. [22] The wavelengths of the progressive wave do suggest a length scale for the consideration of a switch of phase passing through the nodes. The standing wave suggests that the essential length scales of the finger need to be considered. The timescales of the progressive waves can be large, but the standing waves do give a crisp limit on the timescale for the switch. [23] The hyperbolic partial differential equation (3) contains a coefficient ΔτΓ*/2 on the second-order time derivative. As the time step, Δτ → 0, then equation (3) reduces to equation (1), for which the solution is sought. Rubinstein and Rubinstein [1993, theorem 21.1.1] show that the solutions of the 1-D telegraph equation do converge to the solutions of the 1-D heat equation as such a coefficient becomes small. Thus the 1-D numerical scheme with truncation error is consistent with the 1-D heat equation. The result can be carried over to two dimensions for the telegraph-heat equation pair. We have not been able to find a similar result for the nonlinear Richards equation. [24] Monotonic solutions of the Richards equation have been developed in one spatial dimension. Monotonicity is easily understood as a 1-D concept, and there does not seem to be any universally agreed upon equivalent concept in a higher number of dimensions. [25] The numerical solutions of Eliassi and Glass [2001] show numerical waves in the finger. These can arise from the discretization and the structure of the resultant truncation error. Their numerical scheme leads to the telegraph equation in linear form on the plateau of the finger. This telegraph equation admits to both standing and travelling wave solutions in the numerical error. These can be controlled by using a higher-order finite difference representation for the time derivative in the Richards equation. [26] Models of hysteresis usually incorporate a change between wetting and drying phases. These switches can be artificially triggered by the presence of numerical waves, especially where sensitive switches near machine accuracy are used. The properties of the standing and travelling waves can provide alternate switching mechanisms.
No takes yet. Share an insight, caveat, or question.
Braddock et al. (2003) studied this question.
Synapse has enriched 5 related papers. Consider them for comparative context: