the standard deviation (ff). We determine a standard deviation value per hydrological year based on a resampled daily time series of a representative column in the recharge area of the Judean Mountains (i.e., near the Bar Giora site, in the district Jerusalem). This column was selected since this region with an average annual precipitation input of610 mm largely contributes to the overall recharge (i.e., it is relevant to the overall recharge dynamics) and at the same time exhibits a fairly thick vadose zone of about 335 m. We measure the annual signal standard deviation at all 66 unsaturated nodes of the vertical column.
The signal intensity of the recharge flux at the control plane groundwater table is largely controlled by the intensity of the precipitation input (see Fig. 3.11a). However, the thick vadose zone substantially contributes to the dissipation of input signal depending on the system state (see Fig. 3.11b), leading to high time-varying dynamics within the vadose zone. These results highlight the need for process-based modeling approaches since linear transfer functions are limited to time-invariant and linear responses.
As illustrated in Figure 3.10c, the vadose zone presently comprises circa 45 %of the dynamic storage volume, while under predevelopment conditions (i.e., before the 1950s when large-scale groundwater abstraction started), the vadose zone storage comprised circa 38 %of the dynamic storage. The increase in stored vadose water results from the substantial decline in groundwater levels since the 1950s and the consequent increase of the vadose zone depth. Dynamic storage refers to the stored water above the natural outlet (Taninim at 4 masl.), i.e., the water available to spring discharge. The dynamic storage is relevant for aquifer management since this water eventually discharges at its natural outlet and represents the available storage for groundwater abstraction. Figure 3.10c illustrates the storage dynamics after the exceptionally wet rainy season, with the vadose storage rapidly increasing and a long-term recession over roughly two years, despite the absolute vadose storage decreasing very quickly to levels before the wet year. This is a secondary effect of increasing water tables (i.e., mainly fed by pathways of fast infiltration) and a consequent decrease in the volume of the vadose zone. The results demonstrate the importance of accounting for vadose flow processes for managing semi-arid karst aquifers with thick vadose zones. Shifted water budgets because of climate change and increased anthropogenic consumption may further accentuate the relevance of the vadose zone for short-term storage and the dissipation of the hydraulic signal.
10
410
310
2i
(m
3s
1)
300 250 200 150 100 50 0
Depth below surface (m)
a)
1991/92 1998/99
10
210
110
0i
/
0(-)
b)
1991/92 1998/99
200 300 400 500 600 800 1000 1200
Pr ec ipi ta tio n he igh t ( m m a
1)
Figure 3.11: Dissipation of the input signal per hydrological year (i.e., from September to August) measured by a) the standard deviation (ffi) and b) the normalized standard deviation (ffi/ff0) at various depths below the surface. The standard deviation is normalized by the standard deviation
(ff0) of the input signal at the soil level.
unknown parameters increases, making it more challenging to obtain parameters via field measurements and define them unambiguously. In contrast, lumping physical processes or geological structures risks having oversimplified model structures that do not account for relevant physical processes or local heterogeneity in carbonate aquifers. Structural uncertainty may also be caused by insufficient consideration of the spatial and temporal variation of natural processes, such as rainfall variability (Sapriza-Azuri et al., 2015).
The problem to be solved determines the appropriate level of model complexity. For karst aquifers with very thick vadose zones and semi-arid climates, considering the duality of flow and vadose flow is imperative when assessing the effects of climate change or dealing with temporally resolved issues. This study partly alleviates the problem of parametric ambiguity by calibrating a lumped hydraulic conductivity and obtaining the properties of the second continuum by measurable hydraulic properties at the core scale (see Section 3.3.1).
However, nearly all numerical models incorporate a simplified representation of physical processes (i.e., process abstraction). Depending on the spatio-temporal scale of the individual problem, process abstraction is applied at a subordinate scale that may have considerable effects on the problem scale. The required degree of process consideration depends on the impact on the problem scale. It is clear that also the proposed model exhibits a simplified process representation that shall be described in the following. Unsaturated flow in fractured and karstified media occurs on various length and temporal scales, complicating an appropriate mathematical representation of the physical processes. A
complex interplay of capillary, gravitational, and viscous forces shape unsaturated flow (Dragila& Weisbrod,2004;Kordilla et al., 2017;Kordillaet al., 2021;Lehmann et al., 2012; Or, 2008; Shigorina et al., 2021) that may not be discretely resolved.
Unsaturated flow in fractured/cavernous media is particularly challenging, as highly dynamic processes govern flow, such as turbulent flow regimes, i.e., inertia-driven infiltration (e.g., droplet, rivulet, and film flow). This effect is negligible for laminar flows in porous media because Reynolds numbers are typically low. However, for conduits and fractures, turbulent flow regimes are the dominant regime of fluid motion. However, bulk-effective approaches, such as the Richards’ equation with VG parameterization, are subject to process abstraction by describing laminar flow based on capillarity and assuming instant hydrostatic equilibration.
Although this study assumes that the Richards’ equation is valid for the second continuum, which is a model built-in necessary simplification of the occurring processes, the proposed model considers physical processes with unprecedented detail. Volume-averaged dual- permeability approaches provide a sweet spot for the required computational resources and data demand for catchment scale problems. The application of discrete fracture/conduit numerical models is impractical presently on large-scale catchments, such as the WMA, as it requires the discrete definition of conduits and fractures. The proposed model circumvents the inaccurate process representation of the second continuum by applying a higher minimum relative conductivity that truncates the saturation-permeability curve by setting a lower permeability threshold. This represents the effect of the lower capillarity, as the pore space is nearly ubiquitous in the second continuum. Water retention in the ample void space of the second continuum is caused only by capillary barriers imposed through fracture intersections (Dragila &Weisbrod,2004;Wood et al.,2005) and surface roughness (Tokunaga &Wan, 1997; Tokunaga et al.,2000). Under inertia-driven conditions, the VG model cannot accurately model the occurring processes of advancing wetting fronts.
However, it is valid to approximate these dynamics at the lower saturation limit (i.e., dry conditions) with a minimum hydraulic conductivity. Applying composite VG parameters and a minimum relative conductivity is a heuristic approximation of the occurring physical processes.
It should be noted thatKordilla et al. (2012) applied a very similar approach to model groundwater flow in a karst aquifer. In contrast to this study, a simpler model geometry (i.e., 2D ), temperate climate with moderate variability of precipitation inputs (i.e., moderate change in state variables), and lower vadose zone thickness constitute the main differences. The VG material model as applied byKordillaet al.,2012leads to the inability of preferential infiltration through the second continuum due to high negative pressures of the matrix continuum under dry conditions. Kordillaet al.,2012 circumvented this shortcoming by directly adding a predefined fraction of the net infiltration to the groundwater table. In contrast, this study applies a composite VG material model for both continua. In
addition, we truncate the VG material model of the second continuum using a minimum relative conductivity to allow for preferential infiltration through the second continuum under dry conditions. Consequently, the applied method still describes the flow in terms of capillary and gravitational forces. However, we approximate inertia-driven infiltration heuristically with the minimum relative conductivity, which is a necessary abstraction of the occurring processes, given the scale of investigation and unknown dimensions of conduits and fractures.
Also, other bulk-effective and discrete approaches have been developed recently to deal with fractured-porous systems. Two-dimensional analytical and numerical solutions such as those provided by Moutsopoulos (2021) are suitable to simulate unconfined dual-permeability systems. However, they cannot capture all relevant processes due to the three-dimensional setting and heterogeneous distribution of hydrogeological and hydrological properties within the WMA. Discrete modeling approaches, such as smoothed particle hydrodynamics (SPH) simulations (Kordilla et al., 2017; Shigorina et al., 2021) are a promising approach to improve the physical representation of inertia-driven infiltration processes. Shigorina et al. (2021) combines SPH for inertia-driven flow through fractures with a bulk-effective approach based on Richards’ equation for the unsaturated flow through the matrix. However, those approaches demand high computation power and are currently limited to the scale of individual fractures or a series of a few fractures. In addition, other discrete approaches, such as discrete fracture network models (Hymanet al.,2015) can often not be implemented on a catchment scale due to the missing information about their location, geometry, and extent of individual conduits or fractures. Here, bulk-effective approaches such as dual permeability models reconcile low data requirements with accurate physical representation.
Furthermore, this study applies a dual-medium water balance model to estimate the partitioning of precipitation into evapotranspiration and infiltration. This approach has some advantages compared to the previously commonly applied linear regression models (e.g., Abusaada, 2011; Berger, 1999; Weiss & Gvirtzman, 2007), such as that it heuristically accounts for the duality of infiltration and considers the temporal variability of precipitation. However, the model accounts for the occurring processes in a lumped physical manner, and some processes are simplified. Applying the globally-implicit integrated surface- subsurface approach within HGS would further improve the predictive capabilities since it accounts for depth-integrated surface flow and the natural partitioning of precipitation into runoff, evapotranspiration, and infiltration. However, currently employed solvers and parallelization schemes of HGS due to computational limitations (i.e., maximum speed up of 2:3by addition of parallel resources; see Section 3.2.3) currently do not allow implementing surface processes on top of an already highly complex subsurface flow model. Thus, given
the current parallelization scheme, numerical implementations, and available computational resources, we believe in having reached the present-day maximum conceivable process representation for a catchment of this scale and vadose zone thickness.
Preexisting catchment scale models of the WMA (e.g., Abusaada, 2011; Bach- mat,1995;Bachmat & Wolman,2000;Berger,1999;Dafny et al.,2010;Guttman
& Zukerman,1995) neglect the effect of a several hundred-meter thick vadose zone on infiltration dynamics and apply parameters with lumped physical meaning, e.g., indirect rep- resentation of vadose zone flow and the separating aquitard via highly anisotropic hydraulic conductivities. Also, they apply simple linear correlations of precipitation and recharge, neglecting the spatial and temporal variation of recharge and the storage effect of the vadose zone. Therefore, these models are subject to higher structural uncertainty (e.g., caused by the non-physical representation of the infiltration process) than the proposed model and are more likely prone to parameter degeneration (i.e., non-physical values). Furthermore, as in many Mediterranean karst aquifers, the unsaturated storage capacity may vary significantly with time and space for aquifers with vadose zones of several hundred-meters. Therefore, it is necessary to consider the superposition of rapid recharge via fractures and dissolution features (e.g., dolines, shafts) and long-term diffuse infiltration via the matrix domain, particularly in the light of the imminent changing precipitation patterns due to climate change effects. Investigating the effect of changing rainfall patterns because of shifting climates requires a process-based approach that accounts for as many relevant physical processes as possible (Brouyèreet al., 2004; Cuthbert&Tindimugaya, 2010; Moeck et al., 2018). Further, the dual-permeability model provides more realistic estimates of the global hydraulic parameters due to the lower structural uncertainty (at the expense of computational time).
In this study, we have demonstrated that variably saturated flow in a thick vadose zone has a considerable impact on the dissipation of the hydraulic signal on the catchment scale. The presented model expands the predictive options compared to previously employed approaches by accounting for fully-distributed (i.e., time and space-dependent) vadose storage dynamics. With the ability to simulate rapid infiltration, the change in storage in the vadose and the phreatic zone, and the characteristics of the karst system dynamics, the constructed flow model has considerable advantages compared to currently available models allowing to predict spatio-temporally distributed recharge and infiltration. In addition, the non-linear response of groundwater recharge to transient climatic inputs requires a daily resolution for the analysis. Lastly, the flow model provides insight into the infiltration dynamics at the catchment scale, i.e., mean residence times in the vadose zone characterized by a potential-dependent exchange between slow/diffuse and fast flow systems. Hence, it predicts the ability to control the long-term release of water.
Consequently, the model can predict groundwater flow and infiltration dynamics under changing climates with extremes that expand beyond previously experienced conditions.
For example, climate change simulations project that average long-term precipitation will decrease, yet the intensity and frequency of short-duration extreme rainfall events will increase (Adinolfi et al., 2020; Bucchignani et al., 2018; Hochman et al., 2018).
Together with increased potential evapotranspiration, it will affect the water budget of (Mediterranean) karst aquifers. Accordingly, the proposed model may predict available water resources under climate change and support decision-makers in planning the required scaling of more costly alternative water sources, such as desalinated water and treated greywater. In addition, the model can be used to assess the potential for managed aquifer recharge applications (i.e., infiltration ponds, dry wells) to mitigate the effects of drought years. The proposed model may also assist in developing more elaborate transfer functions for vadose storage dynamics since it provides fully-distributed information about the storage in the vadose zone. Thus, it may aid in demonstrating weaknesses of process abstraction and improving less complex numerical models that are more convenient in their application for water management authorities.
References
Abusaada, M. J. (2011). “Flow Dynamics and Management Options in a Stressed Carbonate Aquifer System, The Western Aquifer Basin, Palestine”. PhD thesis. Georg-August- Universität Göttingen, p. 182. url:http://hdl.handle.net/11858/00-1735-0000- 0006-B2FD-2(visited on 11/10/2017).
Abusaada, M. J. & M. Sauter (2017). “Recharge Estimation in Karst Aquifers by Applying Water Level Fluctuation Approach”. In: Int. J. Earth Sci. Geophys.3.1. doi:
10.35840/2631-5033/1813.
Adinolfi, M., E. Bucchignani & P. Mercogliano (2020). High-resolution climate simulation over Israel performed with COSMO-CLM and evaluation of results. Tech. rep.
Euro-Mediterranean Center on Climate Change, pp. 1–24.
Aguilar-López, J. P., T.Bogaard & H. H. Gerke(2020). “Dual-Permeability Model Improvements for Representation of Preferential Flow in Fractured Clays”. In: Water Resour. Res. 56.8, pp. 1–20.doi:10.1029/2020WR027304.
Allen, R. G., L. S. Pereira, D. Raes & M.Smith (1998).Crop evapotranspiration - Guidelines for computing crop water requirements. Food and Agriculture Organization of the United Nations. isbn: 92-5-104219-5.url:http://www.fao.org/3/X0490E/
x0490e00.htm (visited on 10/12/2022).
Aquanty Inc.(2015).HydroGeoSphere User Manual - Theory.url:https://aquanty- artifacts-public.s3.amazonaws.com/hgs/hydrosphere_theory.pdf (visited on 10/08/2022).
Bachmat, Y. (1995). Hydrologic model of the Western Mountain groundwater basin. Tech. rep. Harvard Middle East Water Project.
Bachmat, Y. & S.Wolman (2000). “Development and Testing of a Hydrological Model to Assess the movement of Pollutants in the Western Basin of the Mountain Aquifer”.
In:Environmental Protection of the Shared Israeli-Palestinian Mountain Aquifer. Final Report. Jan. 1994 – Dec. 1999. USAID, pp. 89–138.
Bank, R. E. & J. Xu (1996). “An algorithm for coarsening unstructured meshes”. In:
Numer. Math. 73.1, pp. 1–36.doi:10.1007/s002110050181.
Banusch, S., M.Somogyvári, M.Sauter, P.Renard& I.Engelhardt(2022). “Stochas- tic Modeling Approach to Identify Uncertainties of Karst Conduit Networks in Carbonate Aquifers”. In:Water Resour. Res.58.8. doi:10.1029/2021WR031710.
Berger, D. (1999).Hydrological model for the Yarqon-Taninim Aquifer. Tech. rep. Mekorot, p. 50.
Braun, M. & F.Hirsch(1994). “Mid Cretaceus (Albian-Cenomanian) carbonate platforms in Israel”. In: J. Iber. Geol. 18, pp. 59–82.
Brooks, R. H. & A. T. Corey(1964). “Hydraulic Properties of Porous Media and Their Relation to Drainage Design”. In: Trans. ASABE 7.1, pp. 26–28.doi:10.13031/2013.
40684.
Brouyère, S., G. Carabin & A. Dassargues (2004). “Climate change impacts on groundwater resources: modelled deficits in a chalky aquifer, Geer basin, Belgium”. In:
Hydrol. J. 12.2, pp. 123–134.doi:10.1007/s10040-003-0293-1.
Bucchignani, E., P. Mercogliano, H.-J.Panitz& M.Montesarchio (2018). “Climate change projections for the Middle East–North Africa domain with COSMO-CLM at different spatial resolutions”. In: Adv. Clim. Change Res.9, pp. 66–80.doi: 10.1016/j.
accre.2018.01.004.
Celia, M. A., E. T. Bouloutas & R. L. Zarba (1990). “A general mass-conservative numerical solution for the unsaturated flow equation”. In: Water Resour. Res. 26.7, pp. 1483–1496. doi:10.1029/WR026i007p01483.
Contractor, D. N. & J. W. Jenson (2000). “Simulated effect of vadose infiltration on water levels in the Northern Guam Lens Aquifer”. In:J. Hydrol. 229.3-4, pp. 232–254.
doi:10.1016/S0022-1694(00)00157-8.
Cuthbert, M. O. & C. Tindimugaya (2010). “The importance of preferential flow in controlling groundwater recharge in tropical Africa and implications for modelling the impact of climate change on groundwater resources”. In: J. Water Clim. Change 1.4, pp. 234–245. doi:10.2166/wcc.2010.040.
Dafny, E., A. Burg& H.Gvirtzman (2010). “Effects of Karst and geological structure on groundwater flow: The case of Yarqon-Taninim Aquifer, Israel”. In: J. Hydrol.389.3-4, pp. 260–275. doi:10.1016/j.jhydrol.2010.05.038.
Dragila, M. I. & N. Weisbrod(2004). “Fluid motion through an unsaturated fracture junction”. In: Water Resour. Res.40.2, pp. 1–11.doi:10.1029/2003WR002588. Dvory, N. Z., M. Kuznetsov, Y. Livshitz, G. Gasser, I. Pankratov, O. Lev, E.
Adar& A. Yakirevich(2018). “Modeling sewage leakage and transport in carbonate aquifer using carbamazepine as an indicator”. In: Water Res. 128, pp. 157–170. doi:
10.1016/j.watres.2017.10.044.
Dvory, N. Z., Y. Livshitz, M. Kuznetsov, E. Adar& A. Yakirevich (2016). “The effect of hydrogeological conditions on variability and dynamic of groundwater recharge in a carbonate aquifer at local scale”. In: J. Hydrol. 535, pp. 480–494.doi: 10.1016/j.
jhydrol.2016.02.011.
Fabryka-Martin, J., P. Dixon, S. Levy, B. Liu, H.Turin & A. Wolfsberg(1996).
Summary report on chlorine-36 studies: Systematic sampling for chlorine-36 in the exploratory studies facility. Tech. rep. 3783AD. Los Alamos National Laboratory.
Farthing, M. W. & F. L.Ogden (2017). “Numerical Solution of Richards’ Equation: A Review of Advances and Challenges”. In: Soil Sci. Soc. Am. J.81.6, pp. 1257–1269.doi:
10.2136/sssaj2017.02.0058.
Fischhendler, I. & A. Frumkin (2008). “Distribution, evolution, and morphology of caves in southwestern Samaria, Israel”. In: Isr. J. Earth Sci. 57.3, pp. 311–322. doi:
10.1560/IJES.57.3-4.311.
Flecker, R. & R. M. Ellam (1999). “Distinguishing climatic and tectonic signals in the sedimentary successions of marginal basins using Sr isotopes: an example from the Messinian salinity crisis, Eastern Mediterranean”. In: J. Geol. Soc. London 156.4, pp. 847–854. doi:10.1144/gsjgs.156.4.0847.
Fleischer, L. (2002).Stratigraphic Table of Israel: Outcrops and Subsurface. url:https:
/ / www . gov . il / BlobFolder / generalpage / stratigraph - map - of - israel / en / maps_28-5_stratigraphic_table.pdf (visited on 10/26/2020).
Frumkin, A. & H.Gvirtzman(2006). “Cross-formational rising groundwater at an artesian karstic basin: the Ayalon Saline Anomaly, Israel”. In: J. Hydrol.318, pp. 316–333.doi:
10.1016/j.jhydrol.2005.06.026.
Gerke, H. H. & M. T.van Genuchten(1993a). “A Dual-Porosity Model for Simulating the Preferential Movement of Water and Solutes in Structured Porous Media”. In: Water Resour. Res. 29.2, pp. 305–319.doi:10.1029/92WR02339.
– (1993b). “Evaluation of a first-order water transfer term for variably saturated dual- porosity flow models”. In: Water Resour. Res. 29.4, pp. 1225–1238. doi: 10 . 1029 / 92WR02467.
Geyer, T., S.Birk, R.Liedl& M.Sauter(2008). “Quantification of temporal distribution of recharge in karst systems from spring hydrographs”. In: J. Hydrol.348, pp. 452–463.
doi:10.1016/j.jhydrol.2007.10.015.
Guttman, J. & H.Zukerman(1995).Yarkon Taninim - Be’er Sheva Basin. Establishment and calibration of a water-flow and solute-transport model. Tech. rep. 01/95/72. Tahal Consulting Engineers Ltd.
Hartmann, A., J.Lange, À.Vivó Aguado, N.Mizyed, G. Smiatek& H.Kunstmann (2012). “A multi-model approach for improved simulations of future water availability at a large Eastern Mediterranean karst spring”. In: J. Hydrol.468-469, pp. 130–138.doi:
10.1016/j.jhydrol.2012.08.024.
Hochman, A., P.Mercogliano, P.Alpert, H.Saaroni& E.Bucchignani(2018). “High- resolution projection of climate change and extremity over Israel using COSMO-CLM”.
In: Int. J. Climatol. 38.14, pp. 5095–5106. doi:10.1002/joc.5714.
Hwang, H.-T., Y.-J.Park, E.Sudicky & P.Forsyth (2014). “A parallel computational framework to solve flow and transport in integrated surface–subsurface hydrologic systems”.
In: Environ. Model. Softw.61, pp. 39–58.doi:10.1016/j.envsoft.2014.06.024. Hyman, J. D., S. Karra, N. Makedonska, C. W. Gable, S. L. Painter & H. S.
Viswanathan(2015). “DfnWorks: A discrete fracture network framework for modeling subsurface flow and transport”. In: Comput. Geosci. 84, pp. 10–19. doi:10.1016/j.
cageo.2015.08.001.
Jeannin, P.-Y. & M. Sauter (1998). “Analysis of karst hydrodynamic behaviour using global approaches: a review”. In: Bull. Hydrol.16.16, pp. 31–48.
Kordilla, J., T.Noffz, M.Dentz, T.Geyer & A. M.Tartakovsky(2017). “Effect of Unsaturated Flow Modes on Partitioning Dynamics of Gravity-Driven Flow at a Simple Fracture Intersection: Laboratory Study and Three-Dimensional Smoothed Particle Hydrodynamics Simulations”. In: Water Resour. Res. 53.11, pp. 9496–9518. doi: 10.
1002/2016WR020236.
Kordilla, J., M.Sauter, T.Reimann & T.Geyer(2012). “Simulation of saturated and unsaturated flow in karst systems at catchment scale using a double continuum approach”.
In: Hydrol. Earth Syst. Sci. 16.10, pp. 3909–3923.doi: 10.5194/hess-16-3909-2012.
Kordilla, J., M. Dentz & A. M. Tartakovsky (2021). “Numerical and Analytical Modeling of Flow Partitioning in Partially Saturated Fracture Networks”. In: Water Resour. Res. 57.4. doi:10.1029/2020WR028775.
Krijgsman, W., F. J.Hilgen, I.Raffi, F. J.Sierro& D. S.Wilson(1999). “Chronology, causes and progression of the Messinian salinity crisis”. In:Nature 400.6745, pp. 652–655.
doi:10.1038/23231.
Laskow, M., M.Gendler, I. Goldberg, H.Gvirtzman& A. Frumkin (2011). “Deep confined karst detection, analysis and paleo-hydrology reconstruction at a basin-wide scale using new geophysical interpretation of borehole logs”. In: J. Hydrol. 406.3-4, pp. 158–169. doi:10.1016/j.jhydrol.2011.06.011.
