draft of

Some Comparisons of New Time Steps for Mass-Conservative Infiltration Models

Abstract

Modern mass-conservative models of unsaturated infiltration of water into dry soil lend themselves to a simple feedback algorithm that adapts the time step size to maintain a given number of iterations to mass balance. Four cross-combinations of methods, fixed and adaptive, first-order and second order Runge-Kutta, are compared for equivalent time step sizes over a range of about 0.1 to 500 s. Measures of the error in total mass infiltrated, the total number of Newton iterations, and the computational efficiency (error times effort) are used to evaluate the accuracy and utility of the methods. The mass infiltration error was found to be roughly equivalent between fixed and adaptive time steps; going from first-order accurate in time to second-order times steps is more effective in reducing the overall cumulative infiltration error. But for a range of small time steps, or maximum iterations, the new adaptive second-order method has an order of magnitude better error times effort figure than the other three methods. There is a similar effect for adaptive first-order over fixed first order, but it is much less pronounced. These results demonstrate that there are still gains to be made in the accuracy and computational efficiency of mass-conservative methods for numerically simulating infiltration.

Introduction

The objectives of this work are 1) to determine if an ad-hoc but intuitively satisfying method of adaptive time stepping has any actual benefit in accuracy or computational efficiency, and 2) to reexamine claims of perfection for one modern mass conservative methods in light of further study.

It may be that anyone who has written a simulation model of the infiltration of water into dry soil has come to the conclusion that the first time steps have to be small, to maintain stability and convergence, and then have to increase to avoid very long run times. This seems to be particularly true of models that iterate to mass balance or convergence using Newton's method, where the method can become unexpectedly ill-posed, degenerating into limit cycles, as a consequence of the "stiffness" of the differential equation. Often, one puts a limit on the number of allowed iterations to convergence to put an end to such useless behavior. But if the limit is too low, then mass balance is lost, never to be regained. It may also seem that the model becomes less stable when a grid cell in a sharp wetting front approaches saturation.

This creates a strong desire to have an "adaptive" time step that finishes in a set number of iterations, or restarts if the iterations fail. We reason that if only a time step could be devised that becomes shorter or longer to finish in a "reasonable" number of iterations, the model must be more efficient. Kirkland, et al., (1992) in their analysis of the mass-conservative modified Picard method (Celia, et al., 1990), mention in their results section that a "more appropriate adaptive time step algorithm ... should be developed for the (modified Picard) algorithm".

Thirty-some years ago, in one's senior project work, there was an old American Society of Automotive Engineers report, the citation long since lost, on a electro-mechanical dithering system for optimizing ignition spark advance. The advance was moved three degrees forward and back, and if the engine RPM responded by increasing, the advance was moved three more degrees in the direction of positive increase.

Now consider a mass-conservative model in infiltration. In place of increasing engine RPM, let us consider what it takes to decrease the number of iterations necessary to achieve mass balance. Unlike the engine, we have a good idea of which direction produces a positive response. Usually, decreasing the time step size also decreases the number of necessary iterations. So instead of dithering to find the right direction, one can simply set a maximum number of iterations to be allowed. If the number is less, the size of the time step can be increased. If the number is more, it can be decreased. The resulting rules form a simple feedback system, driving the number of iterations in a time step towards a desired value.

We will use the mass-conservative modified Picard method (Celia, et al., 1990). This method uses Newton's iteration on the mixed-head-water-content form of Richards' equation [1], in finite differences, to produce the tridiagonal system of equations [2]. We will use the arithmetic interblock conductivity mean, K_i+1/2 = (K_i+K_i+1)/2, for both the gravity and capillary terms. We use it not because it is best (Baker, 1995b, 1998), but because it needs no lengthy introduction and is consistent with Celia's presentation.

[1] , where

[2]

No one should deny the great service that Celia, et al., have done to the groundwater modeling community in pointing out the simple way in which mass balance can be obtained to an arbitrary degree. But claims in the 1990 paper have been made for the benefits of the method on the apparent basis of single graphs without sufficient analysis and proof. Thus the perfection of the method has been somewhat oversold. The impact of this dominance of the field has been so great that one suspects that it accounts in some part for the current conventional wisdom in funding and consulting circles that says, "We just don't need any new models."

Therefore, it is legitimate to examine the claims in the light of further study. Two in particular stand out. On pages 1491-1492, Celia, et al., state: "The influence of maintaining mass balance is illustrated in Figure 5. ... While the errors associated with progressively larger time steps Dt in the h-based approximation led to poor mass balance and a progressively larger underprediction of infiltration depth, the mixed form maintains the correct infiltration depth but the moisture front is more diffuse as Dt increases. Imposition of global mass conservation forces the infiltration front to maintain its correct position, and the larger discretization errors appear as increased diffusion."

The second claim, on page 1492 states: "While the solutions of Figure 3 exhibited strong dependence on Dt, with deterioration of accuracy as Dt increases, the solutions of Figure 6 are insensitive to Dt. Thus in this case the mass conservative scheme serves to greatly reduce the influence of the time truncation errors, with the solutions of Figure 6 dominated by spatial error." Neither of these statements seems to take into account the possibility of simple coincidence, single graphs that simply look as though they demonstrate a benefit. As it turns out, experimental examination of a rather ad hoc adaptive time-stepping method, enhanced with a mathematically sound diagonally-implicit Runge-Kutta method, demonstrate that these claims may be somewhat overstated.

Methods

We follow here the first example of Celia, et al., [3]. But instead of the initial condition at a pressure head of -61.5 cm, which is quite moist for most field conditions (with a relative unsaturated conductivity of 3.88(10^-3)), we choose a much drier initial pressure head of -929.8 cm, corresponding to a relative unsaturated conductivity of 10^-8. This provides a much sharper wetting front, which is generally more difficult to model. Water infiltrates into the soil surface at x = xl until the end-of-run time, t = tl. We arbitrarily choose tl = 1000 s and xl = 120 cm, to provide a larger possible range of space step, dx, and time step, dt, than did Celia, et al.

[3]

We will compare four different combinations of time stepping methods. Two will use first-order accurate time steps known as the Euler-Implicit method, which can also be described as a first-order diagonally-implicit Runge-Kutta ( DIRK1) method. The other two will use a second-order diagonally-implicit Runge Kutta (DIRK2) method. One each of the DIRK1 and DIRK2 methods will use fixed time steps. The other two will use a new adaptive time step sizing method, introduced in this paper. This cross combines the fixed and adaptive time step sizes with the first- and second-order stepping methods. We will compare the methods by error in cumulative mass infiltrated, the total number of iterations per run, and an ad-hoc error-times-effort measure of computational efficiency. Since we will be using the same space step size, dx = 0.12cm, in all comparisons, the error-times-effort measure will be the total number of iterations in a run times the error in mass infiltrated. This measure equalizes runs which merely shift a reduction in error to an increase in effort, rather than offer real improvement.

We begin, just for curiosity, by setting dt = 1s and varying dx from 0.012 cm to 12 cm in step factors of 1, 2, 5 and 10. Figure 1 shows the results. Note how the mass infiltrated begins to drop quite sharply above dx = 1cm. For convenience, we will use dx = 0.12cm, corresponding to 1001 grid points, including boundary conditions, for the comparison runs.

Figure 1: Numerical solutions of [3] using [2] for dt = 1s and dx = 0.012 cm to 12 cm in steps of 1, 2, 5, 10

The first time stepping method corresponds to that used in Celia, et al., first-order Euler implicit fixed time steps of dt = 0.1s to 500s in step ratios of 1, 2, 5 and 10, with dx = 0.12cm (1000 steps). Equation [4] reprises [2], dropping the interation index, m, and replacing d with d.

[4]

We use the same convergence calculations used by Celia, et al (1990). At each Newton's iteration in a time step, the change of mass in equivalent water depth (cm), sumth [5], is calculated and compared to the net inflow, flx [6], as a fractional mass balance, rmb [7]. Iterations are continued until rmb is less than 10^-8, before the time step is allowed to advance. Reportedly, the value used by Celia, et al., was 10^-6. At the end of the run, the accumulated mass in the profile, am [8], and the total net inflow, ft [9], are used to calculate the cumulative or global mass balance fraction, cmb [10]. In some cases, the top boundary flow may be calculated separately from the bottom boundary flow, but the flow at the top boundary is always orders of magnitude larger than the flow at the bottom.

[5] and

[6]

[7]

[8]

[9]

[10]

The second time-stepping method uses a new way to vary the time step, dt or Dt, in the same first-order Euler time step. Equations [3] to [10] still apply, but the time step is varied to limit the number of iterations to a given maximum, maxit, and encouraged to increase to that maximum. Listing 1 demonstrates the general structure of this mass-conservative adaptive time stepping method, which can be adapted to other iterative methods. In the inner-most loop, the time step proceeds to convergence through iteration. If the number of iterations is maxit-2 or less, the size of the time step is increased by a factor of 1.062. If the number of iterations is maxit-1, then the size of the time step is increased by 1.001. If the number of iterations is maxit, then the time step is decreased by a factor of 0.619. If the convergence fails, or takes more steps than maxit, the results are discarded, the time step is decreased by 0.619 and the time step is restarted.

The factors of 0.619 and 1.062 are arbitrary, 0.619 happening to be about the golden section ratio. The decrease is larger than the increase because experience suggests that starting small enough wastes less iteration calculations in general. The same experience suggests that increasing time step size too fast causes unnecessarily large numbers of restarts. The factor of 1.001 keeps the time step increasing just a little for just below the chosen "optimal" number of iterations. Otherwise, it might never change.

If dt is large, but acceptable, the method will reduce the initial time step until the model converges in maxit or less iterations. If dt is small but acceptable, the method will increase it in successive time steps. In this way, throughout the run, it constantly adjusts dt so the model will reach "perfect" mass balance, or whatever other convergence criterion is used, in maxit iterations or less.

Note that this method includes a maximum time step size, dtt. If either maxit or dtt are smaller than optimum, the model will solve in many more time steps and iterations than necessary. The same can be true for maxit or dtt larger than optimum. Optimum is like gas mileage - it varies.

The third time-stepping method employs a second-order diagonally-implicit Runge-Kutta (DIRK2) method (Baker, 1993, 1995a) with a fixed time step, dt. The fourth method uses a DIRK2 method with the adaptive time step of Listing 1. This is an entire family of second-order implicit time stepping methods, of which the well-known Crank-Nicolson method is a special case. The Crank-Nicolson is usually the worst case, producing the most error and oscillation of this family, even for a simple diffusion problem.

The DIRK2 method takes two mass-conservative steps. The first [11] goes to t+p*dt, where 0 [ p [ 1/2 for stability, according to the von Neumann criterion (DuChateau and Zachmann, 1989). The second [12] goes to t+dt, where the Runge-Kutta parameters, a1 and a2 [13], determine the split between the contribution of the t+p*dt and t+dt arrays. Here again, we have deleted the iteration index, m. Each time step, partial and full, is iterated separately to mass balance convergence, according to [5] to [7].

[11]

[12]

[13]

Calculations of boundary flows for mass balance are also apportioned by the Runge-Kutta parameters. For example, the top inflow at the t+p*dt time step, flpt [14], is added to the full time step top inflow by the ratio a1 to a2 to obtain the composite top inflow, flxt [15]. The uncommented listing of the program, tbnewtc3.for, Listing 2, of the fourth and most complicated method, is offered for the reader's information. We arbitrarily pick p = 1/3, without any attempt at optimization. The reader may devolve this program back to the first three time-stepping methods to verify the results, or to play with the different possible combinations of run parameters to see what happens. Values of p closer to 1/2, for example, may work better.

[14]

[15]

Experimental Results and Discussion

Figures 2a and b show the wetting profiles in pressure head and water content for [3] using the first method, a fixed DIRK1 time step, with dx = 0.12 cm, and dt = 0.1s to 500s at steps of 1, 2, 5 and 10. The run for dt = 1000s failed. In Figure 2a, the largest values of dt produce the curves farthermost to the right. In 2b, the largest dts produce the most apparent dispersion in the wetting front. The figure looks as though it supports the contention that the modified Picard method "forces the infiltration front to maintain its correct position, and the larger discretization errors appear as increased diffusion". But we shall see differently below.

Figure 2a: Pressure Head (cm) versus Depth (cm) for DIRK1 Fixed Time Step, dx = 0.12cm, dt = 0.1s to 500s in steps of 1, 2, 5, 10	Figure 2b: Water Content versus Depth (cm) for DIRK1 Fixed Time Step

Figures 3a and b show similar wetting profiles for the second time stepping method, the adaptive DIRK1 time step, with dx = 0.12cm. Since the adaptive time steps can vary over orders of magnitude from the start to the finish of the run, the only way to make a comparison with the fixed dt is to take the average dt. This is defined as the length of the run, tl = 1000 s, divided by the number of time steps taken, nt. These plots show only a few of the runs, for the values of dt as noted in the caption. Although the plots seem to have bunched together, it is only because the largest dt is 111.11s. There is no direct way to force a larger average time step. The larger average time steps again produce the highest dispersion, and if the reader counts back along the curves, it is apparent that there is not much difference between the first two methods.

Figure 3a: Pressure Head (cm) versus Depth (cm) for DIRK1 Adaptive Time Step, dx = 0.12cm, average dt = 0.056, 5.1, 10.5, 20, 45.5 and 111 s	Figure 3b: Water Content versus Depth (cm) for DIRK1 Adaptive Time Step

Figures 4a and b show similar wetting profiles for the third method, a fixed DIRK2 time step, with dx = 0.12 cm. Here dt varies from 0.1 s to 1000 s in steps of 1, 2, 5 and 10. Note how the water content curves show that the 500 s and 1000 s time steps produce some non-physical overshoot above the inflow boundary condition. But the rest of the curves do indeed show less dispersion than either DIRK1 method.

Figure 4a: Pressure Head (cm) versus Depth (cm) for DIRK2 Fixed Time Step, dx = 0.12cm, dt = 0.1s to 1000s in steps of 1, 2, 5, 10	Figure 4b: Water Content versus Depth (cm) for DIRK2 Fixed Time Step

Figures 5a and b show some of the wetting profiles for the fourth method, a DIRK2 adaptive time step, for dx = 0.12cm and dt = 0.081, 5, 12.5, 20.4, 43.5, 111, 167 and 333 s. Only the largest time step produces noticeable dispersion from the other water content curves in 5b. But if the reader counts back from right, the results are comparable with the fixed DIRK2 time step.

Figure 5a: Pressure Head (cm) versus Depth (cm) for DIRK2 Adaptive Time Step, dx = 0.12cm, average dt = 0.081, 5, 12.5, 20.4, 111, 167 and 333 s	Figure 5b: Water Content versus Depth (cm) for DIRK2 Adaptive Time Step

It is apparent from these curves that the DIRK2 method offers some improvement over the standard Euler-implicit DIRK1 method in maintaining the shape and position of the wetting front with increasing dx. But what about the adaptive time step versus the fixed time step? This and previous experience offers notice of some general conditions and benefits.

Given a convergent model and good starting estimate of dt, this approach generally eliminates startup problems, even for infiltration into dry soil, while maintaining "perfect" mass balance. It will fail if the model is not convergent, or if the starting estimate of dt is too large or too small. If the initial dt is too large, then the summation of mass infiltrated, sumth, will be orders of magnitude smaller then the net inflow, and the model will simply bomb on an exponent error in the solution of the tridiagonal system of equations. If dt is too small, sumth, in the first time step will go to zero even when there is significant net inflow. Then, because the method is currently programmed so, it will decrease dt repeatedly until it hits the limit of 10^-100.

For more conclusive information, we look at run statistics. Figure 6 plots the total number of Newton's method iterations against time step and average time step size. Since all four methods use 1000 space steps, this provides a rough measure of model computational effort. The light line is method 1, the fixed DIRK1 time step, the plus signs mark method 2, the adaptive DIRK1 time step, the heavy line is method 3, the fixed DIRK2 time step, and the crosses mark method 4, the adaptive DIRK2 time step. When some of the outliers at the ends are discarded, where only one method is represented or is marginally stable, and straight lines are fit to the log-log data, the methods show approximate relative efforts of 1.0, 1.15, 1.41 and 1.48 times the first method, respectively.

Figure 6: Computational Effort for Four Time-Stepping Methods light line = fixed DIRK1 time step, pluses = adaptive DIRK1 time step, heavy line = fixed DIRK2 time step, crosses = adaptive DIRK2 time step

Figure 7 plots the total mass infiltrated, expressed as equivalent water depth (cm), during a run against the time step and average time steps sizes of each method. Again, the light line, pluses, heavy line and crosses mark the results for methods 1, 2, 3, and 4, respectively. We can see that both DIRK1 methods rapidly lose infiltrated mass above a time step of dt = 10 s. Both DIRK2 methods produce significantly less error, moving the sharp drop in infiltrated mass to time steps above 100 s, about an order of magnitude higher than the DIRK1 methods. At this scale, it also looks as though there is little difference between the fixed and adaptive time steps, the greatest difference being made by the order of the time step.

Figure 7: Mass Infiltration Errors for Four Time-Stepping Methods light line = fixed DIRK1 time step, pluses = adaptive DIRK1 time step, heavy line = fixed DIRK2 time step, crosses = adaptive DIRK2 time step, Correct infiltration mass is the limit as dt goes to 0, about 5.143869 cm equiv. depth

These plots, however, leave little question that the modified Picard method does not force the infiltration front to "maintain its correct position" if that position is indicated by the equivalent depth of water infiltrated. There is no question that the modified Picard method still contains time step discretization error that is reduced by using higher-order time steps.

The curves in Figure 7 look like they converge to the same amount of infiltrated mass for dt approaching 0. Let us assume that the DIRK2 methods provide the smallest error for all time steps. Runs were made for fixed dt = 0.1, 0.2, 0.3, 0.4 and 0.5 s and for the average dt in the adaptive method for 0.08098, 0.4423, 1.3514, 2.5316 and 5s. For each of these sets, two fitted polynomials were used to extrapolate the value of infiltrated mass at dt = 0. The average of these values came to am_o = am(dt=0) = 5.143869 cm. This allows the expression of absolute relative mass infiltrated error, relerr [16].

[16]

The question of whether or not the adaptive time stepping method presented here offers any numerical benefit is partially answered by Figure 8. This plots the previously mentioned measure of computation efficiency, error times effort, expressed as the total number of Newton's iterations in a run times relerr. As before, the light line, pluses, heavy line and crosses represent the results for the four respective methods. For all dt, the DIRK2 methods, even though they require significantly more effort, are computationally more efficient than the DIRK1 methods by this measure.

Figure 8: Measure of Computational Efficiency for Four Time-Stepping Methods light line = fixed DIRK1 time step, pluses = adaptive DIRK1 time step, heavy line = fixed DIRK2 time step, crosses = adaptive DIRK2 time step, measure is total number if iterations per run, tin, times the relative error, relerr

But the adaptive time steps are only better than the fixed time steps for the smaller time steps. We see that the adaptive DIRK1 time step is better than the fixed DIRK1 time step only below about dt = 2 s. The adaptive DIRK2 time step is better than the fixed DIRK2 time step only below about dt = 30 s. But as dt ranges between 0.5 and 10 s, the adaptive DIRK2 time step is about an order of magnitude or more computationally efficient than the other three methods.

Figure 9 shows the error-times-effort measure of both DIRK2 methods plotted against the iteration parameter, maxit. The corresponding values of maxit for this benefit extend from 5 to 17. One should suspect that the "noise" in this plot and the others of adaptive time step data comes from the fact that the size of the time step varies nonlinearly and in an almost pseudo-random fashion when iteration limits are exceeded. There may be a more optimal approach, based on actual error estimation, that would produce smoother results. This at least demonstrates that adaptive time stepping with higher-order time steps has the potential to offer significant improvements in computational efficiency.

Figure 9: Measure of Computational Efficiency vs Maxit for Adaptive Time-Stepping Methods pluses = adaptive DIRK1 time step, circles = adaptive DIRK2 time step, measure is total number if iterations per run, tin, times the relative error, relerr

The ability of DIRK2 time stepping methods to achieve lower infiltration error over DIRK1 methods is real and can be expected in general due to the proven ability of higher-order Runge-Kutta methods to reduce time-step discretization error. But the benefits in computational efficiency shown in this one example for the combination of a DIRK2 method with a new, relatively ad-hoc, adaptive time stepping method should be considered to be provisional, rather than conclusive. Should anyone wish to determine if this will hold for other cases, the model code for the fourth method, the adaptive DIRK2 time step, is given in Listing 2.

Conclusions

The results of these mass-conservative infiltration models demonstrate that the modified Picard method (Celia, et al., 1990) does not, as advertised, remove all time step discretization error, as reflected in total equivalent water depth infiltrated. One suspects that lack of error found by Celia, et al., may have been due to the rather moist initial conditions they chose for their first example, and possibly to a matter of coincidentally favorable graphs. In any case, Celia, et al., did not present any data on the mass infiltrated in their numerical experiments to verify the claim.

These results show that it is more effective to reduce time step discretization error by improving the order of the time step, via Runge-Kutta methods, than it is to use adaptive time steps. When the total mass infiltrated for adaptive time steps are plotted against those for fixed time steps, using the average time step for adaptive steps, the results are virtually indistinguishable. Indeed, for adaptive method presented here, using the mean time step size as comparison to fixed time steps, the computational overhead increases by about 5 to 15%.

But in this example, where the errors are already small due to time step size, the adaptive time step is more computationally efficient than the fixed time step in the computational efficiency measure of error times effort. Above a particular time step size, this effect reverses and the fixed time steps become marginally more efficient. The reversal in efficiency comes for larger time steps for 2nd-order diagonally implicit Runge-Kutta (DIRK2) methods than for the first-order (DIRK1) method, and is more pronounced. For a limited span of time step size, the adaptive DIRK2 method produced computational efficiencies an order of magnitude better than the other three combinations: fixed DIRK1, adaptive DIRK1 and fixed DIRK2. This time step span corresponds to a limit on the maximum number of iterations to mass balance of from about 5 to 17. This suggests that a limit of 15 iterations may be near optimal for this method.

No mathematical rationale is offered at this time for the apparent benefit of this ad-hoc method, which should be considered as provisional until demonstrated in other examples. To that end, the model code for mass-conservative infiltration using the adaptive DIRK2 method, which can be devolved to the other three, is offered for public use. Nevertheless, this adaptive method assures that a time step will converge to mass balance in a given number of iterations, providing at least the comfort that the model has not gone into such aberrations as endless limit cycles.

Given this evidence, one may reasonably conclude that contrary to conventional wisdom modern mass-conservative models are not "perfect". Errors such as time step discretization error are alive and well, even with mass balances to one part in 10⁸. Issues still remain to be resolved and improvements still remain to be made.

References