P. J. Ross

Numerical methods of calculating water flow in soils usually rely on solving Richards' equation, which is obtained by combining the equation of continuity with Darcy's law for water fluxes. In this study, the one‐dimensional equation of continuity was solved using tables of precomputed steady‐state fluxes read from disk. This results in simpler and faster code for applications such as field water balance models, which may use relatively large depth spacings and many repeated calculations. There is no speed penalty for the use of more complex property descriptions. Steady‐state flux tables were computed by integration of Darcy's law, which results in the same solution as Richards' equation, but in principle they could be obtained in some other way. The method is much faster and more stable than conventional iterative methods both because it is noniterative and because good estimates of steady‐state fluxes allow accurate solutions with larger vertical grid spacings, which in turn allow larger time steps. It is also both faster and simpler than other noniterative methods that provide good flux estimates because the fluxes are already available without the extra code and the time needed to compute them.

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.