Higher Order Time Integration and Discontinuous Galerkin Methods for Variably Saturated Groundwater Flow

The numerical simulation of groundwater flow in the vadose zone continues to be a challenge for many problems of practical interest. Under commonly used constitutive relations, the governing equations can be highly nonlinear and produce sharp fronts in the solution variables for problems such as wetting phase infiltration into an initially dry medium. For a number of multiphase flow problems, the use of variable order and variable step size temporal discretizations has shown significant advantages. However, the spatial discretizations commonly used for variably saturated flow are dominated by low-order finite difference and finite element methods. Over the last decade discontinuous Galerkin (DG) finite element methods have received significant attention in a number of fields for hyperbolic PDE's and, more recently, for elliptic and parabolic problems. DG approaches are appealing for modeling subsurface flow since they can lead to velocity fields that are locally mass-conserving without the need for auxiliary variables or alternative meshes. Moreover, DG discretizations are inherently local and so well-suited for unstructured meshes and h-p adaption strategies. While some work has been done recently for multiphase subsurface flow, there are a range of issues related to the performance of DG methods for highly nonlinear parabolic problems that have not been investigated fully, particularly for air-water systems. In this work, we consider the combination of higher order adaptive time integration with DG spatial discretizations applied to variably saturated groundwater flow. We compare this approach to standard low order methods for a series of test problems and consider a number of issues including the methods' relative accuracy and computational efficiency.