Highly heterogeneous flow dynamics characterize karst aquifers in all compart- ments, i.e., the land surface, the vadose, and the phreatic zone. The heterogeneous distribution of hydraulic properties (e.g., hydraulic conductivity, storage) within the com- partments affects the local and global routing dynamics and hence the magnitude and dissipation of the hydraulic signal between the source (precipitation and recharge) and spring outlet (Jeannin &Sauter, 1998;Smart&Hobbs,1986). In particular, under semi-arid climates with pronounced precipitation seasonality (i.e., dry summer and humid winter), a thick vadose zone may provide substantial storage to infiltrating water (Geyer et al., 2008;Liuet al., 2021;Schmidtet al., 2014) and storage may vary significantly in time and space. However, the role of the vadose zone is often neglected or only accounted
for by employing simple (linear) transfer functions in catchment scale numerical simulations.
Furthermore, a complex interplay between non-linear processes with superimposing effects shapes storage in the vadose zone, i.e., the initial state of saturation controls matric suction and storage capacity, requiring process-based approaches to account for the distributed infiltration through the vadose zone.
- 419 Elevation (masl.)
Flow direction Spring Pumping well
Recharge area WMA (study area)
Mediter EMA ranea n Sea
Figure 3.1: Location and extent of the Western Mountain Aquifer (WMA) and the adjacent Eastern Mountain Aquifer (EMA) catchment (Digital elevation data obtained from the data set SRTM 1 Arc-Second Global).
This study demonstrates a complex process-based method, a variably saturated dual-permeability flow model, to assess infiltration dynamics and storage in the vadose zone for the Western Mountain Aquifer (WMA) in Israel and the West Bank. The reasoning for selecting this aquifer is that it provides (1) a highly complex karst environment including multi-domain heterogeneities in the phreatic and vadose zone, (2) a large catchment size with effects of distributive hydrological input, and (3) at the same time the integrated hydraulic response at the Taninim spring which reflects the entire catchment dynamics. In addition, semi-arid climates exhibit a pronounced seasonality of precipitation and intense short-duration rainfalls and, therefore, can provide clear and pronounced recharge input signals not superimposed by summer rainstorms. The unambiguous signal allows for analyzing the contribution of the vadose zone to the delayed dissipation of the hydraulic discharge response to rainstorm events. The large catchment size of the WMA with 9000 km2 ascertains that small-scale heterogeneity in hydraulic properties and spatial variability in the distribution of precipitation are averaged out in the integrated response at the spring, i.e., the robustness of large-scale model results concerning local hydraulic
parameter variations and data uncertainty. Also, the WMA is highly relevant to the water supply of Israel and the West Bank. A process-based approach improves model validity in long-term predictive simulations. Furthermore, it allows for the extraction of information about the time-variant moisture distribution within the thick vadose zone, hence providing a tool to improve groundwater management strategies for highly stressed aquifer systems.
For instance, Figure 3.2 shows the effect of very wet years (e.g., 1973/74, 1991/92) on the long-term recession of spring discharge at the Auja spring (West Bank) and groundwater hydrographs in the recharge area. Despite the Auja spring being located outside the WMA in the Jordan Valley and belonging to the adjacent Eastern Mountain Aquifer (EMA), the hydraulic response of the Auja spring can be considered representative of the WMA since the WMA and EMA share a common recharge area and similar geological units.
However, the signal in the EMA is not perturbed by the effects of intensive groundwater abstraction. High discharge and groundwater tables were recorded for several years even though recharge was considerably less. This observation implies groundwater recharge is stored in either the vadose or the phreatic zone. Therefore, the most appropriate method to investigate the coupled process of infiltration, recharge, and groundwater flow is the discrete simulation of the hydraulically coupled phreatic and vadose zone, considering their respective duality in flow behavior. The superposition of rapid recharge via fractures and dissolution features (i.e., dolines, shafts) and long-term diffuse infiltration via the limestone rock matrix domain must be accounted for when addressing the storage dynamics.
1P (mm a) 1000 608 640 465 765 581 482 636 573 449 785 628 529 814 426 478 449 693 707 527 613 490 1169 653 442 723 579 646 583 322 511 506
a) Seasonal precipitation
1972 1974 1976 1978 1980 1982 1984 1986 1988 1990 1992 1994 1996 1998 2000 0.0
0.2 0.4 0.6
Spring discharge (m3 s1) b) Hydraulic response
290 300 310
Hydraulic head (masl.)
Figure 3.2: Comparison of a) precipitation input and b) observed hydraulic response of the Auja spring (West Bank) discharge and the hydraulic head in the recharge area (Spring discharge data obtained fromSchmidt et al.,
In the following, we briefly discuss the existing groundwater flow and recharge models of the WMA and describe the applications and merits of the presented model.
Presently existing flow models for the WMA calculate fully saturated flow and parameterize
the model by employing a single bulk-effective continuum (i.e., a lumped matrix/fracture- conduit continuum) (e.g., Abusaada,2011;Berger, 1999;Dafnyet al.,2010). However, a saturated single continuum flow model can not account for the storage effect of a several hundred-meter thick vadose zone and preferential infiltration pathways. Moreover, a dual-permeability model can better reproduce the temporal dynamics of event-based spring discharge and local groundwater table responses, characterized by a rapid rise after rainfall and gradual and extended recessions. For instance,Abusaada&Sauter(2017) represents the vertical aquifer layer structure spatially lumped by applying an accentuated hydraulic anisotropy, where the low vertical hydraulic conductivity reflects aquitards and the vadose zone. The appropriate choice of process consideration depends on the study objective, prevailing flow processes, and site-specific characteristics. This study aims to provide a tool to assess distributed vadose and phreatic storage dynamics on a regional scale and therefore requires the consideration of vadose flow processes. Vadose flow processes gain further importance when the numerical model is applied to assess the effects of climate change and unsustainable groundwater abstraction, likely scaling up the temporal delay because of a thicker vadose zone. In addition, recharge to the WMA is quantified in most previous studies by linear regression models between annual precipitation or the precipitation in the recharge effective winter months and annual recharge (Abusaada,2011;Berger, 1999;Weiss& Gvirtzman,2007), neglecting the intra-annual variability of rainfall. For instance,Abusaada (2011) computed annual recharge (i.e., defined as input at the level of the groundwater table) based on a linear regression between precipitation during the recharge-effective winter months and subsequently accounted for the transfer through the vadose zone by an empirically obtained normalized recharge curve on a monthly time step.
The recharge curve is the median seasonal water table fluctuation over multiple years.
Dafnyet al. (2010) define the recharge flux at a fixed fraction of the annual precipitation.
While the above methods have their merits (e.g., simplicity, parsimony of applied parameters, and lesser demand on data), they do not account for the spatio-temporal distribution of precipitation, recharge, and infiltration. Therefore, they can not be employed to predict future recharge and storage dynamics for different climatic conditions. Sheffer et al.
(2010) provide an approach to assess distributed recharge to the WMA based on a daily soil water balance model for three soil types (i.e., Terra Rossa, Rendzina, and desert) coupled with a monthly groundwater flow model (fully saturated). However, this approach neglects the storage and delaying effect of the vadose zones and the dual-domain flow dynamics. Furthermore, the recharge area of the WMA is characterized by a large proportion of exposed bare carbonate rock, providing potential pathways for rapid infiltration into the consolidated formations. Dvory et al.,2018,2016 have considered the vadose flow processes at the local scale by employing a 2-dimensional dual-permeability infiltration model for a smaller sub-catchment of the WMA (i.e., the Soreq catchment).
This study accounts for infiltration and vadose storage dynamics at the catchment scale by employing a variably saturated dual-permeability flow model. We discretize the rock matrix and secondary porosity (i.e., conduits and fractures) as separate overlapping continua. Flow is computed via the bulk-effective Richards’ equation with van Genuchten (VG) parameters to model flow in karstified limestone rock material, with rapid flow through conduits and slow flow through the rock matrix in the vadose and the phreatic zone. The individual continua transfer water based on a first-order exchange term across a coupling interface with a coupling length and a conductance as controlling parameters. The applied model accounts for short-term recharge, essential in climate zones characterized by pronounced seasonality (i.e., semi-arid climates). Furthermore, preferential flow paths enhance accentuated recharge through the vadose zone. Climate simulations for Israel and West Bank suggest an increased frequency of extreme rain events and an overall decrease in annual precipitation because of climate change (Adinolfi et al.,2020; Bucchignaniet al., 2018; Hochman et al., 2018), further emphasizing dynamic infiltration characteristics, their combined effects on vadose zone storage, and the necessity to account for appropriate modeling techniques.