Quantification of Groundwater Recharge from an Ephemeral Stream into a Mountainous Karst Aquifer

Sustainable groundwater production from karst aquifers is primarily dictated by its recharge rate. Therefore, it is essential to accurately quantify annual groundwater recharge in order to limit overexploitation and to evaluate artificial methods for groundwater enrichment. Infiltration during erratic flood events in karst basins may substantially contribute to aquifer recharge. However, the complicated nature of karst systems, which are characterized in part by multiple springs, sinkholes, and losing/gaining streams, impede accurate quantification of the actual contribution of flood waters to groundwater recharge. In this study, we aim to quantify the proportion of groundwater recharge accrued during runoff events in a karst aquifer. The role of karst conduits on flash flood infiltration was examined during four flood and controlled runoff events in the Soreq creek near Jerusalem, Israel. We distinguished between direct infiltration, percolation through karst conduits, and diffuse infiltration—the latter of which is most affected by evapotranspiration. A water balance was calculated for the 2014/15 hydrological year using the Hydrologic Engineering Center-Hydrologic Modelling System (HEC-HMS). Simulations show that 6.8 to 19.2% of the annual recharge volume was added to the aquifer from infiltration of runoff losses along the creek through the karst system.


Introduction
Sustainable groundwater production from karstic aquifers is primarily dictated by their recharge rates.Therefore, in order to limit overexploitation, it is essential to accurately quantify groundwater recharge [1,2].Infiltration during erratic floods in karstic basins may contribute substantial amounts to aquifer recharge [1,[3][4][5].However, the complicated nature of karst systems presents a large obstacle to accurately predicting the true contribution of flood waters to groundwater recharge.
Traditionally, karst systems created over time by karstification processes are characterized in part by multiple springs, sinkholes, and losing/gaining streams [6,7].The complex nature of karst systems and their structural heterogeneity subsequently create complicated hydrodynamics unique to each system's storage and drainage structures [8].These intricacies heavily influences the surface water-groundwater interactions, particularly during flooding events [9,10].
Water 2018, 10, 79 2 of 16 Typically, the knowledge available about karst system structure (cavities, storage), their distribution within the groundwater system, their hydraulic conductivities, and storage coefficients, proves inadequate when modeling karst hydrodynamics [11].Therefore, many researchers have employed loss models, such as the: "initial and constant-rate loss" model, "deficit and constant-rate" model, "SCS curve number (CN) loss" model [12], and "Green and Ampt loss" model [13].Loss models are group of models that are used to estimate the volume of runoff, given the precipitation and properties of the watershed.In such models, losses can also occur by infiltration [12].In order to account for fast infiltration processes using these models, several studies employed spatially distributed infiltration coefficients for karst conduits [14,15].Other studies used different methods to represent karst systems, such as simulating sinkholes with ponds with a small drainage area and a large hydraulic conductivity at the bottom of the pond, or representing losing streams by tributary channels with high hydraulic conductivity in the stream bed [6].Refs.[16][17][18] used a deep aquifer percolation coefficient [19] in order to remove water from the watershed system.Another method used to estimate fast infiltration processes is by performing losing/gaining reaches analysis [11,20].
Another important factor affecting groundwater balance is evaporation [3,21].However, for the most part, studies which previously sought to quantify flood water infiltration from ephemeral streams, misestimated evaporation by insufficiently applying evaporation coefficients, or overlooking the differing influence of evaporation on distinct flow paths.Other studies addressed evaporation by making a valid assumption that evaporation was negligible due to the brevity of the flooding event [11,22,23], or did not differentiate between evaporation and percolation losses [24].In some studies, although evaporation coefficients were used, they were applied for time periods insufficient for conducting a yearly water balance analysis [25].Therefore, the current study employed a combined water balance approach, using multiple sources to quantify basin gains and year-round losses.
The objective of this study was to quantify the proportion of groundwater recharge from runoff events in a karst aquifer.In this study, we examine the role of karst conduits on flash flood infiltration during four flood and controlled runoff (by dam water release) events in the Soreq creek near Jerusalem.The runoff was monitored in short time steps to allow for an accurate estimate of surface flow volumes, which is of particular importance in semi-arid climates, where ephemeral flow could substantially contribute to groundwater sources [26].This paper differentiates between "direct" infiltration, i.e., percolation during runoff events through karst conduits, and "diffuse" infiltration, i.e., precipitation through hillslopes which can be notably affected by evapotranspiration.A water balance was then calculated for the hydrologic year 2014/2015.

The Soreq Creek and Floods Events
The study area is located on the Soreq creek, an ephemeral stream that runs westward of the city of Jerusalem, Israel (Figure 1).At this site, the Soreq creek drains a catchment area of approximately 99 km 2 .The upper part of this catchment area (78 km 2 ) drains runoff water into a reservoir (Beit Zait, Israel) that limits downstream natural surface flow.The channel length from the reservoir up to the study area outlet is 7.25 km, with a mean channel slope of 1.2%.The hillslopes in this part of the catchment have a mean slope of 29%.The mean annual precipitation in this area, according to the IMS Jerusalem meteorological station for 1980-2010, is 537 mm.Over 70% of the annual rainfall occurs during December-February.Although the stream has no base flow, the area experiences intense rainfall events which can cause the reservoir dam to overflow, which then develops into floods downstream.Such events occur every few rainy seasons, e.g., according to readily available data, between 2004/05 and 2014/15, the dam overflowed five times in three hydrological years.The floods are characterized by narrow flood peaks with short time lags.During the spring season, water is regularly released from the dam downstream.The study area is located within the recharge area of the fractured carbonate Western Mountain aquifer (Yarkon-Taninim) [27].The aquifer extends from south of the Carmel Mountains, in the north, to the Sinai Desert in the south, and from the Judea and Samaria Mountains in the east, to the Mediterranean coastline in the west.This is a fresh water aquifer, and it is mainly used for drinking water.The aquifer naturally discharges into two springs: Yarkon and Taninim.However, today discharge occurs mainly by wells.Due to economic reasons, there are no production wells in the study area.The main production is from the confined area of the aquifer-at the Judea and Samaria Mountains foothills.The local geology of the study area is composed of the carbonated Judea group in which epikarst features are common, and contain a weathered zone of enhanced porosity near the surface or the soil/bedrock contact [28,29].The groundwater level is deep, ranging from a few tens of meters to two hundred meters below ground level [27,30,31].The general flow direction is westward.The flow characteristics in the unsaturated zone were described by [32,33], who demonstrated a source-response mechanism between precipitation and groundwater level.The vegetation in study area is defined as Maquis and forests [34].The dominant plant cover is evergreen shrubs, with the main species being Palestine oak and Pistacia palaestina.
The current research examines the annual groundwater recharge which occurred during the hydrological year of 2014/2015 in three sub-basins (EK4, Sataf, and Matash) downstream of the Beit Zait reservoir (Table 1).During the study year, four runoff and infiltration events were monitored.The first event took place between 26 November and 28 November 2014 in the sub-basins below the Beit Zait Dam as a result of intense precipitation.The subsequent two events (the second and third events) occurred due to the dam overflow between 11 January and 13 January 2015, and 19 February and 23 February 2015, respectively.The last event took place between 6 June and 14 June 2015 as a result of a controlled discharge of the dam/reservoir water.In total, during this hydrological year, 0.99 Mm 3 (million m 3 ), 1.11 Mm 3 , and 0.43 Mm 3 of runoff water flowed into EK4, Sataf, and Matash sub-basins, respectively.The study area is located within the recharge area of the fractured carbonate Western Mountain aquifer (Yarkon-Taninim) [27].The aquifer extends from south of the Carmel Mountains, in the north, to the Sinai Desert in the south, and from the Judea and Samaria Mountains in the east, to the Mediterranean coastline in the west.This is a fresh water aquifer, and it is mainly used for drinking water.The aquifer naturally discharges into two springs: Yarkon and Taninim.However, today discharge occurs mainly by wells.Due to economic reasons, there are no production wells in the study area.The main production is from the confined area of the aquifer-at the Judea and Samaria Mountains foothills.The local geology of the study area is composed of the carbonated Judea group in which epikarst features are common, and contain a weathered zone of enhanced porosity near the surface or the soil/bedrock contact [28,29].The groundwater level is deep, ranging from a few tens of meters to two hundred meters below ground level [27,30,31].The general flow direction is westward.The flow characteristics in the unsaturated zone were described by [32,33], who demonstrated a source-response mechanism between precipitation and groundwater level.The vegetation in study area is defined as Maquis and forests [34].The dominant plant cover is evergreen shrubs, with the main species being Palestine oak and Pistacia palaestina.
The current research examines the annual groundwater recharge which occurred during the hydrological year of 2014/2015 in three sub-basins (EK4, Sataf, and Matash) downstream of the Beit Zait reservoir (Table 1).During the study year, four runoff and infiltration events were monitored.The first event took place between 26 November and 28 November 2014 in the sub-basins below the Beit Zait Dam as a result of intense precipitation.The subsequent two events (the second and third events) occurred due to the dam overflow between 11 January and 13 January 2015, and 19 February and 23 February 2015, respectively.The last event took place between 6 June and 14 June 2015 as a result of a controlled discharge of the dam/reservoir water.In total, during this hydrological year, 0.99 Mm 3 (million m 3 ), 1.11 Mm 3 , and 0.43 Mm 3 of runoff water flowed into EK4, Sataf, and Matash sub-basins, respectively.

Monitoring
Hydrological monitoring was performed at four gauging stations in order to monitor runoff infiltration of surface flow.Data loggers (uncertainty of ±0.05%) were used to measure the water level in each gauging station in short time steps (4 min).The water level data was then translated into discharge using calculated and calibrated rating curves from the Israel Hydrological Service.The upper-stream gauging station, Runoff Observation (RO) Dam, is located 400 m downstream from the reservoir.The RO Dam station monitors surface flow from dam overflow, controlled drainage from the reservoir, and runoff from the adjacent sub-basin, located between the dam and the station (3.3 km 2 ).The other three monitoring stations, RO EK4, RO Sataf, and RO Matash, are located below each of the studied sub-basins, with catchment areas of 6.4 km 2 , 5.0 km 2 , and 6.3 km 2 , respectively.The inflow toward the above-mentioned stations is composed of both the water that flows from the reach above, and the drainage of the respective catchment area.Precipitation rates, which were measured in 10 min intervals, were obtained from the Israel Meteorological Service (IMS) Tzuba and Jerusalem Stations (Figure 1).Daily calculated evapotranspiration rates were obtained from Tzuba station, based on the Penman-Monteith method [36].

Flow Model
A lumped continuous model was built using the Hydrologic Engineering Center-Hydrologic Modelling System (HEC-HMS) in order to simulate surface flow and infiltration toward the groundwater.The HEC-HMS is designed to model the precipitation-runoff processes of dendritic watershed systems [12].The program includes a variety of separate mathematical sub-models to represent each component of the runoff process, including models that compute rainfall losses, runoff generation, base flow, and channel routing [37].For the present investigation, we distinguished between diffuse groundwater recharge that occurs on the hillslopes over the majority of the study area and direct groundwater recharge which results from surface flow infiltration through highly developed karst conduits that are distributed along the Soreq creek.Each of the distinct groundwater recharge processes were quantified with a specific set of HEC-HMS sub-models.

Hillslopes "Diffuse" Groundwater Recharge
Diffuse groundwater recharge and excess precipitation were calculated for each sub-basin by a "deficit and constant" loss function and by a "simple canopy model" [12].A kinematic wave function was used to transform the excess precipitation into discharge volumes.The measured discharge volumes were used to calibrate the model.
The "deficit and constant" loss function uses a single soil layer to account for the continuous changes in soil moisture due to excess precipitation (pe t (L)), and allows for a reduction in soil moisture during periods without rain due to evapotranspiration [12]: where P t (LT −1 ) is the precipitation intensity during one-time step (t), P i (L) is the cumulative precipitation depth, and I c (L) is the canopy interception.I a (L), the soil storage deficit, represents the maximum precipitation depth that can fall without runoff, and is used as an initial condition to represent interception, depression storage, and soil moisture deficit (I a and I c are emptied between events).f c (LT −1 ) is the maximum rate of percolation from the soil layer to deeper layers, which occurs once the soil moisture reaches field capacity.
R Di f f (LT −1 ) is the diffuse groundwater recharge rate that percolates into deeper layers once the soil moisture reaches field capacity: Canopy interception and evapotranspiration were calculated with the simple canopy method [12], in which all precipitation is intercepted until the canopy storage capacity is filled.Once filled, all further precipitation is stored in the soil storage using the deficit and constant functions.The canopy storage is emptied in periods without rain, according to the potential evapotranspiration values.Once the canopy storage is exhausted, all unused potential evapotranspiration empties the water that is stored in the soil layer.
The kinematic wave flow routing function was used for transforming the excess precipitation from the hillslopes into runoff.This function is based on the Saint-Venant equation, which is derived from the conservation of mass and conservation of momentum equations.Using the assumption that the energy slope is parallel to the channel slope and the flow is steady and uniform, the complex Saint-Venant equation can be reduced to the kinematic wave approximation [38]: where h is the depth of ponding (L), q is the water discharge per unit width (L 2 T −1 ), R is the precipitation rate (LT −1 ), I is the infiltration rate (LT −1 ), t is time (T), and x is the distance along the slope (L).The depth-discharge relationship is expressed as: where m = 5/3, the parameter α is computed according to Manning's hydraulic resistance law [12] as: where S is the surface slope (L −1 ), and n is the Manning's roughness coefficient for overland flow.

"Direct" Groundwater Recharge along the Creek
The direct groundwater recharge rate R Dir (L 3 T −1 ) results from the infiltration rate in the riverbed through sinkholes and karst systems.It was calculated using a constant loss function with two parameters, constant flow rate loss and fraction loss.This loss function reduces the discharge, Q cl (L 3 T −1 ), from the combined outflow of each sub-basin, and then multiplies the remaining discharge by the fraction loss parameter, f , which results in the sub-basin outlet discharge, Q out (L 3 T −1 ).The function can be specified using the following equation: where ) is inflow that is added from upstream, and Q hs (L 3 T −1 ) is runoff developed at the sub-basin hillslopes.The direct recharge is the difference between reaches inflow and outflow at each sub-basin: Runoff along the Soreq creek was simulated with the kinematic wave transform function and the transmission loss sub-model.Runoff water within the Soreq creek was routed using the same transform function that was used for the hillslopes, but with different geometrical attributes that represent the Soreq channel.

Model Setup
The flow model setup was based on both basin and meteorological models.The basin model included the following hydrologic elements: sub-basin, reach, junction, reservoir, and the mathematical hydrologic sub-models that corresponded to each element, as was described earlier.The basin and sub-basin boundaries, as well as creek networks, were delineated using digital elevation model (DEM) data (Figure 1), which was obtained from a 5 m interval contour map.The geometry of each sub-basin is characterized by two rectangular planes with equivalent lengths that spread along both sides of the reach.The plane's equivalent length and area were derived from the 5 m interval contour map.Main channel geometry was determined based on field observations.Meteorological data were obtained from the Tzuba and Jerusalem meteorological stations.Precipitation rates were calculated for each sub-basin by the Thiessen polygon method.Evaporation rates were taken from the Tzuba station for all sub-basins.

Model Calibration and Sensitivity Analysis
Simulations were performed for the 2014/15 hydrological year with a time step interval of 10 min.Observations of runoff discharge from 2012/13 and 2013/14 were used to calibrate the model using the trial and error method to fit the simulated water level in the channel to the observed level in each of the four gauging stations.The following parameters were estimated: precipitation loss ( f c ), Manning's roughness coefficient (n), constant flow rate loss (Q cl ), and flow fraction loss ( f ).The maximum soil storage parameter (I a ) was based on [35] who calibrated those values for terra rossa soil in the Judea mountain outcrops in an area which included the present study location.Maximum canopy interception (I c ) was determined from the National Land Cover Database (NLCD).The rest of the parameters, which are related to field geometry, were set based on field observations and measurements.
Sensitivity analysis was carried out to assess the effect of the hydraulic parameters and soil properties on simulated outlet discharge.The analysis included varying a single parameter from its best fit value, which was determined from the calibrated model, and comparing the relative influence of the change on the simulated runoff peak, the total flow volume and the Nash-Sutcliffe coefficient [39].

Runoff Volumes and Flow Dynamics
The runoff discharge varied between each of the flow events (Figure 2).The four flow events are categorized by three differing levels of intensity.The first and the fourth events were considered as moderate and low intensity events, respectively.The second and the third events were considered high intensity events.The high intensity events occurred during the second part of the rainy season, when the surface soil layer was already at field capacity and the water level at the Beit Zait Reservoir was high (about 1 m below the overflow level).These events occurred during continuous rainstorms, with maximum daily rainfall depth of 49.2 mm in January (second event) and 50.6 mm in February (third event) causing dam overflow.
influence of the change on the simulated runoff peak, the total flow volume and the Nash-Sutcliffe coefficient [39].

Runoff Volumes and Flow Dynamics
The runoff discharge varied between each of the flow events (Figure 2).The four flow events are categorized by three differing levels of intensity.The first and the fourth events were considered as moderate and low intensity events, respectively.The second and the third events were considered high intensity events.The high intensity events occurred during the second part of the rainy season, when the surface soil layer was already at field capacity and the water level at the Beit Zait Reservoir was high (about 1 m below the overflow level).These events occurred during continuous rainstorms, with maximum daily rainfall depth of 49.2 mm in January (second event) and 50.6 mm in February (third event) causing dam overflow.The moderate intensity event, which occurred in November, was caused by a short but intense storm, with a daily rainfall depth of 131.6 mm.Despite the relatively high precipitation rates during this event, the dam did not overflow because the reservoir was nearly empty (more than 10 m below the overflow level) prior to the event.The last (fourth) runoff event was an artificial event, meaning The moderate intensity event, which occurred in November, was caused by a short but intense storm, with a daily rainfall depth of 131.6 mm.Despite the relatively high precipitation rates during this event, the dam did not overflow because the reservoir was nearly empty (more than 10 m below the overflow level) prior to the event.The last (fourth) runoff event was an artificial event, meaning it was not affected by precipitation, and occurred when the surface soil layer was below field capacity.
Water 2018, 10, 79 8 of 16 Total runoff volumes along the creek also varied during each event between the RO stations (Table 2).In both high intensity (second and third) events, similar patterns of total flow volume fluctuation were observed along the creek (Figure 3).In both events, the total flow increased between the RO Dam and RO EK4 stations, increasing 7.2% in the second event and 27.1% in the third event.
In the second event, a reduction in total flow volume was observed further downstream, at the Sataf and the Matash stations.Compared to the previous station, a 64.0%reduction in total flow was measured at the Sataf station, and a further reduction of 51.4% at the Matash station.During the third event, there was also an observed decrease in total flow volume downstream, where the total flow decreased 57.3% and 40.0%, from EK4 to the Sataf, and from the Sataf to the Matash stations, respectively.In both high intensity (second and third) events, the flow at the RO Dam gauging station was higher than the flow at the RO Sataf station, which implies high reach losses.it was not affected by precipitation, and occurred when the surface soil layer was below field capacity.
Total runoff volumes along the creek also varied during each event between the RO stations (Table 2).In both high intensity (second and third) events, similar patterns of total flow volume fluctuation were observed along the creek (Figure 3).In both events, the total flow increased between the RO Dam and RO EK4 stations, increasing 7.2% in the second event and 27.1% in the third event.In the second event, a reduction in total flow volume was observed further downstream, at the Sataf and the Matash stations.Compared to the previous station, a 64.0%reduction in total flow was measured at the Sataf station, and a further reduction of 51.4% at the Matash station.During the third event, there was also an observed decrease in total flow volume downstream, where the total flow decreased 57.3% and 40.0%, from EK4 to the Sataf, and from the Sataf to the Matash stations, respectively.In both high intensity (second and third) events, the flow at the RO Dam gauging station was higher than the flow at the RO Sataf station, which implies high reach losses.During the first (moderate intensity) event, the total flow increased 548.2% between the RO Dam and the RO EK4 stations.Downstream, similar to the pattern identified in the high intensity events, a 61.9% reduction in overall flow was documented at the Sataf sub-basin outlet, compared to the EK4 outlet; an additional reduction of 39.6% was documented at the Matash sub-basin outlet During the first (moderate intensity) event, the total flow increased 548.2% between the RO Dam and the RO EK4 stations.Downstream, similar to the pattern identified in the high intensity events, a 61.9% reduction in overall flow was documented at the Sataf sub-basin outlet, compared to the EK4 outlet; an additional reduction of 39.6% was documented at the Matash sub-basin outlet compared to the Sataf sub-basin outlet.However, unlike the high intensity flow events, the total flow volume was Water 2018, 10, 79 9 of 16 lower at the RO Dam station than at the RO Sataf station, implying a relatively high contribution of runoff water from the hillslopes into the creek.
The fourth (low intensity/artificial) event reveals a different flow pattern in which total flow volume reduced by 82.4% between the RO Dam and the EK4 gauging stations.There was no flow measured at the two lower gauging stations during this event.

Model Calibration
Calibrating the flow model to the observed discharge in the RO Dam, RO EK4, RO Sataf, and RO Matash gauging stations yielded good agreement between measured and observed discharge rates.The Nash-Sutcliffe coefficient for the runoff discharge was 95.2%.Calibrated model parameters are presented in Table 1.
The model simulates diffuse recharge from the hillslopes, and direct recharge along the creek through fast flow paths.The model predicts, quite well, the discharge rates of the first and moderate intensity events at all of the monitoring sites (Figure 4).On the other hand, relatively high differences were detected in peak discharges at the RO Dam and RO Matash gauging stations.The model described the flood discharge during the high intensity events fairly well.Peak discharge was underestimated in all events, with the exception of the RO Matash gauging station during the January (second) event.Generally, the calculated recession limb was longer than that measured in all sites, except at the Sataf site during the second event.Additionally, during the late period of the January (second) event, a relatively small flow episode developed at the EK4 sub-basin.This flow episode was poorly predicted for both the RO Dam and RO EK4 gauging stations dataset.The total discharge at RO Dam station during the artificial low intensity event was predicted very well.However, the flow was overestimated at the RO EK4 site.
Water 2018, 10, 79 9 of 16 compared to the Sataf sub-basin outlet.However, unlike the high intensity flow events, the total flow volume was lower at the RO Dam station than at the RO Sataf station, implying a relatively high contribution of runoff water from the hillslopes into the creek.The fourth (low intensity/artificial) event reveals a different flow pattern in which total flow volume reduced by 82.4% between the RO Dam and the EK4 gauging stations.There was no flow measured at the two lower gauging stations during this event.

Model Calibration
Calibrating the flow model to the observed discharge in the RO Dam, RO EK4, RO Sataf, and RO Matash gauging stations yielded good agreement between measured and observed discharge rates.The Nash-Sutcliffe coefficient for the runoff discharge was 95.2%.Calibrated model parameters are presented in Table 1.
The model simulates diffuse recharge from the hillslopes, and direct recharge along the creek through fast flow paths.The model predicts, quite well, the discharge rates of the first and moderate intensity events at all of the monitoring sites (Figure 4).On the other hand, relatively high differences were detected in peak discharges at the RO Dam and the RO Matash gauging stations.The model described the flood discharge during the high intensity events fairly well.Peak discharge was underestimated in all events, with the exception of the RO Matash gauging station during the January (second) event.Generally, the calculated recession limb was longer than that measured in all sites, except at the Sataf site during the second event.Additionally, during the late period of the January (second) event, a relatively small flow episode developed at the EK4 sub-basin.This flow episode was poorly predicted for both the RO Dam and RO EK4 gauging stations dataset.The total discharge at RO Dam station during the artificial low intensity event was predicted very well.However, the flow was overestimated at the RO EK4 site.

Cumulative Groundwater Recharge
Cumulative groundwater recharge from both diffuse and direct infiltration varied between each sub-basin (Figure 5).The highest total recharge depth (both diffuse and direct infiltration) was calculated for the Sataf sub-basin at 305.7 mm.In this sub-basin, the total recharge was the result of all four flow events (Table 3).As a result of the first, moderate intensity, flow event the total groundwater recharge was calculated to be 89.9 mm; 4.9% of this amount was contributed to by direct infiltration.The total recharge at the end of second and the third (both high intensity) events, was calculated to be 130.4 mm and 75.3 mm, respectively.In these events, the direct infiltration accounted for 10.4% of the total infiltration during the first event, and 40.7% during the second event.Direct infiltration was greater in the high intensity events compared to the moderate intensity event.During the fourth low intensity and artificial event, only 10.1 mm percolated into the aquifer by direct infiltration; due to the absence of precipitation, there was no diffuse infiltration from the hillslopes.

Cumulative Groundwater Recharge
Cumulative groundwater recharge from both diffuse and direct infiltration varied between each sub-basin (Figure 5).The highest total recharge depth (both diffuse and direct infiltration) was calculated for the Sataf sub-basin at 305.7 mm.In this sub-basin, the total recharge was the result of all four flow events (Table 3).As a result of the first, moderate intensity, flow event the total groundwater recharge was calculated to be 89.9 mm; 4.9% of this amount was contributed to by direct infiltration.The total recharge at the end of second and the third (both high intensity) events, was calculated to be 130.4 mm and 75.3 mm, respectively.In these events, the direct infiltration accounted for 10.4% of the total infiltration during the first event, and 40.7% during the second event.Direct infiltration was greater in the high intensity events compared to the moderate intensity event.During the fourth low intensity and artificial event, only 10.1 mm percolated into the aquifer by direct infiltration; due to the absence of precipitation, there was no diffuse infiltration from the hillslopes.An intermediate level of the total recharge depth of 253.1 mm was calculated for the Matash sub-basin.Infiltration into this sub-basin resulted from the first three events only.The first (moderate intensity) flow event contributed 82.7 mm to the aquifer recharge.In this event, 2.5% of the total depth percolated into the aquifer by direct infiltration.The total infiltration depths from the two high intensity flow events were 118.7 mm and 51.7 mm, respectively.Similar to the Sataf sub-basin, the portion of the direct infiltration rose during the high intensity events, to 4.0% and 20.1%, respectively.A relatively low total recharge depth of 223.4 mm was calculated for the EK4 sub-basin.
Groundwater recharge was calculated for all of the four events.The first event added a recharge depth of 68.3 mm, where 1.1% of which was contributed by direct infiltration.As in the other sub-basins, the total recharge depths from the two high intensity events were higher, and reached 107.7 mm and 39.7 mm, respectively.Direct infiltration accounted for 3.1% and 18.3% of said recharge.A recharge depth of 7.7 mm was calculated for the last flow event, during which the only contribution was from direct infiltration.
Overall, direct infiltration during floods contributed the most to groundwater recharge during the second high intensity (and third overall) event.During the controlled runoff (fourth) event, the water percolated into the subsurface by direct infiltration only, in the upper sub-basins, before it reached the lower Matash sub-basin.Direct infiltration accounted for 6.8% of the total annual recharge in the Matash sub-basin, 8.5% in the EK4 sub-basin, and 19.2% in the Sataf sub-basin.

Sensitivity Analysis
In order to examine the sensitivity of the flow model to parameter variations, parameters were varied within ±20% of their calibrated values.A sensitivity analysis was performed for the Manning's roughness coefficient (n), the flow fraction loss coefficient ( f ), the maximum soil storage parameter (I a ), and the precipitation loss parameter ( f c ).The results of the sensitivity analysis show that the model is most sensitive to changes in the flow fraction loss coefficient, and least sensitive to changes in the Manning's roughness coefficient, the maximum soil storage and the precipitation loss parameters (Figure 6).
Variation of the f parameter resulted in the largest changes in both the flood peak and the total volume of the runoff.The deviation in this parameter is related to the direct infiltration, and therefore, it is related to the contribution of direct infiltration to the annual recharge volume.By reducing the f parameter by 20%, the simulated runoff peak could increase up to 87%, and the total event flow volume by 89%.Practically, this means that less water would reach the groundwater.On the other hand, the impact of increasing the f parameter has a smaller impact; a 20% increase in the f parameter results in a reduction of both the simulated runoff peak and the total event flow volume by 21%.
Therefore, the effect of the f parameter variation on the surface runoff configuration illustrates the important influence of the karst conduits and their occurrence along the creek on the potential water capacity and infiltration into the subsurface during flow events.
Water 2018, 10, 79 12 of 16 total event flow volume by 21%.Therefore, the effect of the f parameter variation on the surface runoff configuration illustrates the important influence of the karst conduits and their occurrence along the creek on the potential water capacity and infiltration into the subsurface during flow events.
Figure 6.Results of the sensitivity analysis.The impact of changing each parameter on the simulation results, as expressed in percent variation with respect to the respective optimum simulated value.

Proportion of Groundwater Recharge during Runoff Events in Relation to the Annual Recharge for Karst Aquifers
The results of the runoff model simulations support the assumption that direct infiltration through karst conduits along the creek can have a considerable contribution to annual groundwater recharge.The calculated annual recharge volumes for all three sub-basins lay within the expected recharge range for the study area (Figure 7).In comparison to the other two sub-basins, EK4 had a lower annual precipitation rate, which resulted in a lower volume of annual recharge.The simulation results are in good agreement with earlier studies.These studies proposed the relationship between recharge and precipitation.In all sub-basins, the calculated recharge from diffuse infiltration fell between the expected recharge/precipitation relation values, according to [40,41].The total groundwater recharge that was calculated as a result of the direct infiltration in the EK4 and Matash sub-basins fell between the [41,42] proposed recharge/precipitation relation.However, the contribution of direct infiltration to the groundwater recharge in the Sataf sub-basin, which has a well-developed karst system along the creek, suggests higher total infiltration.In this case, the total annual recharge and the recharge/precipitation relation is closer to the suggested fit from an earlier study [43].However, it is proposed that this exception relates to the existence of local karst features along the Soreq creek, and although it does not represent the regional conditions of the aquifer, it occurs spatially over the aquifer outcrops [44][45][46].In addition, it is also important to emphasize that, due to lack of information about the spatial variability of karst system properties in Figure 6.Results of the sensitivity analysis.The impact of changing each parameter on the simulation results, as expressed in percent variation with respect to the respective optimum simulated value.

Proportion of Groundwater Recharge during Runoff Events in Relation to the Annual Recharge for Karst Aquifers
The results of the runoff model simulations support the assumption that direct infiltration through karst conduits along the creek can have a considerable contribution to annual groundwater recharge.The calculated annual recharge volumes for all three sub-basins lay within the expected recharge range for the study area (Figure 7).In comparison to the other two sub-basins, EK4 had a lower annual precipitation rate, which resulted in a lower volume of annual recharge.The simulation results are in good agreement with earlier studies.These studies proposed the relationship between recharge and precipitation.In all sub-basins, the calculated recharge from diffuse infiltration fell between the expected recharge/precipitation relation values, according to [40,41].The total groundwater recharge that was calculated as a result of the direct infiltration in the EK4 and Matash sub-basins fell between the [41,42] proposed recharge/precipitation relation.However, the contribution of direct infiltration to the groundwater recharge in the Sataf sub-basin, which has a well-developed karst system along the creek, suggests higher total infiltration.In this case, the total annual recharge and the recharge/precipitation relation is closer to the suggested fit from an earlier study [43].However, it is proposed that this exception relates to the existence of local karst features along the Soreq creek, and although it does not represent the regional conditions of the aquifer, it occurs spatially over the aquifer outcrops [44][45][46].In addition, it is also important to emphasize that, due to lack of information about the spatial variability of karst system properties in places other than the Soreq creek, the current model did not include other potential contributions to fast flow infiltration into the groundwater, such as sinkholes that are located outside the creek bed.Therefore, the actual groundwater recharge volume is likely to be higher than predicted by the model.places other than the Soreq creek, the current model did not include other potential contributions to fast flow infiltration into the groundwater, such as sinkholes that are located outside the creek bed.Therefore, the actual groundwater recharge volume is likely to be higher than predicted by the model.
Figure 7. Computed annual recharge of total and diffuse infiltration at the EK4, Sataf, and Matash sub-basins over the expected recharge range (grey) at the study area [32].The recharge range is between fast flow paths along conduits and slow flow paths through the matrix.EK11 and EK12 are sites that represent these two different types of flow paths, respectively.

Implications for Groundwater Enrichment Potential
Groundwater recharge from both diffuse and the direct infiltration varied in space, depending on the sub-basins, and in time for different flow events (Table 3).Assessing the differences in potential groundwater recharge in each basin may provide important information for improving watershed management practices and optimizing groundwater recharge.
During the 2014/15 hydrological year, more than 43% of the total annual recharge in all sub-basins occurred as a result of the first high intensity (second overall) flood event.However, only up to 4.5% of the total recharge volume from that event was from direct infiltration through karst systems along the Soreq creek.The moderate intensity (first overall) flow event contributed about 30% of the total annual recharge in all the sub-basins.In this event, the contribution of direct infiltration to the annual water balance was lower than 1.5%.The second high intensity (third overall) flow event contributed 18%, 25%, and 20% to the annual recharge in the EK4, Sataf, and Matash sub-basins, respectively.In this event, direct infiltration contributed more to the annual water balance than the two previous events, reaching up to 3.2% at EK4 sub-basin to 10.0% at the Sataf sub-basin.The relatively high contribution of direct infiltration is the result of high discharge rates from along the creek and the long duration of the flood event in comparison with the previous events.In comparison, all the recharge during the artificial flow event was from direct infiltration.Remarkably, although this flow event only occurred in two sub-basins, it contributed more than 3.3% to the annual recharge volume, which is about 40% and 20% of the annual direct infiltration in the EK4 and the Sataf sub-basins, respectively.In total, during the artificial event, 0.10 Mm 3 percolated into the sub-surface over the course of nine days.The total flow loss at the Matash outlet during the two high intensity events was 0.22 Mm 3 , which accounted for 88% of the total annual water loss in this outlet.The more in-depth understanding of the water balance dynamics in the basin provided by the model can also help improve watershed management practices.During the winter of 2014/15, the interval between major precipitations events exceeded 20 days, therefore, by Computed annual recharge of total and diffuse infiltration at the EK4, Sataf, and Matash sub-basins over the expected recharge range (grey) at the study area [32].The recharge range is between fast flow paths along conduits and slow flow paths through the matrix.EK11 and EK12 are sites that represent these two different types of flow paths, respectively.

Implications for Groundwater Enrichment Potential
Groundwater recharge from both diffuse and the direct infiltration varied in space, depending on the sub-basins, and in time for different flow events (Table 3).Assessing the differences in potential groundwater recharge in each basin may provide important information for improving watershed management practices and optimizing groundwater recharge.
During the 2014/15 hydrological year, more than 43% of the total annual recharge in all sub-basins occurred as a result of the first high intensity (second overall) flood event.However, only up to 4.5% of the total recharge volume from that event was from direct infiltration through karst systems along the Soreq creek.The moderate intensity (first overall) flow event contributed about 30% of the total annual recharge in all the sub-basins.In this event, the contribution of direct infiltration to the annual water balance was lower than 1.5%.The second high intensity (third overall) flow event contributed 18%, 25%, and 20% to the annual recharge in the EK4, Sataf, and Matash sub-basins, respectively.In this event, direct infiltration contributed more to the annual water balance than the two previous events, reaching up to 3.2% at EK4 sub-basin to 10.0% at the Sataf sub-basin.The relatively high contribution of direct infiltration is the result of high discharge rates from along the creek and the long duration of the flood event in comparison with the previous events.In comparison, all the recharge during the artificial flow event was from direct infiltration.Remarkably, although this flow event only occurred in two sub-basins, it contributed more than 3.3% to the annual recharge volume, which is about 40% and 20% of the annual direct infiltration in the EK4 and the Sataf sub-basins, respectively.In total, during the artificial event, 0.10 Mm 3 percolated into the sub-surface over the course of nine days.The total flow loss at the Matash outlet during the two high intensity events was 0.22 Mm 3 , which accounted for 88% of the total annual water loss in this outlet.The more in-depth understanding of the water balance dynamics in the basin provided by the model can also help improve watershed management practices.During the winter of 2014/15, the interval between major precipitations events exceeded Water 2018, 10, 79 14 of 16 20 days, therefore, by performing a controlled drainage of the Beit Zait reservoir and diverting runoff waters into the Soreq sub-basin, the water could have been used to supplement the aquifer.

Conclusions
A lumped continuous model (HEC-HMS) was applied to simulate runoff flow and karst aquifer recharge in three sub-basins along the ephemeral Soreq stream.Two primary infiltration mechanisms were identified during the study.The first is diffuse groundwater recharge, which occurs at hillslopes over the majority of study area.The second infiltration mechanism is direct groundwater recharge, which happens as a result of surface flow infiltration through highly developed karst conduits along the Soreq creek.
Results of simulations show that such an approach can be used to describe the spatial distribution of groundwater recharge.The calculated recharge volumes align with results from previous studies, however, the present methodology also allowed for the quantification of the annual recharge via direct infiltration through the creek bed.In this highly developed karst system, direct infiltration accounted for 19.2% of the annual recharge volume.Moreover, the present methodology allows for the assessment of potential infiltration from the creek bed, and subsequently, the user is able to evaluate watershed management practices that could potentially improve artificial enrichment of the aquifer.
The main limitation of the applied model is the lack of information about the spatial variability of karst system structure.The model did not include other potential contributors to fast flow infiltration into the groundwater, such as sinkholes, other than from the Soreq creek bed.Therefore, actual groundwater recharge volume is expected to be higher.Another limitation is the assumption that water percolated straight into the groundwater.In places where there is a thick unsaturated zone, infiltration is not necessarily immediate, and could potentially take several hydrological years to create a hydraulic response at the water table.

Figure 3 .
Figure 3.Total flow volume during the four runoff events at each monitoring site along the Soreq creek.

Figure 3 .
Figure 3.Total flow volume during the four runoff events at each monitoring site along the Soreq creek.

Figure 7 .
Figure 7.Computed annual recharge of total and diffuse infiltration at the EK4, Sataf, and Matash sub-basins over the expected recharge range (grey) at the study area[32].The recharge range is between fast flow paths along conduits and slow flow paths through the matrix.EK11 and EK12 are sites that represent these two different types of flow paths, respectively.

Table 2 .
Total measured runoff (m 3 ) at the monitoring stations in the study site during all studied flow events.

Table 2 .
Total measured runoff (m 3 ) at the monitoring stations in the study site during all studied flow events.

Table 3 .
Computed groundwater recharge depth (mm), and in parentheses, the percentages from the total recharge in each sub-basin during the study period.