Published in Agron. J. 95:1352-1361 (2003).
© American Society of Agronomy
677 S. Segoe Rd., Madison, WI 53711 USA
SYMPOSIUM PAPERS
Modeling Soil Water and Solute Transport—Fast, Simplified Numerical Solutions
P. J. Ross*
CSIRO Land and Water, Long Pocket Lab., 120 Meiers Rd., Indooroopilly, Brisbane, QLD 4068, Australia
* Corresponding author (Peter.Ross{at}csiro.au).
Received for publication August 26, 2002.
 |
ABSTRACT
|
|---|
Modeling of water and solute transport in soils is increasingly important in hydrology and agriculture, but there remain many gaps and unresolved issues. One of these, the speed and robustness of numerical solutions, is important in large-scale and stochastic modeling. A fast method was developed for solution of the Richards equation for water transport. Brooks–Corey soil hydraulic property descriptions were used with water content as the dependent variable in unsaturated regions and pressure head in saturated regions. Central time weighting was used in unsaturated conditions to improve accuracy and fully implicit weighting in saturated conditions to improve stability. Water fluxes were calculated using matric flux potentials combined with a novel spatial weighting scheme for the gravitational component. Flow across soil property interfaces was calculated by equating fluxes to and from the interface. Results on a test problem involving rainfall, surface ponding, evaporation, and drainage over 400 h showed the new method to be an order of magnitude faster and more robust than an iterative scheme typifying current practice. Execution time for the test problem with <3% error was 0.006 s on a 166-MHz personal computer. Solutions for solute transport described by the advection–dispersion equation were obtained concurrently without substantial increase in execution time by assuming average water fluxes over several steps of the water transport solution. The methods presented here should also be applicable to two- and three-dimensional problems.
|
INTRODUCTION
|
|---|
MUCH EFFORT HAS BEEN EXPENDED in recent years on modeling soil water and solute transport, and much progress has been made, but many gaps and unresolved issues remain. Many of these are due more to difficulties in obtaining information on the composition and structure of the soil itself than to limitations of transport modeling. There seems little point in assuming a knowledge of soil properties that cannot be obtained, either because of experimental limitations or expenditure of labor, except perhaps for exploratory studies. Particularly for large-scale modeling, for which there is increasing pressure due to requirements of environmental management, obtaining even minimal soil transport properties is often a problem. Then there is the well-known but difficult area of nonequilibrium flow due, for example, to soil structural features (cracks, channels, etc.), water repellence, fingering of flow, or (for solute particularly) normal spatial variation of soil properties (Bouma, 1981; Glass et al., 1988; Kluitenberg and Horton, 1990; Ogden et al., 1992). For solute, there are also nonequilibrium chemical exchanges with the soil. Nor can transport properties always be taken as permanent soil characteristics because they can change, particularly near the soil surface, in response to biology, climate, and management (Gupta et al., 1991). Uptake of water and solutes by plants is another difficult area, and ongoing research will change the modeling of these processes.
Not all progress, however, depends on better soil characterization and research into soil and plant processes. There are many reports of modeling results in the literature, but it is often difficult to assess their general implications because of differing system parameters and input variables and because a thorough analysis of sensitivity to uncertainties in parameters and inputs is often not given. Of course, such uncertainty analysis is not easy, especially when it concerns patterns such as climatic variation (Janssen et al., 1994). In addition, many processes such as those mentioned above are ignored through lack of information, and resulting uncertainties are difficult or impossible to quantify. Some framework for assessing modeling results would be helpful. This could include a classification of climatic scenarios as well as better uncertainty analysis of parameters and processes in the model for the different scenarios.
Although past development of efficient numerical methods combined with ever faster computers has removed concerns about computing time for many purposes, and efficient software is readily available (e.g., Simunek et al., 1998; Verburg et al., 1996), there are still important areas where such concerns may prevent the use of numerical solutions of the Richards equation for water transport and the advection–dispersion equation for solute transport. One of these is in large-scale modeling, mentioned above, and another is in stochastic modeling, which is helpful in uncertainty assessment (e.g., Bresler, 1987). Computational time can also become a concern in two- and particularly three-dimensional modeling. This paper presents methods to greatly increase the speed of solution of the equations and thereby increase their range of application. Although the soil hydraulic properties chosen are simple in form, in many applications, detailed knowledge is unavailable or not particularly relevant. Uptake of water and solute by plants has not been included because it is not an essential aspect of the methods, but for many applications, it would be necessary to include it and perhaps to combine modeling of water and solute transport with modeling of plant growth.
|
MATERIALS AND METHODS
|
|---|
Soil hydraulic properties are described by the Brooks and Corey (1964) relations:
 | [1] |
where S is degree of effective saturation; 
is volumetric water content (m3/m3); K is hydraulic conductivity (m/s); S and KS are water content and conductivity at natural saturation, respectively; r is residual water content; h is pressure head; and he is pressure head at air entry (negative). The parameters S, r, KS, and he are scaling parameters while dimensionless parameters 
and 
define the shapes of the water retention and conductivity curves and are commonly taken to be related, e.g., = 2/ + 3 (Brooks and Corey, 1964; Campbell, 1974). The parameter r depends on soil hysteresis (Haverkamp et al., 2002), which is not considered here, and will be taken as zero. The matric flux potential 
(Gardner, 1958; Campbell, 1985) is given by
 | [2] |
where e = KShe/(1 - ). It is also known as the Kirchhoff transform and is frequently used in both analytical and numerical work because it linearizes the matric potential term in Darcy's law for vertical downward flow:
 | [3] |
where q is water flux (m/s) and z is depth, allowing simple calculations of flux to be accurate over larger distances than if h is used. Note that for saturated soil (h 
he), is linearly related to h and gives the same fluxes as using h. Applying mass balance to Darcy's law yields the Richards equation. Similarly, applying mass balance to the expression for solute flux yields the advection–dispersion equation.
To numerically solve the Richards equation and the advection–dispersion equation, both space and time must be divided into discrete intervals 
x and t. The differential equations can then be converted to a finite set of algebraic equations whose solution provides an approximation to that of the differential equations. The most common method of conversion for the Richards equation leads to a set of nonlinear algebraic equations that must be solved iteratively (e.g., Ross and Bristow, 1990). Iteration is also needed for the advection–dispersion equation if nonlinear chemical exchange isotherms are used. The methods followed here yield sets of linear equations that can be solved without iteration.
Discretisation
The soil is divided conceptually into n horizontal layers of thicknesses xi, i = 1,2,...,n (Fig. 1)
. Water volume or solute mass per unit area for layer i is denoted by Qi, and flux from layer i to layer i + 1 is denoted by qi. The mass balance equation for layer i is
 | [4] |
View larger version (14K):
[in this window]
[in a new window]
|
Fig. 1. Discretisation of the soil profile into n layers of thickness xi with centers at zi, i = 1,2,...,n with water or solute fluxes qi between them.
|
|
For water, Qi = xii = xi[ri + (Si - ri)Si], where i and Si are the mean values of water content and degree of saturation for layer i. For solute, Qi = xisi, where si (kg/m3) is the mean mass of solute per unit volume of soil for layer i. So far, no error has been incurred by the spatial discretisation. However, fluxes calculated from a knowledge of i, Si, and si are approximate and introduce one of the two main sources of error. The other source arises from dividing time into steps t. The mass balance equation can then be written
 | [5] |
for the time step from t to t + t where Qi is the change in Qi and the fluxes are to be evaluated a fraction 
through the step. This equation is also exact if exact fluxes and the exact value of required are calculable, but alas they are not. However, mass balance is accurately maintained in the sense that Eq. [5] is satisfied whether or not the fluxes are a good approximation. For water, Qi = xii = xi(Si - ri)Si while for solute, Qi = xisi.
Water
The water fluxes at a fraction through the time step are evaluated by linearizing qi in terms of Si and Si+1, giving
 | [6] |
where the superscript 0 indicates the beginning of the time step. The flux q0 at the soil surface and the drainage flux qn at the bottom of the profile must also be known or calculable. Substituting Eq. [6] into Eq. [5], we obtain a tridiagonal set of n linear equations
 | [7] |
where a1 and cn are zero because S0 and Sn+1 do not exist. Such a set of equations can be solved very efficiently using the Thomas algorithm (see Campbell, 1985). For the results reported in this paper, the value of was taken as 1 when any part of the profile was saturated (to improve stability) and 0.5 otherwise (to improve accuracy). Note that solving for the Si ensures that the mass balance equation is satisfied; if the i were used instead, the equation would no longer be exactly satisfied in unsaturated regions. Note also that no iteration is needed to obtain the Si. In the approach usually employed, the fluxes are nonlinearly related to the i or hi, and a solution of the nonlinear equations is found by iteration based on the linearized equations. Because the iterations can fail to converge, this procedure is less robust than that described here, and as will be seen later, there is not enough gain in accuracy to offset the extra work.
Special Cases
The above equations apply only where the soil is unsaturated because Si is zero when layer i is saturated (hi hei). In saturated layers, h or, equivalently, the matric flux potential is used in place of S. Because conductivity and water content are constant in saturated regions, the fluxes given by the linear equations will be exact there.
If ponding occurs on the soil surface, another equation for the pond height h0 must be included:
 | [8] |
where qprec and qevap are the rates of precipitation and potential evaporation, respectively.
Choice of Time Step
The larger the step size t is, the faster the solution and the greater the error. Specifying a maximum allowable change Smax for the degree of saturation in any layer allows control of the step size and the error. A relative change can also be used, but the absolute change seems to perform a little better. A safety limit of some maximum fraction by which S can be decreased (e.g., 0.9) could conceivably be needed sometimes. Calculation of the step size involves an initial estimation
 | [9] |
using the rate of change of greatest magnitude at the beginning of the time step. This is followed by checking of the actual Si after solution of the equations. If any |Si| is greater than (1 + E1)Smax, where E1 is some allowable excess (e.g., 0.25), then the time step is recalculated as
 | [10] |
and the equations solved again after recalculating the bi coefficients. Note that the fluxes do not need to be recalculated. This procedure takes only a relatively small part of the calculation time for the step and so is an efficient way to adjust the step size on the few occasions when this is necessary.
The initial estimate of t may need to be reduced to limit the step size to some desired maximum or to avoid overshooting the finish time. If the profile is fully saturated, the initial estimate cannot be calculated because the maximum dS/dt is zero. An estimate can be obtained by comparing the input rate qprec - qevap with the drainage rate qn. If the input rate is greater, then the profile will remain saturated, and time can be advanced to the finish time, or the maximum step size can be used (but note that the equations must include a surface pond for there to be a unique solution). If the drainage rate is greater, then the profile will desaturate when the surface pond has infiltrated, and the time increment for this to occur can be estimated as
 | [11] |
where h0min is some slight negative value for acceptable overshoot (e.g., -0.02). Even when the profile is not fully saturated, it is worthwhile estimating the time when the pond will disappear and adjusting t if necessary to
 | [12] |
Overshoot of h0min can be checked after solution of the equations, and if the change h0 is too large, the step size is reduced to
 | [13] |
Usually h0min can be taken small enough to be negligible, but if a lower limit of zero is required for cosmetic reasons, the required amount of water can be extracted from the surface layer after the values of S have been updated.
Another adjustment of step size may be needed to avoid unacceptable overshoot of S above unity when a layer saturates. If a maximum acceptable value Smax (e.g., 1.001) is exceeded, the step size should be reduced to
 | [14] |
Again, water could be moved around after updating to reduce S to 1 for cosmetic reasons; this was not done here. Adjustment of the time step for overshoot of when a layer desaturates does not seem to be necessary or desirable; it is enough to designate the layer unsaturated with S (1 
S < Smax) left alone.
A further adjustment is associated with the onset of ponding. Although the step size is reduced to avoid overshoot of S when the surface layer saturates, ponding does not occur until the pressure head rises from he1 to zero. To avoid errors in the pond height (or in runoff if ponding is not allowed), it is advisable to check for this rise by monitoring 1 and to reduce the step size to a small value when it occurs because the matric flux potential values in saturated regions can change very quickly and the current ones are only estimates made from information at the beginning of the previous time step.
Calculation of Fluxes
The fluxes qi are normally equated to estimates of steady-state fluxes from layer to layer, i.e., changes in water content are ignored. This procedure incurs errors when changes in water content are significant and particularly when layer thicknesses become large (compared with he). In wet soils, the estimates can be based on pressure head differences. For uniform soil, the flux q from Point 1 to Point 2 is often taken as
 | [15] |
where z is the distance between the points and g is the cosine of the angle between the direction from 1 to 2 and the downward vertical (for the work reported here, g = 1, but for completeness, it will be included in the equations). In drier soils, these estimates are poor because conductivity varies greatly spatially and it is difficult to estimate an appropriate average value. The matric flux potential incorporates conductivity and provides a better option:
 | [16] |
where the weight w is often taken as 0.5 as for the pressure head estimate. Note that in saturated soil, this equation is identical to the previous one, whatever the value of w. For horizontal flow (g = 0), only the first term occurs, the other two terms estimate the flux due to gravity, and the relation gives the exact steady-state flux. The matric flux potential is not continuous across boundaries between different soil types, and special treatment is needed for gradational soils, but its use can be effective nonetheless (Ross and Bristow, 1990) as long as hysteresis is not present in the relation between conductivity and pressure head. A previously unresolved difficulty can occur during flow from a drier to a wetter region when layers are thick (z is of the same order of magnitude as he) and a fixed weight w is used. The flux due to gravity can then be greatly in error, and the net flux may even be in the wrong direction. This can lead to instability in the solution. Figure 2 shows an example of the flux q1 from Layer 1 to Layer 2 (g = 1) plotted against an increasing pressure head in Layer 2. The pressure head in Layer 1 is 3he and z = -he. The exact flux is zero when h2 = 2he (zero hydraulic head difference) and then reverses direction, but the estimate never reverses and gets worse as Layer 2 approaches saturation (h2 approaches he).
View larger version (13K):
[in this window]
[in a new window]
|
Fig. 2. Exact water flux from Layer 1 to Layer 2 below compared with that calculated from Eq. [16] with a fixed weighting of w = 0.5 for the gravity component.
|
|
This problem can be resolved by calculating the weight w so that the flux is zero when the pressure head at Point 1 (assumed to be higher than Point 2) is h2 - gz. Then
 | [17] |
where and K are expressed as functions of h2 - z. To avoid extra power function evaluations, which take time, and to avoid trouble with round-off errors for large negative h2 (32-bit arithmetic is adequate for the work reported here), a simple approximation can usually be used when Layer 2 is unsaturated (h2 < he):
 | [18] |
where is the shape parameter in the K(h) relation. This is exact for = 2 and = 3, with a maximum absolute error of 0.007 in between. For >3, the absolute error is below 0.01 provided v < 4/( - 3). Because <5 for most field soils, and normally v < 2, the approximation can be used most of the time. Note that for small v, w 
0.5. When Layer 2 is saturated, Eq. [17] can be used. The weight w approaches zero as h2 - gz approaches he and so is set to zero when h2 - gz he. Figure 3 shows the great improvement obtained with this weighting scheme (WS1). However, the flux calculated when either layer is saturated is less accurate than in unsaturated conditions. A much more accurate value (WS2 in Fig. 3) can be obtained by treating saturated and unsaturated flow separately. Assuming initially that the saturated front is at a known position to calculate w, and that w remains constant, equating fluxes in saturated and unsaturated regions leads to a quadratic equation for a new estimate of the front position (see Appendix). Then a new estimate of w can be calculated and the process repeated until the front position changes by less than a desired amount. Convergence is so fast that it is seldom necessary to calculate the position more than twice.
View larger version (15K):
[in this window]
[in a new window]
|
Fig. 3. Water fluxes as in Fig. 2 compared with: WS1, flux calculated from Eq. [16] with weighting w for the gravity component calculated from Eq. [17]; WS2, flux calculated by equating fluxes through unsaturated and saturated regions.
|
|
The calculation of fluxes at the soil surface in the presence of ponding was described above. In the absence of ponding, a maximum possible flux toward the soil surface is calculated by assuming that = 0 there and calculating the flux from the center of Layer 1 over the distance x1/2 in the usual way. If this flux is greater than the potential evaporative flux qevap less the precipitation rate qprec, then the flux q0 is set to -(qevap - qprec); otherwise, it is set to the negative of the maximum flux. The flux out of the base of the soil profile is taken equal to the conductivity of the bottom layer on the assumption that there is a unit hydraulic gradient there. A fixed potential can also be assumed and the flux calculated in a similar way to the maximum evaporative flux.
When layers i and i + 1 have different soil hydraulic properties, the above method for calculating qi must be modified. Because in steady-state flow the flux qi1 toward the interface in layer i is equal to the flux qi2 from the interface in layer i + 1, these two fluxes can be calculated in terms of the unknown interface pressure head and equated to obtain the unknown head. In practice, some modification of this simple procedure is desirable to avoid having to solve the resulting nonlinear equation from scratch every time step. Because linearization is to be used, it is preferable to work with the matric flux potential rather than with the pressure head h. Although is discontinuous across the interface, that on the upper side (i1) can be used to calculate h, and because h is continuous across the interface, this enables on the lower side (i2) to be calculated. Because i1 is known initially (or can be found or approximated), and can be updated at each time step, it is always available for calculating the flux qi1 toward the interface in layer i; and because i2 can be found from it, the flux qi2 away from the interface in layer i + 1 can likewise be found. To update i1, Newton iteration is used:
 | [19] |
where i1c is the current value of i1. This is continued until the change i1 is small enough when a final calculation of qi1 is made. The method is effective because a single update without iteration usually suffices. A rather loose criterion (e.g., |i1| less than half the mean of the old and new values of i1) seems to be adequate. To keep i1 positive in case of a poor initial estimate, i1 can be limited to >-0.5i1.
Expressions for the derivative f'(i1) and the partial derivatives 
qi/Si, etc., required for the tridiagonal set of equations are given in the Appendix.
Solute
Solute fluxes are affected by changes in soil water content as well as changes in solute concentration s; therefore, for solute we have
 | [20] |
The changes are known from solution of the Richards equation, and so we have a tridiagonal set of linear equations similar to Eq. [7] for water but with the coefficient di including the terms in . The fluxes are given by
 | [21] |
where 
is dispersivity (m), qw is average soil water flux, and cw is the solute concentration in the soil water (kg/m3). The water flux is known from average boundary fluxes and changes (assumed uniform) in water content over the time period in question. The contribution of diffusion to solute dispersion has been ignored. Layer thicknesses should not be too much greater than dispersivity, or predicted concentrations may oscillate; z < 2 is desirable, but if the dispersivity is of the order of one-tenth the size of the flow domain, this is not a significant constraint. The solute concentration s is related to the concentration cw in the soil water through the (possibly nonlinear) exchange isotherm fx(cw) (kg/kg) by s = cw + 
bfx(cw), where b is the soil bulk density (kg/m3). Changes in cw are therefore related to changes in s and by
 | [22] |
allowing evaluation of the derivatives in Eq. [20] as given in the Appendix. Some nonlinear isotherms may require iteration to update cw after s has been updated, but Newton's method as in Eq. [19] is normally very efficient for this. For boundary conditions, we assume that solute concentrations of any input water fluxes are known while concentrations in output fluxes (except evaporation) are taken to be equal to those in the soil water. Any number of solutes may be included, as long as they behave independently. Time stepsize can be estimated similarly to that for water (Eq. [9]), or if the solution period is tied to the time steps for water, a preset number of steps (perhaps one) can be used.
|
RESULTS AND DISCUSSION
|
|---|
Computer code for calculating water and solute transport using the numerical methods presented here was written and verified by comparison with the SWIM program (Verburg et al., 1996). The solution of the Richards equation by the new method (referred to as the fast method) was then compared with the solution by the iterative method of Ross and Bristow (1990) (referred to as the standard method). Note that, unlike the fast method, the standard method is not limited to Brooks–Corey soil hydraulic properties. Both methods conserved mass, the fast method by updating water content and the standard method by ensuring (through iteration) that updated pressure heads were in accord with the soil water retention curve. The standard method positioned nodes at layer boundaries, including the soil surface and the bottom of the profile, rather than at the center of the layers. It calculated fluxes by the -based Eq. [16] using w = 0.5, and the time step size was halved if solution of the nonlinear equations could not be obtained in 20 iterations. The time discretisation was fully implicit (equivalent to choosing = 1). For comparison with the fast method, it was modified to use time steps chosen with Eq. [9]. A minor modification to limit the change in h at the surface on desaturation was found necessary to avoid failure of the standard method. The fast method was used with parameters such as Smax set to the example values given in the previous section. Fluxes were calculated as discussed previously, using both weighting schemes WS1 and WS2 (Fig. 3) for comparison.
Most current numerical solutions use similar iterative methods although for reasons of flexibility, is seldom used. For example, the HYDRUS code (Simunek et al., 1998) uses h to calculate fluxes according to Eq. [15] while the SWIM code (Verburg et al., 1996) uses an inverse hyperbolic sine transform. Because the aim here is to examine improvements over existing methods and because is more efficient than h or other transforms, it seems reasonable to regard the method of Ross and Bristow (1990) as the standard for comparison. However, a solution using fluxes calculated with the h-based Eq. [15] was also obtained. This is referred to here as the h fluxes method. (Note that this is quite different from solving the h-based Richards equation, where the /t term is replaced by C(h)h/t, with C(h) as the water capacity. Numerical solution of this equation can produce substantial mass balance errors unless time steps are small. All the methods discussed here maintain accurate mass balance, as discussed in connection with Eq. [5].)
Comparison of Results of Standard and Fast Methods
A test problem was chosen with rainfall at 25 mm h-1 for 10 h and potential evaporation of 0.5 mm h-1 for 400 h applied to a soil profile consisting of 40 cm of sandy soil overlying 160 cm of clayey soil (Fig. 4)
at an initial matric head of -1000 cm. The sandy soil, with he = -10 cm, was divided into four 10-cm layers and the clayey soil, with he = -40 cm, into four 20-cm layers immediately below the interface and a further two 40-cm layers at depth, giving 10 layers altogether. Solution with layers twice as thick was possible for the fast method (with twice the error), but the standard method gave ridiculous results due to errors in calculating the evaporative flux. Sensible results, but with large errors, could be obtained by implementing the improved weighting scheme. Solution by the standard method using 256 layers was used as a reference (referred to in the figures as the accurate solution). In all cases, time steps were made small so that they did not significantly affect the results, which therefore reflect mostly discretisation errors (very small for the accurate solution).
Results for the depth of ponding (which is potentially runoff) over the first 20 h are shown in Fig. 5
. Both standard and fast methods using the matric flux potential gave results very close to the standard while the h fluxes solution was inferior. Evaporation and drainage over the 400 h are shown in Fig. 6
. Again the fast method was excellent while the standard -based method overestimated evaporation but not as much as the h fluxes method, which had unacceptably large errors. The overestimate for the -based standard method may have been due to evaporation from storage associated with the node at the surface; with a 10-cm-thick surface layer, evaporation occurs freely from the surface 5 cm of soil. If so, improvement could be effected by modifying the calculation of evaporation to avoid this. It may then be possible to obtain better results from layers twice as thick using the improved weighting scheme for flux calculation. All methods estimated drainage well, but for the h fluxes method, this could have been fortuitous because the large overestimate of evaporation resulted in less water in the profile before drainage started.
There was no significant difference between weighting schemes WS1 and WS2 for this test problem. Further tests using a range of soil hydraulic properties reinforced this conclusion, the greatest difference found being 1.9 mm in 56 mm of evaporation. Presumably when flow occurred from unsaturated to saturated regions, or vice-versa, fluxes were being controlled by conditions elsewhere in the profile.
Errors Associated with Time Step Size
A major advantage of the new fast method is its speed, which arises partly as a result of not needing the several iterations required to solve a set of nonlinear equations at each time step and partly from better approximations such as the use of = 0.5 in unsaturated soil. As the time step size was increased using larger values of the control parameter Smax, the time discretisation errors increased. Values of Smax up to 0.4 were used for the fast method, but for the standard method, failure of iterations to converge and consequent halving of the step size eventually increased the number of steps as Smax was made larger. The drainage error was the greatest and so was used for comparison of the two methods. Figure 7
shows errors relative to results for large numbers of very small time steps plotted against the total number of evaluations of soil hydraulic properties and fluxes, which largely determines the execution time. For the fast method, this is just the number of time steps while for the standard method, it is the product of the number of steps and the average number of iterations per step. The filled symbols are for 10 layers; for comparison, a few results for 20 layers are included (open symbols). Because these lie on almost the same lines, the results appear to be general. The slopes of the lines indicate an error proportional to t for the standard method, as expected from use of = 1, and an error almost proportional to t2 for the fast method where = 0.5 is used for much of the time. In addition, there is the gain from not iterating; it seems that solution of the nonlinear equations does little to provide a more accurate solution. In the error range 0.01 to 0.1, the improvement ranges from a factor of 45 to a factor of 20. It seems likely that modifying the standard method to use = 0.5 for unsaturated conditions would improve it considerably.
View larger version (32K):
[in this window]
[in a new window]
|
Fig. 7. Relative errors in drainage for standard (triangles) and fast (squares) solutions of the test problem as a function of the total number of iterations through the flux calculations. Filled symbols are for 10 soil layers and unfilled are for 20.
|
|
Applications
The fast method should prove useful in applications requiring many solutions of the Richards equation, such as modeling the surface hydrology of large areas or obtaining statistical parameters for populations of soil profiles. Although it is as yet untested in such applications, expected execution times can be estimated. Execution time per time step per soil layer was approximately 6 µs on a 166-MHz Pentium processor. Assuming 100 time steps with 10 soil layers to give a relative drainage error of 0.02 for the test problem, 167 such problems could be solved per second, or approximately half a million per hour. Because 166 MHz is slow by current standards, this is a considerable underestimate. The fast method appeared to be very robust, with no failures, but it has not yet been thoroughly tested. The standard method sometimes did not achieve convergence at larger step sizes and needed to halve the step size, particularly when boundary conditions changed. Failures in applications involving millions of solutions could cause problems.
Solute Transport
The methods described in this paper are applicable to solute transport described by the advection–dispersion equation. Because iteration is used to solve this equation only when nonlinear exchange isotherms are present, substantial improvement in speed due to the fast methodology could only be expected then. However, because this equation is usually solved in conjunction with the Richards equation, reduction in the number of time steps used to solve the latter could speed up solution of the former provided that larger time steps do not introduce substantial errors. This raises the question of further decoupling the solutions of the two equations by averaging water fluxes over longer times and controlling errors associated with time step size independently or by relating the time step size for solute transport to that for water transport, for example by using a fixed number of steps ns per nw steps of the water solution. To investigate this further, the advection–dispersion equation was solved to estimate solute transport during the test problem with 100 time steps for the water solution and many small time steps for the solute solution (ns large to minimize error). A noninteracting solute was initially concentrated in the 0- to 10-cm layer; dispersivity was assumed to be 20 cm. Figure 8
shows that even with only four periods (nw = 25), results differ little from those where solute was treated every time step (nw = 1). Even averaging over the entire 400 h gave sensible, though not very accurate, results. Figure 9
shows results for varying (ns, nw) pairs. There is no significant error in using one solute step every five water steps and very little in 1 every 10. Two solute steps every 20 water steps is slightly less accurate. Interacting solutes move more slowly than the noninteracting solute used here, so even greater gains are possible for them. It appears that speed can be gained by decoupling solutions of the water and solute transport equations over substantial time periods. Because the execution time per step per soil layer was only 2 µs, one-third of the time for the water solution, the gain in speed is not great with only one solute; but for several, when some of the interacting ones have nonlinear isotherms, the gain could be substantial.
View larger version (18K):
[in this window]
[in a new window]
|
Fig. 8. Solutions for solute concentration (arbitrary mass units per cm3 of soil) at 400 h for the test problem with varying numbers of time periods for averaging the water fluxes. There were initially 1000 mass units of solute in the top 10-cm layer.
|
|
View larger version (15K):
[in this window]
[in a new window]
|
Fig. 9. Solutions for solute concentration (arbitrary mass units per cm3 of soil) at 400 h for the test problem with varying pairs (ns, nw) of ns solute time steps per nw water time steps with water fluxes averaged over the nw water steps. There were initially 1000 mass units of solute in the top 10-cm layer.
|
|
|
CONCLUSIONS
|
|---|
Many difficult problems in modeling soil water and solute transport remain. To improve speed in large-scale and stochastic modeling, a new fast method for solving the Richards equation for water transport and the advection–dispersion equation for solute transport was developed and tested. For the Richards equation, speed increases of 20 to 45 times were achieved for a test problem involving infiltration, surface ponding, drainage, and evaporation compared with an iterative method from the literature, taken to typify current practice. Simple Brooks–Corey type soil hydraulic property descriptions enabled saturated conditions to be treated efficiently while using degree of saturation as the main variable in unsaturated conditions enabled mass balance to be achieved in the absence of iteration. A new spatial weighting for the gravity component of the water flux increased robustness as did the fast approach. A time-centered weighting under unsaturated conditions ensured accuracy with large time steps while a fully implicit weighting when saturated regions were present ensured stability. Calculation of evaporative flux was superior to that of the iterative method. It seems likely that iterative methods could be improved substantially by employing the same spatial and temporal weighting and by employing better calculations of evaporation. Water fluxes calculated using matric flux potentials gave better results for the test problem than fluxes calculated using pressure heads. Solutions for solute transport were obtained very efficiently by decoupling the solution from that of water transport over substantial time periods, assuming average water fluxes over several steps of the water transport solution. There appears to be no reason why the methods presented here could not be extended to solution of the two- and three-dimensional transport equations.
|
APPENDIX
|
|---|
Mixed Saturated/Unsaturated Flow
If flow is occurring between Points 1 and 2, with 1 above 2, and zS0 is the initial estimate of the length of the saturated region between the points (e.g., zS0 = z/2), then by equating fluxes in saturated and unsaturated regions, we obtain a quadratic equation a(zS)2 + bzS + c = 0, where for Point 1 unsaturated
 | [23] |
and for Point 2 unsaturated
 | [24] |
In the equation, w should be calculated from Eq. [17] or Eq. [18] using z - zS in place of z and he in place of h2 if Point 1 is unsaturated. Because zS is unknown, zS0 is used instead. A better estimate is then obtained by solving the equation and the process repeated until successive estimates differ by less than a desired amount (e.g., 0.01z). The flux is then given by (e - 2)/zS + KS for Point 1 unsaturated and by (1 - e)/zS + KS for Point 2 unsaturated. Differentiating the quadratic equation with respect to u, where u can be 1 or 2 or S1 or S2, gives zS/u, which allows derivatives of q to be found:
 | [25] |
where S = e - 2 for Point 1 unsaturated and 1 - e for Point 2 unsaturated. It has been assumed for calculating derivatives that w is constant.
Expressions for Other Derivatives
Partial derivatives are obtained by differentiating the equations for the fluxes. We have
 | [26] |
Strictly, the variation of w with Si+1 should be included, but except in wet conditions, it is negligible, and its inclusion seems to effect little improvement. Expressions for the derivatives d/dS, etc., for the Brooks–Corey soil hydraulic properties are
 | [27] |
In saturated regions, is used in place of S, and we have
 | [28] |
Flow across a soil interface is a little more complicated. Suppose the matric flux potential is i1 just above the interface in layer i and i2 below in layer i+1, and let the fluxes to and from the interface in the downward direction be qi1 and qi2, respectively. Then, given i1 (and therefore i2), qi1 and its partial derivatives with respect to Si (or i, if layer i is saturated) and i1 can be calculated in the usual way, using a value xi/2 for z. Similarly, qi2 and its partial derivatives with respect to i2 and Si+1 (or i+1 if layer i + 1 is saturated) can be calculated using a value xi+1/2 for z. In a computer program, the same subroutine is used as for other fluxes and derivatives. To find qi, the flux from layer i to layer i + 1, and its derivatives, it is necessary to allow for the changes in i1 and i2 that occur to maintain qi1 = qi2. We have
 | [29] |
where Si is regarded as the independent variable, Si+1 as fixed, and the flux as a function of Si and i1. The partial derivatives are those mentioned above while using the interface pressure head hi1 = hi2, we have
 | [30] |
where Ki1 and Ki2 are the conductivities at the interface in layers i and i + 1, respectively (these are also needed for the calculation of the fluxes qi1 and qi2). The derivative dqi/dSi is the partial derivative qi/Si required, i.e., the rate of change of qi with respect to Si with Si+1 held constant. Similarly, we find qi/Si+1 by holding Si constant and evaluating
 | [31] |
Equations for rather than S are the same except that S is replaced by . Equations for the Newton iteration to adjust i1 are simpler. We have
 | [32] |
where the derivatives are as above.
For solute, from Eq. [21] and [22],
 | [33] |
with the term di in the tridiagonal set of equations given by
 | [34] |
|
REFERENCES
|
|---|
- Bouma, J. 1981. Soil morphology and preferential flow along macropores. Agric. Water Manage. 3:235–250.
- Bresler, E. 1987. Modeling of flow, transport, and crop yield in spatially variable fields. Adv. Soil Sci. 7:1–51.
- Brooks, R.H., and A.T. Corey. 1964. Hydraulic properties of porous media. Hydrology Paper 3. Civil Eng. Dep., Colorado State Univ., Fort Collins.
- Campbell, G.S. 1974. A simple method for determining unsaturated conductivity from moisture retention data. Soil Sci. 117:311–314.
- Campbell, G.S. 1985. Soil physics with BASIC. Elsevier, New York.
- Gardner, W.R. 1958. Some steady-state solutions to the unsaturated flow equation with application to evaporation from a water table. Soil Sci. 85:228–232.
- Glass, R.J., T.S. Steenhuis, and J.-Y. Parlange. 1988. Wetting front instability as a rapid and far-reaching hydrologic process in the vadose zone. J. Contam. Hydrol. 3:207–226.
- Gupta, S.C., B. Lowery, J.F. Moncrief, and W.E. Larson. 1991. Modeling tillage effects on soil physical properties. Soil Tillage Res. 20:293–318.
- Haverkamp, R., P. Reggiani, P.J. Ross, and J.-Y. Parlange. 2002. Soil water hysteresis prediction model based on theory and geometric scaling. Geophys. Monogr. 129:213–246.
- Janssen, P.H.M., P.S.C. Heuberger, and R. Sanders. 1994. UNCSAM: A tool for automating sensitivity and uncertainty analysis. Environ. Software 9:1–11.
- Kluitenberg, G.J., and R. Horton. 1990. Effect of solute application method on preferential transport of solutes in soil. Geoderma 46:283–297.
- Ogden, C.B., R.J. Wagenet, H.M. van Es, and J.L. Hutson. 1992. Quantification and modeling of macropore drainage. Geoderma 55:17–35.
- Ross, P.J., and K.L. Bristow. 1990. Simulating water movement in layered and gradational soils using the Kirchhoff transform. Soil Sci. Soc. Am. J. 54:1519–1524.[Abstract/Free Full Text]
- Simunek, J., M. Sejna, and M.Th. van Genuchten. 1998. The HYDRUS-1D software package for simulating the one-dimensional movement of water, heat, and multiple solutes in variably-saturated media. U.S. Salinity Lab., Riverside, CA.
- Verburg, K., P.J. Ross, and K.L. Bristow. 1996. SWIMv2.1 user manual. Div. Rep. 130. CSIRO Div. of Soils, Australia.
TOP
ABSTRACT
INTRODUCTION
