Time step and stability control for a coupled model of surface and subsurface flow

We discuss several numerical issues that arise in the integration of a coupled model of surface runoff and groundwater flow. The model is based on the combination of a finite element discretization of the 3D Richards' equation and a 1D diffusion wave formulation solved by a Muskingum-Cunge scheme with variable parameters. The effects of catchment topography and river network structure on stormflow response are incorporated through a DEM-based conceptual drainage network. Special treatment is used for the modeling of topographic depressions such as lakes and reservoirs. The surface routing and subsurface flow phenomena are physically coupled: overland flow and surface ponding are determined by and fed back into the subsurface-controlled partitioning of atmospheric forcing (rainfall and evaporation) into infiltration/exfiltration and runoff. This coupling is nonlinear, but the time-explicit Muskingum-Cunge discretization allows the construction of a noniterative algorithm for the solution of the resulting system of PDEs. The exchange of information in the coupling algorithm occurs through the source terms of the two PDEs. The resulting numerical algorithm raises several interesting numerical issues that will be examined in this study. For example, the highly nonlinear character of Richards' equation requires the use of an iterative Newton-like linearization scheme. The time step size control algorithm is crucial in guaranteeing global convergence and thus a successful simulation. This in turn is linked to the fact that the surface routing and subsurface flow processes act at different time scales, so that different time step sizes can be used for the time integration of the two equations. The Muskingum-Cunge scheme is subject to a stability constraint on both time step and DEM cell size that influences the performance of the coupled model in a complex manner. Numerical examples based on a real basin and with real data will be used to show the properties of the proposed numerical approach.