Distributed modeling of groundwater recharge in a semi-arid carbonate aquifer using the integrated surface-subsurface flow simulator HydroGeoSphere

Present methods to quantify recharge in karst aquifers in many cases rely on spatially and temporally aggregated precipitation values, neglecting the highly erratic, non-linear nature of infiltration dynamics that give rise to a dual-domain behavior with a slow diffuse and fast focused recharge component. Here, we demonstrate the applicability of integrated surface-subsurface flow models to simulate diffuse and preferential infiltration within the large scale Western-Mountain-Aquifer (WMA) in Israel and the Palestinian territories. A semi-arid climate region with a highly pronounced seasonality of precipitation and intense short-duration rainfalls, such as the Mediterranean region, emphasizes the importance of understanding and accounting for the complex dynamics of dual-domain infiltration and partitioning of the precipitation input signal via spatially discretized overland flow processes.

We apply HydroGeoSphere as a dual-continuum flow simulator for transient variably-saturated water flows, discretizing the rock matrix and secondary porosity (i.e., conduits and fractures) as separate overlapping continua. Flow is respectively computed via the Richards' equation with volume-averaged van Genuchten parameters, assuming that the Richards' equation is valid for both domains. The presented model accounts for surface flow via the two-dimensional Saint-Vénant equations under nonexistent inertial forces. We apply precipitation directly to the overland flow continuum and naturally account for the partitioning into Horton overland flow and percolating water. However, modeling of unsaturated flow through the conduit/fracture continuum with the van Genuchten parameterization is often limited, as the Richards' equation describes flow solely in terms of capillary forces, leading to high matric suction in the matrix continuum as a result of the smaller pore spaces (and hence constant exchange from the fracture continuum to the matrix system). In a natural system, non-linear transfer processes govern the transfer between fracture/conduit and matrix flow, such as inertia-driven infiltration (i.e., droplet, rivulet, and film flow) that initially retains itself from equilibration of capillary pressure heads and avoids instant matrix imbibition. This study demonstrates parametrization strategies to allow for infiltration through the fracture/conduit continuum using small-scale process-based simulations. Further, we offer procedures that help to achieve convergence of complex catchment-scale variably-saturated simulations.