Fill and Spill Hillslope Runoff Representation With a Richards Equation-Based Model

Integrated surface-subsurface hydrological models (ISSHMs) are well established numerical tools to investigate water flow and contaminant transport processes over a wide range of spatial and temporal scales. However, their ability to correctly reproduce the response of hydrological systems to natural and anthropogenic forcing depends largely on the accuracy of model parameterization, including the level of detail in the representation of the bedrock. This latter is typically incorporated in some way via the bottom boundary of the model domain. Issues of bedrock topography, variable soil depth, and the resulting hillslope storage distribution representation in ISSHMs are vitally important but to date have received little attention. A standard treatment of the bottom boundary, especially in large catchment and continental scale applications, is to model it as a flat or inclined (e.g., parallel to the surface) impermeable base (sometimes with some simple leakage term). This approach does not allow the model to correctly reproduce bedrock-controlled threshold responses such as the fill and spill process, as observed across many hillslope and catchment scale field studies. It is still unclear whether Richards equation-based numerical models are actually able to generate such responses. Here we use a Richards equation-based model (CATHY) to simulate internal transient subsurface stormflow dynamics observed at the well-characterized Panola experimental hillslope in Georgia (USA). Soil and bedrock properties were calibrated starting from values reported in previous studies at the site. Our simulation results show that the model was able to reproduce threshold mechanisms, which in turn affected both the integrated and distributed hydrologic responses of the Panola hillslope. We then developed a set of virtual experiments with modified boundary conditions and base topography at the soil-bedrock interface to explore the bedrock boundary control on transient groundwater flow patterns. Our results show that accurate representation of the lower boundary is crucial for ISSHM simulations of hillslope-scale storm runoff and for connectivity of transient groundwater. We summarize our findings with the development of a new bedrock topographic wetness index that takes into account the unsaturated infiltration dynamics. The index is able to help represent the spatial variability of water table response over the bedrock surface compared to standard surface topography-based indices. This new index may be useful in larger-scale ISSHM applications where an exact bedrock topography representation is not feasible or possible.