Hydrologic response in the karstic and basaltic hydro-geological units of Lake Kinneret watershed

Journal of Hydrology 583 (2020) 124572

Contents lists available at ScienceDirect

Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol

Research papers

Hydrologic response in the karstic and basaltic hydro-geological units of Lake Kinneret watershed

T

⁎

Matan Ben Yonaa,b, Alon Rimmerc,1, Eylon Shamird, , Iggy Litaora a

Tel-Hai College, Upper Galilee 1220800, Israel The Soil Erosion Research Station, Soil Conservation and Drainage Division, Ministry of Agriculture and Rural Development, Bet Dagan 50250, Israel c The Yigal Alon Kinneret Limnological Laboratory, P.O. Box 447, Migdal 14950, Israel d Hydrologic Research Center, 11440 West Bernardo Court, Suite 375, San Diego, CA 92127, USA b

A R T I C LE I N FO

A B S T R A C T

This manuscript was handled by Marco Borga, Editor-in-Chief

The hydrologic responses of two similar in size watersheds in the Lake Kinneret basin are examined: Karstic Hermon watershed and basaltic Meshushim watershed. The streamflow properties at these watersheds are different in their surface flow to baseflow ratio, runoff coefficients and response to extreme rainfall events. The HYdrological Model for Karst Environment (HYMKE) model, developed to simulate daily flow in karst formation, was modified to address the basaltic formation at the Meshushim watershed. Water fluxes simulations were used to compare the hydrologic processes between the watersheds. The twice-higher annual precipitation in the Hermon watershed, generates annual streamflow that is five times more than in the Meshushim. This amplification of the rainfall-streamflow ratio between the watersheds is because larger fraction of the rainfall is lost to evapotranspiration at the basaltic Meshushim watershed and larger portion of the rainfall percolates to the karstic aquifer through preferential flow and appear as baseflow at the Hermon watershed. During large rainfall events when soil water content is near saturation, lower rainfall amount over the Meshushim watershed generates larger peak flow than the Hermon watershed. In both watersheds, surface flow was generated only when the soil water content is at or near saturation conditions. The streamflow dependency on soil water content, suggests high sensitivity of streamflow generation to the intra-annual variability of precipitation. Understanding the hydrologic response of the different sources of Lake Kinneret improves the ability to forecast its inflow.

Keywords: Lake Kinneret basin Karst hydrologic response Basalt hydrologic response Runoff coefficient

1. Introduction Lake Kinneret (known also as Sea of Galilee), located in the northcentral part of the Jordan Rift Valley, is the largest freshwater lake in Israel. The surface flow into Lake Kinneret comes from the Upper Jordan River watershed (~1600 km2) and streams that directly drain into the Lake (~965 km2). Lake Kinneret watershed consists of four distinctive geologic formations that vary in their hydrologic response to precipitation (Rimmer and Givati, 2014; Sade et al., 2016): (1) The Jurassic mountainous karst of Mount Hermon characterized by faults and semi impermeable layers; (2) The basalt plateau of the Golan Heights; (3) The Cenomanian-Turonian carbonaceous karst of the Eastern Galilee Mountains; and (4) The flat alluvial and peat formations of the Hula Valley (Fig. 1). The development of a hydrologic model for the entire watershed to simulate inflow into Lake Kinneret is an on-going research challenge. Recognizing the hydrologic complexity of the Lake’s watershed, the

⁎ 1

HYdrological Model for Karst Environment (HYMKE) (Rimmer and Salinger, 2006) was developed to address the hydrologic complexity caused by the karstic formation of the Hermon Mountain, which its tributaries provide the majority of the perennial contribution to the Upper Jordan River that feeds into Lake Kinneret. Since its development, HYMKE was used in many hydrological studies of the Lake Kinneret watershed (e.g. Samuels et al., 2009; Samuels et al., 2010; Rimmer et al., 2011; Hartmann et al., 2013; Gilboa et al., 2015; Rimmer and Hartmann, 2014; Sade et al., 2016). In addition, the model has been used as an operational tool by the Israeli Hydrological Service for hydrological nowcasting of daily streamflow in the tributaries of the Hermon Mountain. Although HYMKE provides for a well-founded hydrologic modeling strategy, it requires enhancements in order to adequately represent the other hydrologic formations over the Lake Kinneret basin. In this study, we compare the hydrologic processes in two relatively small watersheds with similar sizes located in the Mount Hermon and the Golan Heights hydrologic units.

Corresponding author. Deceased.

https://doi.org/10.1016/j.jhydrol.2020.124572 Received 24 August 2019; Received in revised form 20 November 2019; Accepted 10 January 2020 Available online 20 January 2020 0022-1694/ © 2020 Elsevier B.V. All rights reserved.

Journal of Hydrology 583 (2020) 124572

M. Ben Yona, et al.

Fig. 1. The Upper Jordan River basin (Top), elevation gradient of the Meshushim watershed and the Hermon watershed (bottom).

hydrologic unit, is the largest source of water to the Upper Jordan River with an average annual flow at the outlet of the Hermon stream of 111 × 106 m3/year (Fig. 1). The Hermon Mountain range is an elongated anticline of mostly karstic limestone of Jurassic age with thickness greater than 2000 m (Michelson, 1979; Gilad, 1980; Gilad and Schwartz 1978). The karst formation with its large cavities and interconnected network of channels is characterized by preferential flow pathways that convey surface water to the deep aquifer (e.g. Gilad, 1980). The soil in the Hermon watershed is classified as Terra Rossa that is typically overlaid on hard limestone rock (Singer, 2007). The land surface of the watershed is pitted with karst features such as dolines and sinkholes (Frumkin et al., 1998). Land cover in the lower elevation (below 1300 m) is classified as drought resistant evergreen shrublandforest dominated by Palestine oak (Quercus calliprinos). The higher part of the watershed (1300–1900 m) is covered by a Xero-Montane open forest composed of winter deciduous trees and shrubs that can withstand low temperatures and strong wind. The highest elevation (above 1900 m) is classified as Subalpine Tragacanthic belt typical of arid alpine regions and composed of shrubs and short herbaceous species

Following a detailed analysis of the watersheds’ responses, we modified the HYMKE to simulate streamflow at the basaltic formation of the Golan Heights and used it to examine the hydrologic differences between the two watersheds. Our objective in this study is to understand the dominant hydrologic processes differences between the two watersheds and to develop a robust and reliable hydrologic model that captures the nuance differences in the hydrological processes between the two hydrologic units. This understanding serves the broader objective of developing a hydrologic model for the Kinneret basin that appropriately captures the hydrologic response of the different contributing tributaries. 2. Experimental setup 2.1. Study area The study was conducted in the perennial streams of the Hermon and Meshushim watersheds, located at the Hermon Mountain and Golan Heights hydrologic units, respectively. The Hermon watershed (147 km2, 200–2200 a.m.s.l), which is located in the Hermon Mountain 2

Journal of Hydrology 583 (2020) 124572

M. Ben Yona, et al.

Table 1 Data that were used for the model development.

Rain Streamflow Potential evapotranspiration

Meshushim

Hermon

Source

Average of five stations at the Golan Heights (Fig. 1) Streamflow data from the Meshushim station Pan evaporation data from three stations (Fig. 1)

From Rimmer and Salingar, 2006) Streamflow data from the Hermon station

Israeli Meteorological Service Israeli Hydrological Service Israeli Meteorological Service

over Northern Israel (Ziv et al., 2014). In general, the rainfall spatial distribution over Lake Kinneret watershed can be characterized by its north–south decreasing gradient attributed to increasing distance from the main Mediterranean cyclone track; west–east decreasing gradient attributed to depleted moisture further away from the Mediterranean Sea, the source of the air moisture; and increase rainfall with elevation attributed to orographic enhancement (e.g. Kutiel, 1987; Shamir et al., 2016; Shentsis et al. 2018).

(Auerbach and Shmida, 1993; Danin, 1988). In addition to intermittent contribution of surface runoff, the main contribution of flow to the Hermon perennial stream is from a series of springs to include the Banias, Sa’ar, Kezinim, and Si’on. During the winter, a snowpack is formed at the high elevation (above 1400–1900 m a.m.s.l) to completely ablate during the spring (March–May) (e.g. Shmida, 1977; Gilad and Bonne, 1990; Sade et al., 2014). The Meshushim watershed (160 km2, 150–1000 a.m.s.l) drains the central Golan Heights into Lake Kinneret (Fig. 1). The Golan Heights is a basalt-covered plateau that gently slopes from north to south, with a steep slope in the western edge of the plateau as it approaches the Hulla Valley in the Upper Jordan River and Lake Kinneret. The watershed is covered by basalt of Pliocene-Pleistocene age (Mor, 1986; Heimann, 1990) from the Bashan geological group (Mor, 1986). The shallow channel at the Meshushim’s headwater deepens as it flows downstream towards the lake Kinneret. The average flow at the outlet of the Meshushim stream is 27.5 × 106 m3/year. Most of the flow at the Meshushim stream is from winter precipitation events that trigger surface runoff. However, the perennial flow at the Meshushim is attributed to a relatively shallow aquifer at the headwater that seeps into the channel at the interface with the channels’ walls (Dafny et al., 2003, 2006). The streamflow at the outlet of the Meshushim watershed represents a natural flow regime with little disturbance by upstream reservoirs and agricultural practices (Shtober-Zisu and Inbar, 2005). The Vertisols soil at the lower portion of the watershed had originated from eroded soil of the higher portion of the watershed. The upper portion of the watershed has shallower soil that consists of Alfisols (Brown and Red Mediterranean soils) (Singer, 2007). The basalt soil in the Golan Heights has low hydraulic conductivity (0.1–0.3 m/ hour), causing low infiltration rate and surface runoff generation that is highly dependent of rainfall intensity, which also yield a high erosivity level of the soil (Singer, 1987). The higher parts of the Meshushim watershed (higher than 500 m a.s.I.) are covered with oak woodland that is dominated by Palestine oak. The lower part of the watershed is mainly barren and rocky with shrub-steppe cover dotted with Tabor oak (Quercus ithaburensis) (Danin, 1988). Some cultivated areas exist mainly in the lower portion of the basin where the soil is deeper. The low hydraulic conductivity (less than 1 m/day) of the basalt rocks yields a slow recharge rate of the groundwater aquifer (Dafny et al., 2003). In addition, the basaltic rock column is composed of repeated layers of basalt rocks from different eruptions separated by clayey paleosols material. This composition causes the vertical hydraulic conductivity to be substantially lower than the horizontal hydraulic conductivity (anisotropy of low hydraulic conductivity) to form many small springs that are fed from the limited perched aquifers (Dafny et al., 2003). The climate of the Lake Kinneret watershed is classified as Mediterranean with average annual rainfall ranges from 1200 mm/year in the north to less than 500 mm/year in its southern part. Most precipitation events, which occur during the winter season (November–March), are triggered by the so-called Cyprus Low synoptic system. These are low-pressure systems that are formed near Cyprus that move eastward towards the Eastern Mediterranean shore. The passage of these systems is accompanied by westerly winds that transport moist air over the warm Mediterranean Sea, where it becomes increasingly moist and conditionally unstable before it is transported

2.2. Observations Daily mean areal rainfall was estimated for each watershed using the observation network of the Israeli Meteorological Service (IMS) (Table 1). For the Hermon watershed, we adopted the 1976–2001 interpolated time series from Rimmer and Salingar (2006), who used a seasonally varying linear regression with elevation as the predictor to estimate the mean areal precipitation over the entire basin. Such interpolation is required since the north most rain gauge is located at the Golan Heights at an elevation of 960 a.m.s.l, and no reliable observed precipitation record is available for the high elevation of the Hermon Mountain (e.g. Shamir et al., 2016). Using a denser network of rain gauges for the Golan Heights, the mean areal rainfall for the Meshushim watershed was estimated for 1976–2013 as the average of the qualitycontrolled records from five observation stations (Fig. 1). Daily streamflow discharge data are available for the outlet of both watersheds from hydrometric stations operated by the Israeli Hydrologic Service. Daily climatological potential evapotranspiration (PET), the amount of evapotranspiration that would occur if sufficient water source are available, were calculated by assigning a Sine function that describes the seasonal variability to long-term (1970–2000) daily pan evaporation observations as suggested by Rimmer and Salingar (2006). 2.3. Hydrologic response analysis The differences in the hydrologic regimes between the two watersheds are exemplified in Fig. 2. The cumulative distributions of the total water year (1 October–30 September) rainfall and the mean average discharge over the Hermon watershed are significantly higher than the Meshushim watershed (Fig. 2a, b). Note, that throughout the manuscript the term ‘annual’ denotes water year. The average differences between the watersheds are noteworthy since the Hermon is two to five times larger than Meshushim when comparing the annual rainfall and daily average annual discharge, respectively. To further understand rainfall-runoff response we analyzed the cumulative distributions of the daily rainfall and daily streamflow annual maxima series (i.e. the series of maximum daily values during each of the water years) (Fig. 2c, d). As expected, the extreme daily rainfall values in the Hermon are larger than the Meshushim. However, the discharge maxima series in the Meshushim are larger than the Hermon. This suggests that smaller rainfall events trigger larger flow events in the Meshushim. The differences in the transformation of rainfall to streamflow in the two watersheds, are further examined in Fig. 3 where the total observed flow and the estimated baseflow component are shown. The baseflow component, the slowly changing discharge that is contributed to the streamflow with a considerable delay following 3

Journal of Hydrology 583 (2020) 124572

M. Ben Yona, et al.

baseflow exhibits a distinct seasonal cycle ranging between 1.1–4.6 m3/ sec during March and November, respectively. During large flow events in the Hermon, the surface flow is up to about 7 times of its baseflow (Fig. 3b). The baseflow as measured at the Meshushim outlet is only about 17.5% (4.8 × 103 m3/year) of the total streamflow. The streamflow at the hydrometric station is mainly attributed to surface flow that occurred after the soil water content has reached saturation. There is almost no observed variation in baseflow during the winter and the daily baseflow ranged between 0.08 and 0.26 × 103 m3/sec in March and November, respectively. The surface flow increased by 250 times more than before the flooding events. (Fig. 3a). The annual runoff coefficient (RC) – the ratio of water year flow and mean areal rainfall volumes – are shown for the watersheds in Fig. 4 as a function of the water year rainfall. Compared with the Hermon watershed the RC for Meshushim is smaller and has a noticeable positive dependency on the annual rainfall. The annual rainfall explains ~60% and ~40% of the RC inter-annual variability for the Meshushim and Hermon, respectively. These results of RC dependency on the rainfall, imply that at the Hermon streamflow generation is dominated by multiyear climate conditions, while at the Meshushim the rainfall pattern of the current year is the dominant factor. A key reason for the annual rainfall being a limited predictor of the RC is the watersheds capacity to absorb and retain water in the soil. At the beginning of the water year (October 1st), after the long warm and dry summer season (June–September), the soil is completely dry. In both watersheds, surface flow events were not observed before, on average, 250 and 220 mm of rainfall occurred since October 1st, which is about 18% and 34% of the annual totals in the Hermon and Meshushim, respectively (Table 2). The large variability of rainfall amount that is required for the first surface runoff event is attributed to the characteristics of the sequence of rainfall events at the beginning of the wet season. During dry spells between rainfall events, evapotranspiration and percolation processes deplete the water content of the soil and therefore the frequency and duration of these dry spells control the surface runoff generation.

Fig. 2. Meshushim and Hermon watersheds cumulative distributions of (a) annual total rainfall (mm/year); (b) annual average daily discharge (m3/sec); (c) annual maxima of daily rainfall; and (d) annual maxima of average daily discharge (m3/sec). All the cumulative distributions were calculated for 1976–2001 water years.

rainfall events, was estimated herein using a recursive digital filter method to separate high from low-frequency signals. We used the Eckhardt’s (2005) method, in which the streamflow is described as follows: ∗ Q∗j = QBj + QSj∗ .

(1)

Where Q* is the total observed streamflow, Q*S is the direct surface runoff, Q*B is the baseflow and the j subscript represents the time steps. To calculate the baseflow we used a robust filter with two empirical parameters (α and β) (Eckhardt, 2005; Rimmer and Hartmann, 2014): ∗ QBj =

∗ ∗ [(1 − β ) αQBj − 1 + (1 − α ) βQj ]

(1 − αβ )

∗ ; QBj ⩽ Q∗j .

(2)

The β is a scaling parameter and α corresponds to the recession constant as seen in Eq. (3) for the special case that β is set to zero: ∗ ∗ QBj = αQBj −1

3. Hydrologic model We used the HYMKE for the Hermon watershed as implemented in Rimmer and Salingar (2006) (the HYMKE formulation is in Appendix A). The HYMKE consists of four modules: the surface layer (0), the vadose-zone layer (1), groundwater aquifer (2) and surface flow (3). Because of the large difference in depth to groundwater between the watersheds (greater than 2000 m in the Hermon watershed (Gilad, 1980) and 100–600 m at the Meshushim watershed (Dafny et al., 2006)) in addition to the differences in percolation characteristics, the

(3)

For the Hermon, we set the α and β values as recommended by Rimmer and Salingar (2006). For the Meshushim, the α parameter was sequentially estimated by optimizing Eq. (3) and then selecting a β parameter as outlined in Rimmer and Hartmann (2014). In the Hermon watershed the baseflow component comprises of about 75% of the total streamflow (84 × 106 m3/year) and the daily

Fig. 3. Observed streamflow and baseflow at the Meshushim (a) and Hermon (b) watersheds. 4

Journal of Hydrology 583 (2020) 124572

M. Ben Yona, et al.

Fig. 4. Annual runoff coefficient (RC) as a function of the annual rainfall.

experiments, the calibration was carried in two-steps. First, we selected the parameters that optimize the surface flow and second; fixing the surface flow parameters, we selected parameters that optimize the simulation of baseflow. We used the following two goodness-of-fit coefficients to evaluate the model performance: (1) The coefficient of determination (R2), which represents the linear association between two time series, and (2) the normalized benchmark simple model (BSM) efficiency coefficient, which indicates the model's ability to predict the streamflow relative to a simple baseline model (Schaefli and Gupta, 2007; Rimmer and Hartmann, 2014). The modified BSM substitutes the climatology component in the coefficient’s denominator with two simple models (benchmark models) one for baseflow and one for surface flow (Eq. (4)).

baseflow recession at the Meshushim watershed is better represented by eliminating the vadose zone layer and implementing two parallel groundwater linear reservoirs (Rimmer and Hartmann, 2012). Thus, the groundwater module for the Meshushim was modified to simulate the sum of two linear reservoirs (Fig. 5 a). The first linear reservoir represents the regional aquifer and the second represents perched aquifers that lay above the paleosols clay layers. Parts of the Meshushim watershed feed the Golan Heights side springs that are downstream of the hydrometric station and do not contribute water to the baseflow of the Meshushim stream (Dafny et al., 2006). The model coefficients that distribute the surface water entry to the two linear groundwater reservoirs (Table 3) are summed to 90%, while the remaining 10% is assumed to feed the side springs. A modified procedure to calculate the daily actual evapotranspiration (EAj mm) was introduced in this HYMKE version. In this modified procedure the EAj is calculated as a function of the water content of the upper soil. The formulation of this modification is summarized in Appendix A.

n

BSM = 1 −

∑ j = 1 (Qobsj − Qsimj )2 n

∑ j = 1 (Qobsj − Qbmj )2

(4)

where Qobsj and Qsimj are the observed and simulated flow at the j’s time step. The Qbmj is the benchmark model simulation, calculated as described in Eqs. (5) and (6). The BSM efficiency values range from −∞ to 1 and a BSM value of 1 corresponds to a perfect match between the simulated and observed data. An efficiency of 0 (BSM = 0) indicates that the model predictions are as accurate as the alternative simple baseline model and efficiency values of less than zero (BSM < 0) indicate that the simple model serves as a better predictor than the model. The benchmark model for the baseflow consisted of a cyclic function (Sine) that calculates the baseflow as a function of the day in the year according to the following equation:

3.1. Parameters estimation We used a split sample methodology in which the parameters were calibrated using the 1976–1996 data and evaluated on the remaining time series (1997–2013 for Meshushim and 1997–2001 for the Hermon). A description of the parameters and their selected values are given in Table 3. To estimate the parameters we used a multistep automatic calibration approach, an approach that combines the strengths of the manual and automatic calibration strategies (Hogue et al., 2000). The observed discharge data was separated to baseflow and fast flow components, as shown in Eq. (2). The model parameters were divided into parameters that dominate the simulations of the baseflow and the fast flow. The parameters of each group were estimated by minimizing the error between the simulated and observed time series for both the baseflow and surface flow. We used the Microsoft Excel Solver (Premium Solver Platform, 2010) to select parameter sets that minimized the sum of the squared residuals. Using a comprehensive and iterative trial-and-error

Qsimj = b [1 + a ∗ sin(λ ∗ (JD + ω))]

(5)

Were JD is the day of the year, a and b are scaling coefficients (1000 m3), λ is the angular frequency (rad), and ω is the phase shift. The benchmark model for surface flow calculates the flow with the daily rain, according to the following equation:

Table 2 Inter-annual average and range in parenthesis of rainfall and streamflow characteristics. Surface runoff event is defined as streamflow event greater than 0.57 m3/s, (~50,000 m3/day).

Hermon Meshushim

Number of days from October 1st to the first surface runoff event.

Cumulative rainfall (mm) from October 1st to the first surface runoff event.

Cumulative annual rainfall (mm)

63 (36–98) 90 (34–159)

252 (134–407) 220 (30–327)

1383 (940–2142) 645 (382–1140)

5

Journal of Hydrology 583 (2020) 124572

M. Ben Yona, et al.

Fig. 5. Schematic description of the conceptual hydrologic model: (a) Meshushim model and (b) Hermon model. Table 3 The model parameters.

Table 4 Surface flow and baseflow simulation evaluation as measured by the coefficient of determination and benchmark simple model (BSM) efficiency coefficient.

Surface Layer Parameter

Symbol

Units

Meshushim

Hermon

Hydraulic conductivity

KD (θ)

m ∗ day−1

0.00004

0.004

Residual moisture content

θr

m−3/ m−3

0.05

0.16

Saturation moisture content

θs Z N

m−3/ m−3 M –

0.35

0.4

0.55 3

1 4

Thickness of the topsoil layer Empirical constant (Mualem and Dagan, 1976)

Hermon

Calibration Validation

Meshushim

Calibration Validation

Surface Flow Fraction contributed to surface runoff Decay factor of surface flow

as

–

0.8

0.45

K3

day

1.6

7

BSM

0.65 0.93 0.48 0.61 0.64 0.86 0.55 0.75

0.33 0.44 0.22 0.42 0.35 0.21 0.31 0.16

0.72 0.96 0.68 0.65 0.44 0.82 0.32 0.61

4.1. Model evaluation

ain1/ ain2

–

0.2 / 0.7

K1

day

339.3

K2

day

57.2

K2

day

α β

– –

A summary of the model evaluation results for both the calibration and validation durations are given in Table 4. As mentioned above, both watersheds were calibrated on 1973–1993 data and evaluated for 1976–2015 and 1976–2001 in the Meshushim and Hermon, respectively. It is seen that the simulations for the Hermon performed better than for the Meshushim and the simulated baseflow outperforms the simulated surface flow. The BSM indices indicate that the model simulations outperformed the simulations of the baseflow and surface flow models in both watersheds.

8.25

107.8

Baseflow separation parameters Baseflow filter Baseflow filter

R2 (benchmark model)

4. Results

The vadose zone layer and groundwater aquifer Fraction distribution of entry to groundwater aquifer Decay factor of groundwater flow (regional aquifer) Decay factor of groundwater flow (perched aquifers) Decay factor of flow in the vadose zone

Surface flow Baseflow Surface flow Baseflow Surface flow Baseflow Surface flow Baseflow

R2 (new model)

0.995 0.2

0.995 0.8

4.2. Hydrologic response

Qsimj =

Q¯obs ∗ A ∗ Rj − lag A ∗ ¯R¯ obs

The hydrologic response differences between the two watersheds were further examined using the hydrologic model simulations. The average fraction of the simulated water fluxes of evapotranspiration, surface flow flux and are summarized in Table 5. The sum of these three fluxes in each water year is about equal to the annual rain. The small deviation of the total fluxes from the annual rainfall is attributed to model’s numerical error of –0.001% and −0.036% in the Meshushim and Hermon, respectively. In this region, because we assume that by the end of the long and dry summer (~June–September) the soil is dry, year-to-year differences in soil water storage at the beginning of the water year (October) are negligible. The differences in the partitioning of the water fluxes between the two watersheds are apparent. As expected, since the rainfall in the

(6)

where Q¯obs and Robs are the average annual observed streamflow and rainfall, respectively, j represents daily time steps, Rj − lag is the daily rain, lag represents the duration (days) between rain events to peak flow, and A is the drainage area of the watershed. In both watersheds, the response of the streamflow to rainfall is within a day and therefore the lag parameter was set to one day.

6

Journal of Hydrology 583 (2020) 124572

M. Ben Yona, et al.

5. Discussion

Table 5 Average annual water fluxes. mm Meshushim 1976–2013

Hermon 1976–2001

Rain Evapotranspiration Surface flow Groundwater aquifer Rain Evapotranspiration Surface flow Groundwater aquifer

611 470 111 30 1328 742 144 441

Surface runoff is usually generated when either precipitation rate exceeds the infiltration rate (‘‘infiltration excess“), or precipitation falls on saturated land surface (‘saturation excess’). The daily contribution surface flow to the total flow in the studied watersheds is clearly associated with the seasonal wetting and drying dynamics of the upper soil. At the beginning of the hydrologic year (October), coming out of a long dry and warm summer, when the soil is dry, no surface runoff was observed in both watersheds, regardless of the observed rainfall intensity. Surface runoff was observed in the Hermon and Meshushim watersheds only after mean areal precipitation of about 250 and 200 mm had accumulated since October 1, respectively. This suggests that runoff generation in these watersheds is dominated by saturation excess processes. Saturation access as the major runoff generation process was previously reported in several other Mediterranean climate studies (e.g. coastal Israel (Peleg et al., 2015); plot scale experiments in Jerusalem and the Galilee mountains (Lavee et al., 1998); Galilee Mountains (Ben-Zvi, 1988); a mountainous watershed in the Pyrenees (Gallart et al., 2002; Latron et al., 2008)). We note however, that infiltration excess processes may occur in temporal and spatial scales that are finer than the analysis conducted in our study (e.g. Ries et al., 2017). During the winter months (December-February), the soil is becoming saturated and rainfall events often trigger surface flow. During this period, groundwater levels rise from percolation through preferential flow and matrix flow processes to increase the baseflow in the outlet of the watersheds. During the fall and summer, the soil is drying and surface flow events are infrequent while baseflow is gradually declining. This dependency of streamflow on the soil water content of the upper soil was found to be stronger in the basaltic Meshushim watershed because of the larger fraction of surface flow that is contributed to the total streamflow. In the mostly karstic Hermon watershed, the percolation flux contributes to baseflow regardless of the soil saturation conditions at the upper soil. The implication of such dependency is that the rainfall temporal and intra-annual variability is a key factor that controls the surface flow generation and total streamflow. In seasons with long dry spells between events, the upper soil desiccates by evaporation and percolation processes and surface runoff is being suppressed. This dependency on intra annual variability is demonstrated by Dubrovský et al. (2014) who showed that decrease in the numbers of rain events required larger rain events to saturate the soil and generate surface flow. In addition, the decrease in number of rainy days increased the total evapotranspiration, increased the percolation flux and decreased the total streamflow.

%

77 18.13 4.9 55.9 10.9 33.2

Hermon is more than double of the Meshushim watershed, the absolute values in the Hermon are larger for all fluxes. However, in relative terms the only flux with larger fraction in the Hermon is water percolation into the aquifer. Annual distribution of the flux partitioning (Fig. 6) shows the large inter-annual variability of the surface flow in the Meshushim and the relatively low variability in the Hermon. The evapotranspiration fluxes in both watersheds have small changes from year to year. However, in both watershed – more pronounced in the Meshushim – during dry years the relative fraction of evapotranspiration is higher. The seasonal wetting and drying dynamics of the soil is a key factor that controls both the generation of baseflow and surface flow. Increase water content of the soil increase the groundwater recharge rate due to increase of the percolation rate. It is only when the top soil reaches saturation, the excess water becomes available for percolation to the sub surface and surface flow is generated. During wetter years, the soil water content is higher for longer duration, which yields higher efficiency in transformation of rainfall to streamflow. The dynamic variability in soil water content during the year is similar in the two watersheds because of the strong correspondence in precipitation events between the two watersheds. This strong correspondence is realized in the magnitude, the timing, and the number of rainfall events. The soil water content and the rainfall of two typical dry and wet hydrologic years (1986–1987 and 1985–1986, respectively) are shown for the two watersheds (Fig. 7). Although the dynamic of these two years were different, the two watersheds exhibit a similar soil water content dynamic. This similarity in the soil water dynamic is shown in spite of the Hermon watershed experiencing more rainfall (Fig. 7). The surface flow in both watersheds (observed and simulated) appeared only when the soil had reached saturation and during the period in which the soil had been saturated. In both years and watersheds, rainfall events that occurred before and after the saturated soil period, did not trigger streamflow in the watersheds’ outlets.

Fig. 6. The annual partitioning of the evapotranspiration, groundwater recharge and surface fluxes in the Meshushim and Hermon watersheds. 7

Journal of Hydrology 583 (2020) 124572

M. Ben Yona, et al.

Fig. 7. The wetting and drying dynamic of the soil during a typical wet and dry hydrologic years (1986–1987 and 1985–1986, respectively).

implies that in order to accurately capture the hydrologic response of the watersheds it is important to describe, in addition to the magnitude of the precipitation events, the precipitation intra-annual temporal characteristics. These precipitation features should be considered and analyzed in climate models projections, in order to appropriately projects future water resources inflow to Lake Kinneret.

Future climate study of the region projected a decrease of total annual precipitation, decrease in consecutive rainy days, decrease in the number of wet days and increase in extreme rainfall events (Samuels et al., 2017; Sade et al., 2016; Peleg et al., 2015; Rimmer et al., 2011). Although the impact of the projected future precipitation regime on the water balance of the two watersheds will be manifested differently, it is likely that surface runoff will decrease and the fraction of the evapotranspiration flux will increase in both watersheds. In the Meshushim watershed, a larger fraction of the rainfall will be lost to evapotranspiration, which will substantially decrease the proportion of rainfall that transforms into streamflow. In the Hermon watershed because of the dominance of the percolation flux, the proportion of rainfall that is transformed to streamflow will be less affected. In this study, the watersheds were treated as uniform units without accounting for their internal spatial variability. Their spatial variability can be expressed in the variability of the hydrologic input, and variability in the spatial characteristics of the land–surface. Moreover, the model parameters were selected to optimize the integrated response of the watersheds. Despite these potential sources of uncertainty, the model simulations describe well the dominant processes in the watersheds. This good performance enables the use of the model simulations to diagnose the unique hydrologic responses of the two watersheds and better understand their differences. In this study, the watersheds were treated as uniform units without accounting for their internal spatial variability. Their spatial variability can be expressed in the variability of the hydrologic input, and variability in the spatial characteristics of the land–surface. Moreover, the model parameters were selected to optimize the integrated response of the watersheds. Despite these potential sources of uncertainty, the model simulations describe well the dominant processes in the watersheds. This good performance enables the use of the model simulations to diagnose the unique hydrologic responses of the two watersheds and better understand their differences. Our work highlights the strong dependency of streamflow generation on the soil water content of the watersheds. This dependency

6. Conclusions We examined the hydrologic response differences between two similar in size watersheds in a karstic (Hermon) and a basaltic (Meshushim) geologic formations. These watersheds, which are located in Lake Kinneret basin, represent different hydrologic regimes. The hydrologic response differences between the watersheds are expressed in the surface flow to baseflow ratio, the runoff coefficient and their response to extreme rainfall. The HYMKE model, which was originally developed to simulate daily surface flow at the Hermon Karstic geologic formation was modified to capture the hydrologic response of the basaltic geologic formation at the Meshushim watershed. The simulated water fluxes among the various model compartments were examined to infer on the dominant processes in the watersheds. Although the prevailing precipitation events in both watersheds are commonly generated by the same weather synoptic system, the annual precipitation in the Hermon watershed is about twice as much as in the Meshushim watershed. This difference between the watersheds is magnified to five times when the annual streamflow are compared between the basins. The magnifying of the streamflow ratio between the two watersheds is attributed to a higher fraction of evapotranspiration in the Meshushim watershed and larger portion of rainfall that percolates to the aquifer and contributes to the stream through baseflow in the Hermon watershed. Although lesser rainfall at the Meshushim, during the winter when the soil water content is at or near saturation, surface flow at the Meshushim is larger than the Hermon. This is because the large portion that percolates through preferential flow to the aquifer in the Hermon 8

Journal of Hydrology 583 (2020) 124572

M. Ben Yona, et al.

watershed and the longer duration of the flood hydrograph, which attenuates and reduces the peak flow. In both watersheds, surface flow is generated only when the soil water content is at or near saturation conditions. This dependency of streamflow on soil water content, makes the generation of surface flow highly sensitive to the intra annual variability of rainfall.

(1st author) who was initially mentored and supervised by Dr. Alon Rimmer (2nd author) who passed away in 2018. Deeply saddened, we dedicate this work to Dr. Rimmer our mentor, colleague, and friend. Alon is an internationally recognized scientist in hydrology and physical limnology who published more than 80 peer reviewed articles. He was a leading expert on the complex hydrology and hydrogeology of the Lake Kinneret watershed. Aside from his scientific activity, Alon was a relentless advocate for environmental protection issues across Israel in general and in Lake Kinneret and its watershed in particular. We attempted to complete this work adhering to Alon’s high scientific integrity and standards. We hope Alon would have approved of this manuscript had he had the opportunity to read it. Thanks are extended to Tel Hai College and the Israel Oceanographic and Limnological Research for the financial support of the 1st author.

Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements This study is based on a master research project by Matan Ben Yona Appendix A

In this appendix the formulation of the HYdrological Model for Karst Environment (HYMKE), which was developed by Rimmer and Salingar (2006) is provided. The module of the surface layer is governed by the mass balance equations:

θPj ; if (θPj < θS ) θj = ⎧ ⎨ θS ; if (θPj ⩾ θS ) ⎩ where θPj = θj − 1 +

qINj − qOUTj

(A1)

Δz −3

here θ is the soil water content (m m ), θP the “potential” soil water content, θS indicates saturation, Δz the thickness of the topsoil layer (m), and ‘j’ is the daily index. The daily flux into the surface layer qIN (m) was defined as the daily rainfall: 3

qINj = 0.001 × Rj ;Rj ⩾ 0

(A2)

where Rj is the daily rainfall time series in (mm) and the 0.001 is used to transform the units from mm to m. The daily flux that leaves the surface layer qOUT (m) is comprised of the following two components:

qOUTj = qDARCYj + EAj

(A3)

It is assumed that the daily soil water content from the surface soil to a depth Δz is uniformly distributed. Thus, qDARCY describes a “unit gradient”, in which the vertical flux, defined by Darcy law is reduced to:

qDARCYj = −KD (θj ) KD (θj ) =

⎧ KDsat

(

θj − θR θS − θR

n

);

θR ⩽ θj ⩽ θS

⎨K ; θ = θ S ⎩ Dsat j

(A4) −1

−3

here KD(θ) and KDsat are the daily and saturated hydraulic conductivity of the soil (m day ), θ r (m m ) is the residual soil water content of the topsoil layer, and the value of the exponent n ranges between 3 and 4 (e.g., Mualem and Dagan, 1978). A modified procedure to calculate the daily actual evapotranspiration (EAj mm) is introduced in this study. The EAj is calculated as a function of the water content in the upper soil. 3

AWj − 1 ⎞ E Aj = E Pj ∗ ⎛ ⎝ AWC ⎠ ⎜

⎟

(A5)

where EPj (mm) is the daily potential evapotranspiration, AWj (mm) is the water content in the upper soil layer, AWC (mm) is the maximum water capacity of the upper soil layer and j is the daily index. We define: (A6)

AWC = 1000 ∗ Δz ∗ θS AWj = 1000 ∗ Δz ∗ θj

(A7) 3

−3

3

Where Δz (m) is the thickness of the topsoil layer, θj (m m ) is the state of the volumetric soil water content and θS (m m saturated soil water content. The actual evapotranspiration is then calculated as followed:

θj − 1 ⎞ EAj = EPj ∗ ⎛ ⎝ θS ⎠ ⎜

−3

) is the volumetric

⎟

(A8)

Note that when the soil water content at the surface soil layer reaches saturation (θ S), the difference Δz × (θPj − θS ) in Eq. (A1) forms the excess saturation. It is proposed that only part of this component (0 ⩽ αSk ⩽ 1) is contributed to the surface runoff qs:

qSj = αSk × Δz × (θPj − θS )

(A9) 9

Journal of Hydrology 583 (2020) 124572

M. Ben Yona, et al.

While the residual, qPREFj flows downward as preferential flow to recharge the aquifer relatively quickly with a large mass of water, typical for environment with large natural pores as karst.

qPREFj = (1 − αSk ) × Δz × (θPj − θS )

(A10)

The next modules, surface flow, vadose zone and groundwater modules consist of various combinations of linear reservoirs. The typical differential equation for a linear reservoir:

KQOUTj QINj hj dh = − ; h0 = A dt A K

(A11) 3

here changes in the height of the reservoir’s water level (h (m)) is a function of the inflow to the reservoir QIN (m day and K the storage coefficient with the dimension of time (day). The outflow can then be calculated with:

QOUTj =

−1

), the reservoir area A (m2)

Ahj (A12)

K

The surface flow module takes as input part of the daily pulse of excess saturation (qs Eq. (A7)) and transforms it into the streamflow by a simple linear reservoir operator (Eqs. (A9) and (A10)). The baseflow discharge is calculated as a set of two cascade linear reservoirs. The input of the vadose zone consists of the two fluxes from the topsoil (epikarst) layer: Darcy law qDARCY Eq. (A1) and the preferential flow qPREF Eq. (A10).

QIN = qDARCY + QPREF

(A13)

The output represents vadose zone is the input of groundwater modules and it is calculated as a simple linear reservoir operator (Eq. (A12)). The baseflow is calculated by additional simple linear reservoir operator (Eq. (A12)). Matan Ben Yona: Conducted the research as part of his master degree. Alon Rimmer: Supervised, mentored, and helped with the study design. Eylon Shamir: Mentored and guided the research, assisted with writing, reviewing and editing of the manuscript. Iggy Litaor: Acquisition of funding, management and coordination responsibility for the research activity, reviewing and editing of the manuscript.

9, 407–422. Latron, J., Soler, M., Liorens, P., Gallart, F., 2008. Spatial and temporal variability of the hydrological response in a small Mediterranean research catchment (Vallcebre, Eastern Pyrenees). Hydrol. Process 22, 775–787. Michelson, H., 1979. The geology and paleogeography of the Golan Heights. Ph.D. Thesis, Tel Aviv Univ., 163 pp. (in Hebrew, English abstr). Mor, D. 1986. The volcanism of the Golan Heights. Ph.D. thesis, Hebrew Univ., Jerusalem, 159 pp., and GSI report, GSI/5/86 (in Hebrew, English abstr.). Mualem, Y., Dagan, G., 1978. Hydraulic conductivity of soils: unified approach to the statistical models. Soil Sci. Soc. Am. J. 42 (3), 392–395. Peleg, N., Shamir, E., Georgakakos, K.P., Morin, E., 2015. A framework for assessing hydrological regime sensitivity to climate change in a convective rainfall environment: a case study of two medium-sized eastern Mediterranean catchments, Israel. Hydrol. Earth Syst. Sci. 19, 567–581. https://doi.org/10.5194/hess-19-567-2015. Premium Solver Platform, User guide. (2010). Frontline Systems, Inc., Notes downloaded from the site http://www.solver.com. Ries, F., Schmidt, S., Sauter, M., Lange, J., 2017. Controls on runoff generation along a steep climatic gradient in the Eastern Mediterranean. J. Hydrol. Reg. Stud. 9, 18–33. Rimmer, A., Givati, A., 2014. Hydrology, in: Lake Kinneret – Ecology and Management. Rimmer, A., Givati, A., Samuels, R., Alpert, P., 2011. Using ensembles of climate model to evaluate future water and solutes budgets in Lake Kinneret Israel. J. Hydrol. 410, 248–259. Rimmer, A., Hartmann, A., 2012. Simplified conceptual structures and analytical solutions for groundwater discharge using reservoir equations. Chapter 10 in InTech Open Access book, “Water resources management and modeling”, pp. 978–953. Rimmer, A., Hartmann, A., 2014. Optimal hydrograph separation filter to evaluate transport routines of hydrological models. J. Hydrol. 514. Rimmer, A., Salingar, Y., 2006. Modelling precipitation-streamflow processes in Karst basin: the case of the Jordan River sources, Israel. J. Hydrol. 331, 524–542. Sade, R., Rimmer, A., Litaor, M.I., Shamir, E., Furman, A., 2014. Snow surface energy and mass balance in a warm temperate climate mountain. J. Hydrol. 519, 848–862. Sade, R., Rimmer, A., Samuels, R., Salingar, Y., Danisyuk, M., Alpert, M., 2016. Water management in a complex hydrological basin—application of water evaluation and planning tool (WEAP) to the Lake Kinneret Watershed, Israel, in: Integrated water resources management: concept, research and implementation. pp. 35–57. Samuels, R., Rimmer, A., Alpert, P., 2009. Effect of extreme rainfall events on the water resources of the Jordan River. J. Hydrol. 375, 513–523. Samuels, R., Rimmer, A., Hartmann, H., Krichak, S., Alpert, P., 2010. Climate change impacts on Jordan River flow: downscaling application from a regional climate model. J. Hydrometeorol. 11, 860–879. Samuels, R., Hochman, A., Baharad, A., Givati, A., Levi, Y., Yosef, Y., Saaroni, H., Ziv, B., Harpaza, T., Alpert, P., 2017. Evaluation and projection of extreme precipitation indices in the Eastern Mediterranean based on CMIP5 multi-model ensemble. Int. J. Clim. 1–18. https://doi.org/10.1002/joc.5334. Schaefli, B., Gupta, H., 2007. Do Nash values have value? Hydrol. Process. 21, 2075–2080. Shamir, E., Rimmer, A., Georgakakos, K., 2016. The use of an orographic precipitation

References Auerbach, M., Shmida, A., 1993. Vegetation change along an altitudinal gradient on Mt Hermon, Israel – no evidence for discrete communities. J. Ecol. 81 (1), 25–33. Ben-Zvi, A., 1988. Enhancement from runoff from a small watershed by cloud seeding. J. Hydrol. 101, 291–303. Dafny, E., Burg, A., Gvirtzman, H., 2006. Deduction of groundwater flow regime in a basaltic aquifer using geochemical and isotopic data: the Golan Heights, Israel case study. J. Hydrol. 330, 506–524. https://doi.org/10.1016/j.jhydrol.2006.04.002. Dafny, E., Gvirtzman, H., Burg, A., Fleischerc, L., 2003. The hydrogeology of the Golan basalt aquifer, Israel. Isr. J. Earth Sci. 52, 139–153. Danin, A., 1988. Flora and vegetation of Israel and adjacent areas. Zoogeography Isr. 251–276. Dubrovský, M., Hayes, M., Duce, P., Trnka, M., Svoboda, M., Zara, P., 2014. Multi-GCM projections of future drought and climate variability indicators for the Mediterranean region. Reg. Environ. Chang. 14, 1907–1919. https://doi.org/10.1007/s10113-0130562-z. Eckhardt, K., 2005. How to construct recursive digital filters for baseflow separation. Hydrol. Process 19, 507–515. Frumkin, A., Shimron, A.E., Miron, Y., 1998. Karst morphology across a steep climatic gradient, Southern Mount Hermon, Israel. Zeitschrift für Geomorphologie 109, 23–40. Gallart, F., Llorens, P., Latron, J., Regüés, D., 2002. Hydrological processes and their seasonal controls in a small Mediterranean mountain catchment in the Pyrenees. Hydrol. Earth Syst. 6, 527–537. Gilad, D., 1980. The Hermon as the main water source of Israel. Mt. Hermon nature and landscape (in hebrew). A. Shmida and M. Livne. Tel-Aviv, Kibbutz. Gilad, D., Schwartz, S., 1978. Hydrogeology of the Jordan sources aquifers. Isr. Hydrol Serv. Rep. Hydro/5/78, 58 (in Hebrew). Gilad, D., Bonne, J., 1990. Snowmelt of Mt. Hermon and its contribution to the sources of the Jordan River. J. Hydrol. 114 (1/2), 1–15. Gilboa, Y., Gal, G., Markel, D., Rimmer, A., Evans, B.M., Friedler, E., 2015. Effect of landuse change scenarios on nutrients and TSS loads. Environ. Process. 2 (4), 593–607. https://doi.org/10.1007/s40710-015-0109-z. Hartmann, A., Wagener, T., Rimmer, A., Lange, J., Brielmann, H., Weiler, M., 2013. Testing the realism of model structures to identify karst system processes using water quality and quantity signatures. Water Resour. Res. 49, 3345–3358. Heimann, A. 1990. The Dead Sea Rift and its margins development in Northern Israel in the Pliocene and Pleistocene. Ph.D. thesis, Hebrew Univ., Jerusalem, 83 pp., and GSI report GSI/28/90 (in Hebrew, English abstr.). Hogue, T.S., Sorooshian, S., Gupta, H., Holz, A., Braatz, D., 2000. A Multistep automatic calibration scheme for river forecasting models. J. Hydrometeorol. 1, 524–542. Kutiel, H., 1987. Rainfall variations in the Galilee (Israel). I. Variations in the spatial distribution in the periods 1931–1960 and 1951–1980. J. Hydrol. 94, 331–344. Lavee, H., Imerson, A.C., Sarah, P., 1998. The impact of climate change on geomorphology and desertification along the Mediterranean arid transect. Land Degrad. Dev.

10

Journal of Hydrology 583 (2020) 124572

M. Ben Yona, et al.

1007/978-3-642-66049-8_6. Shmida, A., 1977. A quantities analysis of the tragacnthic vegetation of Mt. Hermon and its relation to environmental factors, Hebrew University of Jerusalem. Ph.D., 291. Shtober-Zisu, N., Inbar, M., 2005. Downstream effects of river impoundment on hydrological and geomorphological aspects of bedrock rivers (Golan Heights, Israel). IAHS Sp. Pub. 299, 219–224. Ziv, B., Shilo, E., Rimmer, A., 2014. Meteorology, in: Lake Kinneret – ecology and management.

model to assess the precipitation spatial distribution in Lake Kinneret watershed. Water 8, 591. Shentsis, I., Inbar, N., Rosenthal, E., Magri, F., 2018. Numerical representation of rainfall field in basins of the Upper Jordan River and of the Yarmouk River. Environ. Earth Sci. 77, 798. https://doi.org/10.1007/s12665-018-7976-3. Singer, A., 2007. The soils of Israel. Teatre Lliure. https://doi.org/10.1525/auk.2013. 130.4.812. Singer, A., 1987. Land evaluation of basaltic terrain under semi-arid to mediterranean conditions in the Golan heights. Soil Use Manag. 3, 155–162. https://doi.org/10.

11