Simulating Unsaturated Flow With a Finite Volume Method

Dragan Vidović1, Milan Dotlić1, Boris Pokorni1, Milenko Pušić2 and Milan Dimkić1




1 Jaroslav Černi Institute for the Development of Water Resources, Jaroslava Černog 80, 11226, Belgrade, Serbia

2 University of Belgrade, Faculty of Mining and Geology, 11000 Belgrade, Serbia





Method used to average the hydraulic conductivity is known to have a decisive influence on the accuracy of numerical schemes for Richards' equation, particularly if coarse grids are used. All standard averages lead to very inaccurate solutions in some cases. Averaging techniques that address this issue available in the literature were mostly formulated only in the one-dimensional case. We propose a framework to deploy a one-dimensional averaging method in a three-dimensional unstructured finite volume code. The method is demonstrated to give accurate results in all flow regimes.

Keywords: Finite volume method, porous medium, Richards’ equation





Water shares the Earth's crust pore space with other fluid phases, such as air or oil. Influence of these phases on water flow is often decisive so they must be taken into account in groundwater flow simulations.

Physical model of multiphase flow consists of a single conservation equation for each phase and correlations for the capillary pressure necessary to close the system.

Assumption that the air pressure is constant reduces the two-phase system to a single equation of water mass conservation. This is known as Richards' equation. It is essentially the same groundwater equation used in the saturated case but with the hydraulic conductivity and the storage depending on the pressure in a strongly non-linear way.

Richards' equation is known to be difficult to solve using straightforward numerical methods that work well with most physical laws. It can be formulated in the pressure head,saturation, or mixed form. Classical approach is to discretize the pressure head form. This requires that the storage is averaged somehow between two time steps. However, any average gives large balance errors if the saturation changes too much from one time step to another, which is the case at the water table unless the time steps are very small. This problem is solved by discretizing the mixed form of Richards' equation using Celia et al. (1990) or Newton method.

Discretizing the mixed form of Richards' equation is not sufficient to obtain an accurate balance. Iterative algorithms used for unsaturated flow in the de facto industry standard modeling software typically do not converge. Groundwater modelers have no choice but to accept that particular cells along the water table oscillate between the saturated and unsaturated state and stop the computation, with a disastrous consequence on the total balance.

Convergence may be hampered for a number of reasons. However, one often overlooked issue which most drastically influences the convergence as well as the accuracy of unsaturated simulations is the method used to average the hydraulic conductivity. The conductivity dependence on the pressure head is modelled using an analytical formula such as van Genuchten or Gardner. Hydraulic conductivity is needed at points where the pressure head is not available. In finite volume methods the pressure head is typically directly available only at cell centers, while the conductivity must be evaluated between cells. Thus some kind of average must be computed. 

All classical averages ‒ arithmetic, geometric, or harmonic ‒ give wrong results in some cases as shown by Haverkamp and Vauclin (1979) which may also lead to convergence problems. A number of papers proposes using so-called Darcian means to average the conductivity: Baker (1999, 2000, 2006), Desbartes (1995), Warrick (1991), and Gasto et al. (2002). These methods give accurate results but are computationally expensive or have other limitations.

The reason why classical averages give so inaccurate results is that they lead to violation of pressure head maximum and minimum principle. Hydraulic head as well as the pressure head satisfying Richards' equation also satisfy the maximum and minimum principles, which means that inside the domain these variables cannot exceed the maximum and minimum values that they attain at the boundary or at the previous time step. Computing the average conductivity in a way that respects the maximum principle, as in Szymkiewicz (2009) and Misiats and Lipnikov (2013), is cheaper and more flexible then using approximate Darcian means. This is the option that we pursue.

Most literature on solving Richards' equation numerically considers the one-dimensional case only. Physical characteristics as well as numerical difficulties typical for Richards' equation are already pronounced in one dimension. However, unsaturated zone may be only a part of a larger model where the flow under the water table is essentially threedimensional. While using a very small time step is an option in one dimension, this is not practical in more complex three-dimensional cases. Thus in three dimensions an efficient numerical method is both necessary and more difficult to construct.

Finite volume methods present a good choice for unsaturated flow simulations because they conserve the mass locally. When they are applied to Richards' equation in the mixed form then the local mass balance error is proportional to the non-linear solver error which can be made arbitrarily small, subject only to the machine precision. Another consequence of local mass conservation is that finite volume methods are good in capturing steep hydraulic head gradients that often arise from Richards' equation.

This paper is organized as follows. First we discuss Richards' equation and other issues related to the physical model, followed by the list of boundary condition types typically used in groundwater engineering. Then we discuss the numerical method, paying a special attention to hydraulic conductivity averaging in one and three dimensions. Finally we provide numerical examples.