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 curves

Simulazioni robuste ed efficienti del flusso dell’acqua in suolo insaturo è ancora un problema computazionale di difficile soluzione, in particolare nei casi di infiltrazione in suolo relativamente secco per la presenza di grandi gradienti nella soluzione approssimata che esaltano le nonlinearità e rendono difficile la trattazione numerica. Gli approcci tradizionali che impiegano discretizzazioni spaziali e temporali uniformi sono inefficienti e computazionalmente assai costosi, oltre che essere spesso affetti da grandi errori nella conservazione discrete della massa. Algoritmi temporalmente adattativi sono stati proposti ma sono ancora computazionalmente inefficienti in particolare in presenza di forzanti esterne ad alta variabilità, quali quelle causate dai flussi di precipitazione. Nella prima parte di questa tesi viene sperimentato l’approccio risolutivo che impiega tali algoritmi tramite il metodo delle linee (MOL) utilizzando solutori di tipo DAE (DifferentialAlgebraic Equation). Si verifica sperimentalmente che questo modo di procedure può essere preferibile a tecniche a passo fisso per alcuni casi test. Tale approccio, basato su solutori temporali di alto grado di accuratezza, quando applicabile, si dimostra efficiente e competitivo rispetto a metodi più tradizionali, mantenendo buone proprietà di robustezza e conservazione della massa, quando si utilizzino approcci mono-dimensionali. Nella seconda parte della tesi ci si focalizza più direttamente nelle proprietà di convergenza e conservazione della massa di metodi recentemente proposti in letteratura applicati a problemi eterogenei e multidimensionali. La moderna applicazione di modelli idrologici vede un arricchirsi dei processi considerati accoppiando l’equazione di Richards a modelli diversi derivanti da applicazioni metereologiche, ecologiche e bio-geo-chimiche, con la conseguente richiesta di metodi numerici più efficienti e robusti. Idealmente, si richiede a questi schemi una convergenza rapida, possibilmente quadratica, e globale in tutte le circostanze (e.g., suoli aridi, simulazioni “storm-interstorm”, grandi eterogeneità, condizioni al contorno nonlineari, eccetera). Anche il tipo di formulazione matematica dell’equazione di Richards nel continuo (ad esempio in forma completamente conservativa oppure semi-conservativa) determina diverse efficienze dei diversi metodi in funzione delle condizioni forzanti. Per questi motivi, viene implementato e testato in condizioni “estreme” un algoritmo basato sull’applicazione annidata del metodo di Newton e proposto recentemente nella letteratura. La verifica viene fatta per confronto con metodi di Newton e Newton-Picard più tradizionali su tre test sintetici ma realistici e probanti. Per ciascun test si verifica la robustezza dei diversi schemi a variazioni di griglia sia temporale che spaziale, e si confrontano i risultati in termini di velocità di convergenza, e quindi efficienza numerica. I risultati suggeriscono che il nuovo metodo “nested-Newton” risulta assai competitivo dimostrando ottime proprietà di convergenza anche nei test più difficili. Infine, l’ultima parte della tesi è volta a verificare l’efficienza computazionale nell’uso negli schemi di Picard e Newton di tecniche di “look-up” che discretizzano le curve caratteristiche delle relazioni pressione-saturazione (responsabili della nonlinearità nell’equazione di Richards) utilizzando il solutore CATHY (CATchmentHYdrology). I risultati mostrano che l’uso di tabelle di “look-up” migliora l’efficienza computazionale dei metodi alla Newton tradizionali in casi test omogenei e eterogenei rispetto ai metodi analitici qualora le curve caratteristiche possano essere accuratamente descritte da un numero di punti limitato