Recharge Estimation and Soil Moisture Dynamics in a Mediterranean, Semi-arid Karst Region

Knowledge of soil moisture dynamics in the unsat-urated soil zone provides valuable information on the temporal and spatial variability of groundwater recharge. This is especially true for the Mediterranean region, where a substantial fraction of long-term groundwater recharge is expected to occur during high magnitude precipitation events of above-average wet winters. To elucidate process understanding of infiltration processes during these extreme events, a monitoring network of precipitation gauges, meteorological stations, and soil moisture plots was installed in an area with a steep climatic gradient in the Jordan Valley region. In three soil moisture plots, Hydrus-1D was used to simulate water movement in the unsaturated soil zone with soil hydraulic parameters estimated by the Shuffled Complex Evolution Metropolis algorithm. To generalize our results, we modified soil depth and rainfall input to simulate the effect of the pronounced climatic gradient and soil depth variability on percolation fluxes and applied the calibrated model to a time series with 62 years of meteorological data. Soil moisture measurements showed a pronounced sea-sonality and suggested rapid infiltration during heavy rainstorms. Hydrus-1D successfully simulated short and long-term soil moisture patterns, with the majority of simulated deep percolation occurring during a few intensive rainfall events. Temperature drops in a nearby groundwater well were observed synchronously with simulated percolation pulses, indicating rapid groundwater recharge mechanisms. The 62-year model run yielded annual percolation fluxes of up to 66 % of precipitation depths during wet years and of 0 % during dry years. Furthermore, a dependence of recharge on the temporal rainfall distribution could be shown. Strong correlations between depth of recharge and soil depth were also observed.


Introduction
In the Mediterranean region, groundwater is the main source for domestic and agricultural water supplies (EUWI, 2007).Knowledge on the quantity of groundwater recharge is a prerequisite for sustainable water resources planning and effective water use.Small-scale differences in climate, geology, land use, topography and soil properties cause a high spatial and temporal variability of groundwater recharge making the assessment and predictions of recharge a challenge (e.g.Zagana et al., 2007).Karst areas are important in this respect, because during high intensity winter storms precipitation may rapidly infiltrate into exposed karst surfaces and induce high recharge rates (De Vries and Simmers, 2002), which are common in the Mediterranean area (Ford and Williams, 2007).A rapidly increasing water demand in the last decades has led to a widespread overexploitation of groundwater resources (EUWI, 2007).Furthermore, the Mediterranean region has been identified as a "hotspot" of current and future climate change (Giorgi, 2006;IPCC, 2013), imposing additional pressure on its limited water resources.Hence, more insights into processes of aquifer replenishment in Mediterranean karst regions are of vital importance.
A large variety of methods suitable for estimating recharge rates were developed in the last decades (De Vries and Simmers, 2002;Scanlon et al., 2002).Infiltration, percolation and recharge quantities in Mediterranean karst have mainly F. Ries et al.: Recharge estimation and soil moisture dynamics been approached from two sides: on the one hand, hydrologists and geomorphologists characterized the surface water balance on small plots by sprinkling experiments or by runoff measurements during natural rainstorms (e.g.Cerdà, 1998;Lavee et al., 1998).Large-scale experiments also included tracers and facilitated statements on runoff generation processes (e.g.Lange et al., 2003).However, these studies quantified infiltration by the difference between artificial/natural rainfall and measured overland flow but did not differentiate between recharge and evapotranspiration.On the other hand, hydrogeologists frequently assessed average recharge rates of entire karst catchments from spring discharge measurements or hydraulic head data.Methods include knowledge GIS-based mapping (Andreo et al., 2008), multiple linear regression (Allocca et al., 2014), conceptual models (e.g.Hartmann et al., 2013a), coupled water-balance groundwater models (Sheffer et al., 2010), and chloride mass balances (Marei et al., 2010;Schmidt et al., 2013).However, these studies treat karst systems as units, including both the unsaturated and the saturated zones, and are limited in temporal and spatial resolution.Studies on cave drips (Gregory et al., 2009;Arbel et al., 2010;Lange et al., 2010) provided insights into the deeper unsaturated zone in terms of water storage, spatial variability of percolation and flow paths.Their data was also used to incorporate variability in recharge modelling (Hartmann et al., 2012).However, it was difficult to distinguish between processes in the unsaturated soil zone and in the underlying epikarst, and uncertainty remains regarding the representativeness of cave drip data with respect to infiltration processes.This is mainly due to the fact that the contributing areas of cave drips are unknown and caves might have developed their own hydraulic environments.Therefore, cave drips are not necessarily representative for the bulk karst vadose zone (Lange et al., 2010).
Only limited knowledge on recharge dynamics is available for the carbonate Mountain Aquifer system shared between the West Bank and Israel, although it is of regional importance.First recharge estimates were based on longterm spring discharge and groundwater well abstraction data (Goldschmidt and Jacobs, 1958).Later, groundwater flow models were used to establish empirical rainfall-recharge relationships (Baida and Burstein, 1970;Guttman and Zukerman, 1995;Zukerman, 1999).Average recharge rates were assessed by a simple water balance approach (Hughes et al., 2008) and by a chloride mass balance (Marei et al., 2010).Sheffer et al. (2010) coupled a water budget model with a groundwater flow model for the entire western part of the Mountain Aquifer and used spring discharge and groundwater level data for calibration.They reported recharge rates ranging between 9 and 40 % of annual rainfall and showed that the temporal distribution of rainfall within the winter season had considerable effects on overall recharge rates.
Observations of soil moisture may offer unique insights into near-surface hydrological processes, because water fluxes are susceptible to conditions and properties of the va-dose soil zone across several scales (Vereecken et al., 2008).Yet, soil moisture is rarely measured in semi-arid areas and is seldom used for recharge estimation purposes.Scott et al. (2000) exemplified the potential of soil moisture time series to calibrate Hydrus-1D soil hydraulic parameters in southeastern Arizona.Their results demonstrated the high inter-annual variability of water fluxes in these environments where considerable percolation only occurs during aboveaverage wet years.
The objective of this study is to investigate the spatial and temporal variability of soil water percolation and hence groundwater recharge rates for an eastern Mediterranean carbonate aquifer.We use continuously recorded soil moisture data to calibrate one-dimensional water flow models (Hydrus-1D) with the Shuffled Complex Evolution Metropolis (SCEM) algorithm.The calibrated models are then used to assess spatial and temporal patterns of soil water percolation in a Mediterranean karst area, which is characterized by strong climatic gradients and variable soil depths.
A common challenge of hydrological research in semiarid and developing regions is the lack of data.At the same time, sound knowledge on the often-limited water resources is of vital importance, especially in karst areas.This situation necessitates compromises.The calibrated soil hydraulic parameters of our model should be treated as effective parameters that represent both preferential and matrix flow components within a single, unimodal pore size distribution.They are site-specific and should not be used to characterize the physics of a porous medium with the given grain size distribution.Despite increasing work on (preferential) water transport in heterogeneous porous media, there is still no convincing integrated physical theory about non-Darcian flow at the scale of interest (Beven and German, 2013).And even if such a theory existed, measurement problems in natural clay soils would restrict its application to laboratory monoliths.From this perspective, the use of a simple model with a minimum number of calibrated parameters seemed to be a valid compromise to infer statements on groundwater recharge from a limited number of measurements in the unsaturated zone.

Study area
Our study area is located on the western margin of the Jordan Rift Valley 25 km northeast of Jerusalem (Fig. 1).Precipitation shows a pronounced seasonality with cold fronts (mainly Cyprus lows) carrying moisture from the Mediterranean Sea during winter season from October to April (Goldreich, 2003).High rainfall intensities can occur mainly in autumn and spring from convective rainfall events originating from the south (Red Sea Troughs).The topographic gradient from the mountain range (highest elevation: 1016 m a.s.l.) in the west to the Jordan Valley in the east results in a strong precipitation gradient and arid conditions in the Jordan Valley (rain-shadow desert).Long-term average annual precipita- tion decreases from 532 mm in Jerusalem (810 m a.s.l.) to 156 mm in Jericho (290 m b.s.l.) (Morin et al., 2009).Mean annual potential evapotranspiration adds up to 1350 mm in the mountains and 1650 mm in the Jordan Valley (Israel Meteorological Service -http://www.ims.gov.il).
Outcropping geological formations consist of carbonate rocks of the Upper Cretaceous Age (Begin, 1975).They are composed of fractured and highly permeable layers of limestone and dolomite alternating with marl and chalk layers of low permeability, often considered partial aquicludes (Weiss and Gvirtzman, 2007).Senonian chalks form outcrops of low hydraulic conductivity in the southeast (Rofe and Raffety, 1963).Soil parent material consists of residual clay minerals from carbonate rock weathering and from the aeolian input of dust (silt and clay fraction) originating from the Sahara (Yaalon, 1997).Predominant soil types are Terra Rossa and Rendzina, both characterized by high clay contents.Rendzina soils contain carbonate in the soil matrix, are thinner and still show recent development, whereas Terra Rossa soils were formed under past climatic conditions (Shapiro, 2006).As a result of the diverse underlying carbonate rock with different degrees of weathering and due to heterogeneous topography, soil depth is highly variable.The slopes are covered by massive bedrock exposures, and loose rock fragments of different sizes alternate with soil pockets of variable dimensions, shapes, and depths (Fig. 2).Soil development is intensified where dissolution cracks and karst fissures provide favourable drainage of the vadose soil zone to the underlying bedrock.In valley bottoms, fine-textured alluvial soils (Vertisols) with soil depths of up to several metres have de- veloped.Shallow Brown Lithosols and loessial Arid Brown Soils dominate in the eastern, low-lying areas receiving less rainfall (Shapiro, 2006).In general, soils in the region have significantly been transformed by human activities such as land cultivation, terracing, and deforestation during the last 5000 years (Yaalon, 1997).
On the hillslopes, annual plants and Mediterranean shrubs (predominantly Sarcopoterium spinosum) are the dominant vegetation types.They are used for extensive grazing by goats and sheep.South-facing slopes show lower vegetation density and higher proportion of bare soil and rock outcrops than the north-facing slopes, where the presence of biogenic crusts was reported (Kutiel et al., 1998).Minor land use types consist of scattered built-up areas, olive plantations on terraced land and rainfed or partly irrigated agricultural land (annual and perennial crops, herbs and vegetables) in valley bottoms.

Hydrometeorological measurements
To capture the spatial variation of rainfall along the strong climatic gradient, we installed a rain gauge network (Fig. 1) consisting of 14 tipping buckets (RG3-M) connected to a HOBO pendant event data logger (Onset Computer Corporation), recording 0.2 mm per tip.Daily cumulative precipitation was calculated from event data.All gauges were calibrated before employment, maintained, and cleaned twice a year before and after the rainfall season.Temperature was measured at four climatic stations (Thies GmbH and Onset Computer Corporation) at 10 min intervals.Additional rainfall and climatic data was obtained from the Israel Meteorological Service database (http://www.data.gov.il/ims) for long-term analyses.Every 20 min, groundwater levels and temperatures were recorded in a well tapping the local perched spring aquifer using pressure transducers (Mini-Diver, Eijkelkamp).Moreover, we measured water levels in several ephemeral streams of Wadi Auja with pressure transducers (Mini-Diver, Eijkelkamp; Dipper-3, SEBA Hydrome-trie) every 5 min.Irrigation experiments (Sohrt et al., 2014) demonstrated that infiltration rates at locations close to the soil moisture plots were considerably higher than measured rainfall intensities during our observation period.

Soil moisture measurements
Seven soil moisture plots were installed, each equipped with four capacitance soil moisture sensors (5TM/5TE, Decagon Devices Inc.) measuring soil moisture and soil temperature every 10 min.We paid attention that the plots did not receive lateral surplus water from upslope overland flow by placing them distant from rock outcrops and at locations with minimum slope.To minimize disturbance, we inserted the sensors vertically into the upslope wall of manually dug soil pits (depth between 50 and 100 cm).After installation, we refilled the pits with the parent soil material and compacted approximately to pre-disturbance bulk density.The probes were connected to data loggers (EM50, Decagon Devices Inc.), which were sealed by plastic bags and buried in the soil to avoid vandalism.We used the internal calibration function for mineral soils with a measurement accuracy of 4 % of the volumetric water content (VWC).Further information on the performance of the employed sensors can be found in Kizito et al. (2008).Due to instrument malfunction and vandalism, we obtained continuous data of the entire measurement period (October 2011-May 2013) from only three locations (SM-1-SM-3).Plot SM-1 is located at a gentle part of a slope, while SM-2 and SM-3 are located on rather flat topography.Further characteristics of the plots are summarized in Table 1.
The dielectric permittivity of water changes with temperature (e.g.Wraith and Or, 1999).Hence, measurement techniques of soil moisture based on the difference of dielectric permittivity between water and soil matrix are affected by this phenomenon.In our case, soil temperature was highly variable and changed by up to 20 • C within 24 h due to a strong radiation input and partly uncovered soil.We corrected our soil moisture data applying multiple linear regressions against soil temperature as described by Cobos and Campell (2007) (2010)(2011)(2012)(2013).b Rainfall at plot SM-2 is estimated by inverse distance weighted interpolation with elevation as additional predictor.c Textural characteristics were determined in the laboratory by sieving (particle size > 0.063 mm) and sedimentation methods (particle size < 0.063 mm).

Soil sampling and multistep outflow experiments
We took 35 undisturbed soil samples (height = 4 cm, diameter = 5.6 cm) with a volume of 100 cm 3 in the surrounding of the soil moisture plots in depths between 5 and 70 cm.They were analysed in the laboratory of the Forest Research Institute of Baden-Württemberg, Freiburg, Germany, by means of multistep outflow (MSO) experiments (Puhlmann et al., 2009).The setup of the MSO-experiments was based on the pressure cell method, where samples were equipped with microtensiometers, placed on porous ceramic plates and gradually saturated.Suctions of up to −500 hPa were gradually applied at the bottom of the ceramic plates.Cumulative outflow as well as the pressure head were continuously monitored and logged.Furthermore, samples were placed in a pressure plate apparatus to obtain points of the retention curves at −900 hPa.The Mualem-van Genuchten parameters were derived by means of an inverse parameter optimization procedure.We compared water retention and conductivity functions from the laboratory MSO experiments with those derived through inverse modelling of our soil moisture plots.

Modelling of the soil zone
Water balance at the plot scale in absence of surface runoff can be described by where ds/dt is the storage change over time, P is the precipitation, L is the percolation at the profile bottom and E a is the evapotranspiration per time interval.E a is composed of the terms E i (evaporation of intercepted precipitation), E s (soil evaporation) and E t (plant transpiration).
For our three soil moisture plots, soil water content and water fluxes were simulated on a daily basis with Hydrus-1D (version 4.16;Šimůnek et al., 2013) for a period of 32 months.Hydrus-1D solves the Richards equation numerically for water transport in variable saturated media.Matricpotential-dependent water retention and hydraulic conductivity were calculated using the Mualem-van Genuchten soil hydraulic model (van Genuchten, 1980).To reduce the effect of non-linearity of the hydraulic conductivity function close to saturated conditions, an air entry value of −2 cm as suggested by Vogel et al. (2001) was used.Interception by the plant canopy was calculated by an empirical equation including the leaf area index and daily precipitation values (see Šimůnek et al., 2013, for more details).Potential evapotranspiration was calculated by the Hargreaves equation (Hargreaves and Samani, 1985).Originally developed for a lysimeter station in California, this method adequately reproduced potential evapotranspiration under semi-arid climates (Jensen et al., 1997;Weiß and Menzel, 2008).Potential evapotranspiration was split into potential evaporation from the soil surface and potential transpiration from plants according to Beer's law based on the time variable surface cover fraction.Both fluxes were reduced to actual values based on a root water uptake model (Feddes et al., 1978) applying plant parameters for grass and an energy balance surface evaporation model (Camillo and Gurney, 1986).In our study area, vegetation cover shows a strong seasonality due to the restricted water availability during the dry season.To account for this, time-dependent plant growth data was implemented into the model with intra-annual variation of surface cover fraction.According to field observations, the start of the growing season was set to mid November and the maximum vegetation density was assumed for February/March shortly after the largest monthly precipitation amounts were observed.The depth from which plants took up water was controlled by a root distribution function.An exponential decrease of root density with soil depth was assumed, observed at the study sites and often reported for the Mediterranean region (e.g.De Rosnay and Polcher, 1998;De Baets  Depth of soil profile 0.5-1 m Measured at experimental plots a Parameter calibrated for each soil material with SCEM algorithm and Kling-Gupta efficiency as optimization criterion.b The upper parameter limit of r and the lower parameter limit of s were obtained from the lowest, respectively, highest measured volumetric soil moisture value of each layer in the respective soil moisture plot.c Rainfall at plot SM-2 is estimated by inverse distance weighted interpolation with elevation as additional predictor.d Maximum and minimum daily air temperature at the soil moisture plots is estimated by calculation of an elevation-temperature gradient based on meteorological stations in the Jordan Valley and the mountains.et al., 2008).Temporal variations of rooting depth and root density were disregarded.With these components, Hydrus-1D continuously computed water content and water fluxes at user-defined observation points (here: depths of the soil moisture probes) and at the lower profile boundary.Model input data, selected parameter values and their ranges, and the corresponding data sources and calculation methods are summarized in Table 2.

Calibration procedure, uncertainty analysis and parameter sensitivity
An increase of clay content and bulk density with depth was observed at all profiles and the individual probes in various depths at our plots differed noticeably.As a result, a particular soil material with singular soil hydraulic properties was independently assigned for each soil moisture probe.Ob-served soil moisture data from two winter and one summer season (October 2011-April 2013) were used for calibration of Hydrus-1D.We individually determined soil hydraulic parameters for every soil material by inverse modelling using the Shuffled Complex Evolution Metropolis optimization algorithm (SCEM; Vrugt et al., 2003) and the Kling-Gupta efficiency (KGE; Gupta et al., 2009) in a modified version from Kling et al. (2012) as the objective function: where r is the correlation coefficient between simulated and observed VWC (Cov so is the covariance between simulated and observed VWC), α is a dimensionless measure for the bias (µ s and µ o are the mean simulated and observed VWC) and β is a dimensionless measure for variability (σ s and σ s are the standard deviations of simulated and observed VWC).SCEM is widely used to efficiently solve global optimization problems (e.g.Vrugt et al., 2005;Schoups et al., 2005;Feyen, 2007;Hartmann et al., 2012) and to find optimal model parameter sets.As algorithmic parameters for SCEM, 24 complexes/parallel sequences were selected (equal to the number of parameters to be optimized), the population size was set to 144 and the number of accepted draws to infer posterior distribution was set to 1000.The SCEM routine was run until the scale reduction score (SR), a convergence criterion defined by Gelman and Rubin (1992), was fulfilled.
As proposed by Vrugt et al. (2003), a SR value of 1.2 was chosen, indicating that the Markov chain had converged to a stationary posterior distribution for all parameters.Predicted soil moisture ranges were used for parameter uncertainty assessment.They were determined by running Hydrus-1D with 1000 parameter sets obtained through the SCEM algorithm after reaching convergence.

Spatial and temporal extrapolation of percolation
To extrapolate our point measurements of soil water balance, we varied soil depth and climatic input parameters (precipitation and temperature) over ranges observed in our study area.We used the calibrated soil hydraulic parameters of our deepest (1 m) soil moisture plot (SM-1), which had sensors at 10, 25, 40 and 80 cm.Moreover, we assumed that the rooting depth was limited to the soil depth with no changes in the vertical root distribution or plant surface cover fraction.We cut off the profile according to the simulated soil depth, which reduced the number of independent soil layers when the depths fell below 60, 32.5 and 17.5 cm.For soil thicknesses exceeding 1 m, we extended the bottom layer.To simulate the range of climatic conditions with elevations between 400 and 1000 m a.s.l., we modified rainfall and air temperature according to calculated mean annual gradients based on observed rainfall and climatic data.We had three seasons of measured climate data, which we analysed separately due to seasonal differences in cumulative rainfall amount and distribution.
Using a 62-year record of rainfall and temperature (1951-2013) available for Jerusalem (Israel Meteorological Service -www.data.gov.il/ims),we assessed the annual variability of water balance components at the location of our three soil moisture plots.Rainfall and temperature data from the Jerusalem station were corrected for elevation differences between the Jerusalem station (810 m a.s.l.) and the three plots based on calculated elevation gradients.

Hydrometeorological conditions
The 3 years of high-resolution measurements of precipitation and meteorological parameters revealed considerable interannual variability and a strong elevation gradient, especially in terms of rainfall.Mean seasonal precipitation at the Kafr Malek station (830 m a.s.l.) situated close to the Mediterranean Sea-Dead Sea water divide was 526 mm (380-650 mm), while mean seasonal rainfall at the Auja Village station (270 m b.s.l.) in the Jordan Valley accounted for 106 mm (97-120 mm) leading to seasonal rainfall gradients between 6.4 and 7.2 % per 100 m elevation difference (Fig. 3).Mean rainfall intensity for the single stations was between 0.8 and 1.5 mm h −1 , while maximum intensities exceeded values of 10 mm h −1 at some stations for only a few time intervals during the complete observation period.Convective rainfall events with high intensities presumably from Red Sea Troughs were observed only during a short time period in spring 2011 with cumulative amounts below 40 mm.The mean annual temperature was 7 • C higher at Auja Village whereas relative humidity, wind speed, and net solar radiation were slightly higher at the more elevated station.Stations from the Israel Meteorological Service with longterm records at locations in Jerusalem and the Jordan Valley showed similar characteristics.Three major runoff events resulted from storms with large precipitation amounts and periods of high intensity.Runoff coefficients were smaller than 5 % for single events and less than 2 % for the entire season.

Soil moisture dynamics
Observed soil moisture at all soil profiles (Fig. 4) showed a strong seasonality where the annual course can be divided into distinct phases.At the beginning of the rainy season, the previously dry (8-17 % VWC) soil profile was stepwise wetting up starting from the upper to the lower sensors.During rainfall events with high amounts and intensities, the soil moisture data showed rapid infiltration of water into the deeper portions of the profile.Particularly at plot SM-1, saturated conditions started from the bottom probe close to the soil-bedrock interface, where these conditions persisted for several hours up and to two days.During the strongest rainfall events also upper soil layers reached saturation, however, for much shorter periods (Fig. 4b).At plot SM-3 we found indications of soil saturation from the bottom up to the surface during two events for periods of 8 and 16 h, respectively.At the end of the rainy season, the soil dried out within a few weeks and the soil moisture content further declined at a low rate during the whole dry summer period.3, and their distributions are illustrated in Fig. 5.All models were generally able to reproduce the observed temporal soil moisture patterns with KGE values between 0.82 and 0.94 (Fig. 6).However, differences in predictive capacities at distinct water content levels could be observed, which varied between the single plots (Figs. 6, 7).In general, the model tended to overestimate water contents close to saturated conditions except for deeper sections at plot SM-1 where an underestimation of simulated water contents was observed.
Parameter uncertainty was assessed by simulation of water contents using parameter sets obtained with SCEM after fulfilling the convergence criterion.The 95 % soil moisture confidence interval showed a narrow band around the optimum model (Fig. 8 exemplary for plot SM-1).At all sensors the difference between simulated volumetric water content for the best parameter set and the 95 % confidence interval remained below 4 %, i.e. less than the measurement error of the sensors.
Water retention and conductivity functions from the laboratory MSO experiments are given in Fig. 9.In comparison with the functions from inversely calibrated parameter sets with Hydrus-1D, they show similar characteristics at lower matric potential with an increasing deviation at higher matric potentials.Residual water contents from the MSO analyses were generally higher than the calibrated Hydrus-1D parameter for our soil moisture plots.
Water temperature in a groundwater well near soil moisture plot SM-3 (cf.Fig. 1) indicated five distinct recharge events lowering the mean groundwater temperature from 19 • C by 0.7-4 • C (Fig. 7).The events coincided with the main peaks of modelled percolation from the soil moisture monitoring sites.During these events, mean daily air temperature was less than 6 • C.Although the well was strongly influenced by irregular pumping for water supply (visible as minor water level fluctuations in Fig. 7), major recharge events induced sudden rises of the piezometric water level.

Plot-scale water balance
Modelled fluxes of the various water balance components showed high temporal variability (Fig. 8) and considerable differences in annual values between single years (Table 4).Evaporation and transpiration started shortly after the first rainfall events of the winter season when the water content in the upper soil layer began to increase.Percolation from the bottom of the soil zone only started after the cumulative rainfall during the winter season exceeded a certain threshold.This threshold was found to be ca.240 mm at plot SM-1, 200 mm at plot SM-2, and 150 mm at plot SM-3.This threshold was not a fixed value but varied from year to year    depending on the precipitation distribution over the winter season.In case of the season 2010/2011 with below-average rainfall, evapotranspiration during dry spells reduced the soil water storage, and rainfall amounts of the following events were too low to exceed field capacity and to generate percolation at SM-3.Interception, soil evaporation and transpiration were highly variable during the winter season and depended on the length of dry spells between rainfall events.Evapotranspiration almost ceased within a few weeks after the last rainfall events of the winter season.Mean overall losses through evapotranspiration and interception accounted for 73 % of rainfall.Values slightly above 100 % for the dry year 2010/2011 resulted from elevated moisture conditions at the beginning of the simulation period.Percolation strongly varied from negligible amounts during the dry year 2010/2011 to values ranging between 28 and 45 % of cumulative rainfall during 2011/2012 and 2012/2013, respectively.
The largest proportion of percolation was calculated during a few strong rainstorms.On all three plots, more than 50 % of the total percolation of the 3-year simulation period occurred within a time period of 5-10 days.

Spatial extrapolation of deep percolation
During the hydrological year 2010/2011, cumulative rainfall was below average with totals ranging between 275 and 425 mm (Fig. 10) and a maximum daily amount below 50 mm.In this season, percolation was only simulated for soils with depths up to 60 and 110 cm, respectively.Modelled percolation increased to a maximum proportion of 40 % for shallow soils with depths of 10 cm receiving the highest rainfall input.For the following above-average wet year 2011/2012, seasonal rainfall ranged between 450 and 725 mm.Then simulated percolation rates reached up  to 69 % of rainfall and declined to values close to 0 % only under conditions of lowest rainfall amount and soil depths greater than 160 cm.The third simulated year can be regarded as a year with average rainfall conditions (sums of 400-600 mm).Percentages of percolation were comparable to the previous year although cumulative rainfall was considerably less.This could be attributed to higher rainfall intensities during 2012/2013 when daily rainfall amounts twice exceeded 80 mm and 4 days of rainfall accounted for almost 50 % of the seasonal amount.

Temporal extrapolation of deep percolation
Modelling water balance components for 62 years  resulted in strong differences of simulated seasonal soil water percolation reflecting the high variability of rainfall input (Fig. 11).Mean annual rainfall was calculated for the three plots to range between 408 and 537 mm (standard deviation: 128-168 mm) and mean percolation fluxes between 82 and 150 mm (standard deviation: 93-141 mm).Percolation at the three plots varied between 0 and 66 % of cumulative seasonal rainfall with an average between 16 and 24 %.Other seasonal fluxes varied much less during the sim-ulation period.The coefficient of determination between seasonal sums of simulated percolation and rainfall ranged between 0.82 and 0.88 on the three plots.

Soil moisture dynamics
The observed seasonal dynamics of soil moisture, dominated by short wetting phases during single rainfall events and a rapid decrease after the rainfall season, are comparable with observations reported in other studies in the Mediterranean region (Cantón et al., 2010;Ruiz-Sinoga et al., 2011).At all soil moisture plots, our soil moisture data suggested fast infiltration into deeper sections of the soil profile during rainfall events with high intensities and amounts (e.g.plot SM-1 in Fig. 4b).The time lag between the reaction of the uppermost and the lowermost probe was often less than 2 h, indicating flow velocities exceeding 840 cm per day, despite of high clay content.These fast reactions suggest concentrated infiltration and preferential flow within the vadose soil zone as reported for the Mediterranean region by e.g.(2008).Brilliant Blue patterns from infiltration experiments conducted in the vicinity of our plots highlighted the influence of outcrops on infiltration by initiating preferential flow at the soil-bedrock interface.In the remaining soil, preferential flow was less distinct, but vertical flow velocities of 0.08 cm min −1 suggested also here macropore flow (Sohrt et al., 2014).Hence, a certain fraction of preferential flow is ubiquitous and may further be enhanced by a high stone content in the soil and by bedrock outcrops in the vicinity, as observed particularly at SM-1.In general, bedrock and stones may have multiple effects on infiltration, water retention and water movement in the soil (Cousin et al., 2003).
A noticeable difference between the plots was observed during rainfall events of high magnitude.At SM-1 (Fig. 4b), the bottom probe suggested soil saturation for periods between 2 and 90 h.Durations were apparently linked to the depth of the event precipitation (24-191 mm) and to the duration of the event (16-72 h).The upper probes showed saturation only during the largest rainfall events and for a much shorter duration.Volumetric soil moisture at 10 cm always remained below 30 %.We observed a similar behaviour at SM-3 but not at SM-2.We hypothesize that these phases of saturation were caused by impounded percolation water due to limited conductivity of the soil-bedrock interface.Differences between our plots could be attributed to the variable permeability of the underlying Cenomanian dolomite (SM-1 and SM-3) and Turonian limestone (SM-2).While both formations are known to have high permeability (Keshet and Mimran, 1993), we observed Nari crust (Dan, 1977) in the vicinity of SM-1, which may have reduced hydraulic conductivity.Sprinkling experiments on the same geological material type had already documented soil saturation and subsequent overland flow generation (Lange et al., 2003).

Simulation of the plot-scale water balance
The cumulative distribution functions of the parameters suggested narrow ranges and hence good identifiability for most model parameters (Fig. 5).Nevertheless, measured soil moisture fell outside the 95 % uncertainty band especially during high and low moisture conditions (Fig. 7).This may indicate limitations of our simplified model, which is based on a unimodal pore-size distribution.By definition, our inversely estimated model parameters are effective parameters that describe both preferential and matrix flow.Compared to values of saturated hydraulic conductivity (K s ) of a clay-rich soil matrix from established pedotransfer functions (e.g.Carsel and Parish, 1988), our K s values are high (Table 3).Radcliffe and Šimůnek (2010) analysed data from the Unsaturated Soil Hydraulic Database (UNSODA; Nemes et al., 2001).They found decreasing K s with increasing clay content but also a significant increase in parameter spread.This was attributed to a larger effect of soil structure.This effect will become more evident when moving from the scale of small soil cores to the plot scale, reflecting a common phenomenon of changing parameter values with changing spatial scale (e.g.Blöschl and Sivapalan, 1995).From this perspective, our estimated effective low alpha values describe the small pores of the soil matrix, while the high effective K s values represent the effect of preferential flow.Although clay content and bulk density slightly increased with soil depth at our plots, no clear pattern of calibrated soil hydraulic parameters could be observed.The expected decrease of K s was apparently compensated by other factors such as the observed increasing stoniness of the soil with depth, which could lead to enhanced preferential flow at the soil-rock interface (Sohrt et al., 2014) or by water uptake by plants that was limited to the upper soil zone.Furthermore, persistent saturated conditions during major rainstorms as discussed in the previous section could not be simulated, as a percolation impounding soil-rock interface was not implemented in the model and a free drainage had to be assumed.Still, the conductivity and retention function derived from the MSO experiments showed an overall good agreement with those calibrated with the help of Hydrus-1D and SCEM (Fig. 9).We believe that this is another independent proof of the reliability of our simplified model.As discussed earlier, an increasing deviation of the respective functions with increasing matric potential could be addressed to the different measure-ment scales, where the MSO experiments represent mainly the soil matrix, while the parameters calibrated with Hydrus-1D comprise also preferential flow pathways at the plot scale.
A bimodal pore-size distribution (Durner, 1994) may better represent the heterogeneous pore structure of our clay-rich soil, but at the cost of a larger number of calibration parameters with presumably reduced parameter identifiability and higher model uncertainties.
Originally, Mualem (1976) set the parameter L to a fixed value of 0.5 for all soil types.Later, the physical interpretation of the parameter L representing tortuosity and pore connectivity was increasingly questioned and L was rather treated as an empirical shape factor for the hydraulic conductivity function in the Mualem-van Genuchten model (Schaap and Leij, 2000).Schaap and Leij (2000) observed that fixed positive values of L can lead to poor predictions of the unsaturated hydraulic conductivity and that L was often negative for fine-textured soils.Peters et al. (2011) analysed persistent parameter constraints for soil hydraulic functions and concluded that the conservative constraint of L > 0 is too strict and that physical consistency of the hydraulic functions is Interception SM−1 SM−2 SM−3 q q q q q q Evaporation SM−1 SM−2 SM−3 q q q q q q q q q Transpiration SM−1 SM−2 SM−3 q q q q q q q q q q q q Percolation SM−1 SM−2 SM−3 given for This constraint ensures monotonicity of the hydraulic functions.The requirement of Eq. ( 3) is fulfilled for all L values of the parameter sets shown in Table 1.
Simulated mean evapotranspiration at our plots over the 3-year simulation period accounted for 73 % of rainfall, i.e. very close to the long-term average calculated by Schmidt et al. (2014) for the same area.Our values also fall into the range of Cantón et al. (2010), who derived annual effective evapotranspiration rates of more than 64 % of annual rainfall based on eddy covariance measurements in southeastern semi-arid Spain.Our simulated percolation rates ranged between 0 and 45 % of precipitation (arithmetic mean: 28 %) indicating strong inter-annual variability and a strong dependency on depth and temporal distribution of precipitation.During the entire 3-year period, more than 50 % of overall percolation fluxes occurred during less than 10 days of strong rainfall.These findings are supported by the response of groundwater temperatures observed in a nearby well indicating the arrival of groundwater recharge flux at the water table (Fig. 7).Tracer experiments in a similar setting demonstrated that percolating water can pass the vadose soil and the epikarst at flow velocities of up to 4.3 m h −1 (Lange et al., 2010).Regarding the initiation of percolation at the basis of the soil profiles, we found seasonal rainfall thresholds of ca. 150 mm for the shallow and 240 mm for the deep soil moisture plots.Cave drip studies in the region (Arbel et al., 2010;Lange et al., 2010;Sheffer et al., 2011) measured sim-ilar thresholds for the initiation of percolation through the epikarst (100-220 mm).
In contrast to humid environments, lateral subsurface flow on rocky semi-arid hillslopes rarely develops, since they consist of individual soil pockets that are poorly connected due to frequent bedrock outcrops.Soil moisture seldom exceeds field capacity given that evapotranspiration exceeds precipitation depth throughout most of the year (Puigdefabregas et al., 1998).Furthermore, highly permeable bedrock favours the development of vertical structural pathways in karst areas (shafts beneath dolines and sinkholes).In the epikarst a lateral concentration of the percolation water from the soil zone toward such highly permeable pathways can take place (Williams, 1983).Despite this secondary concentration we can conclude that one-dimensional modelling of the soil water balance is a reasonable approach to understand percolation fluxes and subsequent groundwater recharge.
Nevertheless, we cannot exclude that frequently outcropping bedrock may affect water redistribution by surface runoff or by preferential infiltration along the soil-rock interface.The importance of these effects on percolation rates and groundwater recharge on the regional scale is subject to current research.During heavy storm events, overland flow generation cannot be excluded (Lange et al., 2003), but surface runoff typically accounts for only a few percent of annual rainfall (Gunkel and Lange, 2012).A second limitation of our investigations of plot-scale percolation fluxes is the assumption of an identical vegetation cover at the single sites along the climatic gradient and a constant vegetation cycle throughout years of different seasonal rainfall depths.Al-F.Ries et al.: Recharge estimation and soil moisture dynamics though different plant species and vegetation cycles may alter soil moisture conditions prior to rainfall events, we could show that the event rainfall amount is the main factor that influences percolation rates.

Spatial and temporal extrapolation of deep percolation
Water balance modelling for variable soil depths and rainfall gradients revealed considerable differences for the three winter seasons.(Schmidt et al., 2014).
A high temporal variability in percolation fluxes is also apparent from the long-term modelling of water balance components (Fig. 11).For the 62-year simulation period, we calculated seasonal percolation rates between 0 and 66 % (average: 20-28 %) for our plots.The highest value was modelled for the extremely wet winter season 1991/92 (5 times the mean annual percolation of 150 mm).For a slightly shorter time period, Schmidt et al. (2014) calculated an average recharge rate of 33 % for the Auja spring catchment applying a conceptual reservoir model.They found that recharge of only 5 individual years accounted for one-third of the total recharge of the 45-year period.In our study 7 individual years provided one-third of the total recharge.Furthermore, we compared seasonal percolation of our sites with recharge estimations from perched aquifers feeding small karst springs (Weiss and Gvirtzman, 2007) and the entire carbonate aquifer (Guttman und Zukerman, 1995) (Fig. 12).Although our results plotted within the range of these largescale recharge estimates, we want to emphasize that our calculations display point percolation fluxes.Even in years with below-average rainfall, a certain rise in the groundwater table and spring flow can be observed (season 2010/2011 in Fig. 7h; EXACT, 1998;Schmidt et al., 2014).Then recharge presumably occurs on areas with strongly developed epikarst and shallow or missing soil cover.
Our long-term point calculations suggest substantial differences in percolation fluxes between years of similar rainfall depths.Simulated percolation for plot SM-1 during the seasons 1976/1977 and 2004/2005 accounted for 16 and 35 % of seasonal rainfall, respectively, although both seasons had very similar above-average rainfall (578 and 569 mm).These results are in line with findings of Sheffer et al. (2010) q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q seasonal rainfall (mm) simulated percolation (mm) q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 200 400 600 800 1000 1200  Zukerman, 1995) and three small karst springs emerging from local perched aquifers (Weiss and Gvirtzman, 2007).
and Abusaada (2011) about the importance of temporal rainfall distribution on groundwater recharge.

Implications for recharge in Mediterranean karst areas
The steep climatic gradient, the hydraulic properties and characteristics of the carbonate rocks, the heterogeneous soil cover and the high temporal variability of precipitation on event and seasonal scales are dominating hydrological characteristics in our study area.Similar settings can be found across the entire Mediterranean region.Despite recent advances in the determination of groundwater recharge in karst areas, the assessment of the spatial and temporal distribution of recharge is still a challenge.Modelling approaches including hydrochemical and isotopic data (Hartmann et al., 2013b) require additional information from springs (time series of discharge and water chemistry) for model parameter estimation, which are rarely available.Moreover, the exact delineation of the contributing recharge area is often a problem.Although simulated percolation fluxes from plot-scale soil moisture measurements cannot be directly transferred to the regional, i.e. catchment scale, they can still provide insights into the various processes responsible for the temporal and spatial variability of groundwater recharge as well as information on the relative importance of different process parameters.

Conclusions
This study contributes to the assessment of percolation rates based on soil moisture measurements along a steep climatic Hydrol.Earth Syst.Sci., 19, 1439-1456, 2015 www.hydrol-earth-syst-sci.net/19/1439/2015/ gradient in a Mediterranean karst area.We showed that point measurements of soil moisture together with numerical modelling of the water flow in the unsaturated soil zone may help to understand dominant percolation mechanisms.We found an accentuated annual variability of percolation fluxes and a strong dependency on soil thickness, temporal distribution and seasonal depth of rainfall.To extrapolate our findings, we varied soil depth and climatic input parameters (precipitation and temperature) over ranges observed in our study area.Furthermore, we used a 62-year time series  of climatic input to run our calibrated models.Although our calculations are based on plot-scale measurements, the results closely match long-term observations and their patterns of event and seasonal variability.They also reflect the thresholds for the initiation of groundwater recharge reported by other studies in the same region based on different approaches.Our results suggest that groundwater recharge is most prominent when single rainfall events are strong enough to exceed field capacity of soil pockets over a wide range of soil depths.Hence, the temporal distribution of rainfall has a strong effect on event and seasonal recharge amounts.
Our results corroborate the statement of De Vries and Simmers (2002) about the dependence of groundwater recharge in (semi-)arid areas on high intensity rainfall events.The use of empirical rainfall-recharge relationships can lead to large errors, since recharge rates are sensitive with respect to highly variable rainfall distributions and characteristics, which are most probably affected by predicted climate change in the Mediterranean (Giorgi and Lionello, 2008;Samuels et al., 2011;Reiser and Kutiel, 2012).

Figure 1 .
Figure1.Study area with location of meteorological stations, rain gauges, soil moisture plots (SM-1, SM-2, SM-3) and isohyets of long-term average annual rainfall ( ≥ 20 years) according to data from ANTEA (1998).Coordinates in the detailed map are in Palestinian Grid format.In the upper slope sections of Wadi Auja a local perched spring aquifer formed which is tapped by an abstraction well.

Figure 2 .
Figure 2. Typical hillslopes in the study area.The image shows the plain of Ein Samia with semi-arid climatic conditions, where the valley bottom is used for partly irrigated agriculture and the hillslopes are used as extensive grazing land for goats and sheep.The Wadi Auja ephemeral stream enters the plain from the left (west) and drains to the right (east) in direction to the Jordan Valley.

Figure 3 .
Figure 3. Correlation between average annual rainfall and station elevation for the individual hydrological years during the observation period 2010-2013.

Figure 4 .
Figure 4. Observed volumetric soil moisture at different depths of the three experimental plots during the complete monitoring period (a) and details on the winter season 2011/2012 for plot SM-1 (b).

Figure 5 .
Figure 5. Observed volumetric soil moisture at different depths of the three experimental plots during the complete monitoring period (a) and details on the winter season 2011/2012 for plot SM-1 (b).

Figure 6 .
Figure 6.Observed volumetric soil moisture at different depths of the three experimental plots during the complete monitoring period (a) and details on the winter season 2011/2012 for plot SM-1 (b).

Figure 7 .
Figure 7. Time series of rainfall (a), simulated and observed volumetric water content for soil moisture plot SM-1 (b-e), Hydrus-1D simulated percolation (f), water temperature (g) and piezometric water levels (h) in a nearby groundwater well.The grey shaded area represents the 95 % confidence interval of soil moisture based on model parameter sets obtained using SCEM after fulfilment of the convergence criterion.

Figure 8 .
Figure 8. Simulated daily water fluxes at the single soil moisture plots for the simulation period 2010-2013.

Figure 9 .Figure 10 .
Figure 9.Comparison of the water retention and conductivity functions of the Mualem-van Genuchten parameter sets derived from MSO experiments with those inversely calibrated with Hydrus-1D and SCEM using observed soil moisture time series.

Figure 11 .
Figure 11.Seasonal sums of simulated water balance components for the period 1951-2013 using the calibrated soil hydraulic parameters of the various plots.Rainfall and temperature data were obtained from the nearby Jerusalem central station (http://www.data.gov.il/ims) and corrected for the single locations by applying a simple elevation gradient-based correction factor.

Table 2 .
Parameters and value ranges for Hydrus-1D modelling.

Table 3 .
SCEM-optimized hydraulic parameter sets for the different plots and probe depths.

Table 4 .
Cumulative sums of the simulated water balance components in millimetres and percentages for the three consecutive hydrological years 2010-2013 at the individual soil moisture plots.