Implementation and testing of techniques for improving the performance of Richards' equation solvers and the handling of heterogeneous soils

Efficient, robust simulation of groundwater flow in the unsaturated zone remains computationally expensive, especially for problems characterized by sharp fronts in both space and time. Standard approaches that employ uniform spatial and temporal discretizations for the numerical solution of these problems lead to inefficient and expensive simulations. Accurate solution of the pressure-head form of Richards' equation is very difficult using standard time integration methods, because the mass balance errors grow unless very small time steps are used in the time integration process. Richards' equation may be solved for many problems more economically and robustly with variable time step size instead of constant time step size. But variable step-size methods applied to date have relied upon empirical approaches to control step size which do not explicitly control truncation error of the solution. In the first part of this thesis, we solve Richards' equation using the method of lines with a finite difference approach. We show how a differential algebraic equation implementation of the method of lines can give solutions to Richards' equation that are accurate, have good mass balance properties, and are more economical for a wide range of solution accuracy. We have implemented the method of lines using the higher order time integration ode solver ode15s to (i) assure robustness for difficult nonlinear problems and computational efficiency; (ii) develop high order adaptive methods for the time discretization taking into account the different time scales that may appear in the process; (iii) investigate the advantage of using higher-order methods in time. The numerical results demonstrate that the proposed method provides a robust and efficientalternative to standard approaches for simulating variably saturated flow in one spatial dimension. In the 2nd part of the thesis, we to investigate the convergence and mass balance behavior in Richards' equation-based solvers. As hydrological models become increasingly sophisticated (e.g., coupling with various meteorological, ecological, or biogeochemical components) and are applied in ever more computationally demanding contexts (e.g., the many realizations that are typically generated in parameter estimation, uncertainty analysis, data assimilation, or scenario studies), the need for robust, accurate, and efficient codes is greater than ever. The Richards equation for subsurface flow is highly nonlinear and requires iterative schemes for its solution. These schemes have been the subject of much research over the past two decades, but an effective all-purpose algorithm has thus far proven elusive. Ideally, rapid (quadratic as opposed to linear) and global (insensitive to initial guess) convergence is sought, in addition to applicability over a range of conditions (dry soils, storm-interstorm simulations, geological heterogeneity, 3D domains with complex boundary conditions, etc). Richards’ equation can be mathematically formulated and numerically discretized in a variety of manners, and the specific form and scheme chosen will affect the mass balance behavior of the model. We implement and test a promising nested Newton-type algorithm for solving Richards' equation. In the current state of the art, the Picard iteration method is widely used for solving the nonlinear equation governing flow in variably saturated porous media because this method is simple to code and computationally cheap. But the convergence is slow and sometimes fails. On the other hand, the Newton method is more complex and expensive than Picard. As a result the Picard method is more attractive than Newton. However, especially for strongly nonlinear flow problems, the robustness and higher rate of convergence makes Newton an attractive alternative in some cases. In this work the Picard and Newton schemes are compared with results of a nested Newton-type algorithm. Three test cases are presented and each problem is solved over a wide range of vertical discretization. The results highlight the different aspects of the performance of the iterative methods and the different factors that can affect their convergence and efficiency, including problem size, spatial and temporal discretization, convergence error norm, mass lumping, time weighting, conductivity, moisture content characteristics, and boundary conditions. It is suggested that nested Newton-type methods can be effectively implemented and used alongside Picard and Newton and numerical models of Richards' equation. In the final part of the thesis, we studied the performance of a lookup table option, as an alternative to analytical calculation, for evaluating the nonlinear soil characteristics needed in the Picard and Newton schemes. This assessment is conducted for the CATHY (CATchment Hydrology) 3D Richards-based subsurface flow solver. The lookup table method can be a cost-effective alternative to analytical evaluation in the case of heterogeneous soils, but it has not been examined in detail in the hydrological modeling literature. Two layered soil test problems are considered, and the robustness and accuracy of the lookup table method are assessed for uniform and nonuniform distributions of lookup points in the soil retention curve