# Modeling the Effects of Land-Cover Change on Rainfall-Runoff Relationships in a Semiarid, Eastern Mediterranean Watershed

Authors: Noa Ohana-Levi; Arnon Karnieli; Roey Egozi; Amir Givati; Aviva Peeters
Journal: Advances in Meteorology
 (2015)

---

## Abstract

Temporal changes and spatial patterns are often studied by analyzing land-cover changes (LCCs) using spaceborne images. LCC is an important factor, affecting runoff within watersheds. The objective was to estimate the effects of 20 years of LCCs on rainfall-runoff relations in an extreme rainfall event. A 1989 Landsat TM-derived classification map was used as input for a Kinematic Runoff and Erosion (KINEROS) hydrological model along with the precipitation data of an extreme rainfall event. Model calibration was performed using measured runoff volume data. Validation of the model performance was conducted by comparing the model results to measured data. A similar procedure was used with a 2009 land-cover classification map as an input to the KINEROS model, along with similar precipitation data and calibration parameters, in order to understand the possible outcomes of a rainfall event of such a magnitude and duration after 20 years of LCCs. The results show an increase in runoff volume and peak discharge between the time periods as a result of LCCs. A strong relationship was detected between vegetation cover and the runoff volume. The LCCs with most pronounced effects on runoff volumes were related to urbanization and vegetation removal.


---

## Body

## 1. Introduction

Land-cover mapping and monitoring can serve as important indicators of landscape and environmental status, distributions, and patterns. A common and effective way to better understand temporal changes, as well as the spatial variability of an area, is through the study of land-cover change (LCC) and its analysis [1–3]. LCCs are caused by human activities and in many cases have local and regional impacts on populations [4]. There is a growing awareness of LCCs importance and their impacts in order to gain a better understanding of their effects on hydrologic, ecologic, climatic, and biologic processes over space and time [5–8]. One of these processes is the effect of LCC on runoff regime within a watershed. The LCCs are often related to urbanization, a shift to agriculture, afforestation, and vegetation removal [1, 9–12].The impact of LCCs on the landscape features of watersheds greatly affects slopes and channel flows [13]. One of the most influential changes is urbanization, as urban areas increase runoff rates and peak discharges due to increased imperviousness and reduced infiltration of precipitation [14–17]. Moreover, vegetation removal, which usually accompanies the process of urbanization, is responsible for the reduction of rainfall interception and storage and, in arid/semiarid regions, the creation of physical crusts [18]. Other LCCs can be attributed to the use of land resources for agricultural activities, leading to the replacement of grasslands and forested areas with croplands [9]. The spatial distribution and patterns of land-cover dynamics are crucial factors in estimating imperviousness along the slopes and quantifying the areas that contribute runoff to the channels that drain the watershed area [19]. Understanding the spatial variability and patterns in LCCs can greatly improve the process of determining the most influential LCCs and other environmental factors on rainfall-runoff relationship [20].Expected worldwide climate-change trends are likely to introduce shifts in the frequency and intensity of extreme climatic events, including increase in mean annual temperatures and increase/decrease in mean annual rainfall, depending on the region in question [21, 22], and will have a strong impact on hydrologic regime in the watershed scale. Future global trends predict an increase in heavy precipitation events [21]. This projection and its probable effects on the hydrologic regime in the watershed scale led to the decision to focus on one storm event, with a return period of 50 years. The problem at issue is dealing with the extensive urbanization processes occurring worldwide, which are most likely to increase runoff and consequently endanger the urban areas along stream channels during future extreme storm events [23]. These expected climate-change trends highlight the importance of quantifying the shifts in land-cover through time since the last extreme storm, in order to predict the outcome of a similar hazardous event in the future and the effects it may have on the well-being of the population within the watershed.Several means for observing LCCs have been developed and documented. Satellite images have a great potential for monitoring and analyzing LCCs both temporally and spatially [24]. In this study, remote-sensing data was integrated with geographic information system (GIS) techniques. In the scope of this work, LCCs as well as patterns and distribution were monitored and analyzed using remote sensing and GIS techniques and then implemented into a hydrological model in order to estimate their effects on the hydrologic patterns caused by a sever rainstorm in a watershed scale [17, 25].The hydrological model that was chosen for this study was KINEROS2, for two main reasons: (1) it was originally designed for arid and semiarid regions and is appropriate for the research area; (2) it is an even-based model and suits the study scope of modelling a single hydrological event [26, 27]. Moreover, this research implements statistical and spatial analyses for studying spatial and temporal LCCs and hydrological trends. In order to quantify temporal LCCs, a statistical ordination approach was used, after which the LCC effects on runoff volume were analyzed using GIS spatial analysis tools. Land-cover patterns effects of runoff volume were also analyzed spatially, using spatial statistics GIS tools. The aim of this work was to estimate the temporal effects of LCCs on rainfall-runoff relations within the Yarkon-Ayalon watershed, Israel, on a subbasin scale. The specific objectives were twofold: (1) to extract and characterize land-cover classes, for each subbasin, using satellite imagery for assessing LCCs between two time periods, 20 years apart; (2) to use the LCCs findings, together with the data of an extreme rainfall event scenario, in order to simulate runoff variables for past and present land-cover status.


## 2. Methodology

### 2.1. Study Area

The research was conducted within the Yarkon-Ayalon watershed, located in central Israel (Figure1(a)), due to its hydrologic vulnerability. This watershed contains the most densely populated region of the country: Tel Aviv and some of its satellite cities [28]. The watershed is divided by the border between Israel and the West Bank, Palestinian Authority (PA), hence the limitations on conducting extensive field research and the inevitable need for a remote sensing tool and LCC monitoring methods. This area has experienced an extensive urbanization process during the past decades, at the expense of agricultural and natural areas [14, 29].Figure 1
Yarkon-Ayalon watershed: (a) watershed location within Israel and the West Bank (WB), Palestinian Authority (PA); (b) watershed elevation map.
(a)
(b)The Yarkon-Ayalon watershed has a total drainage area of 1,805 km2, originating in the Samaria Mountains, located in the PA in the east and draining west towards the Mediterranean Sea. There is a large variability in elevation along this watershed, with its headwaters defined as mountainous and reaching altitudes of about 800–900 meters above sea level (Figure 1(b)) [30, 31]. The watershed area is located in a semiarid climate zone with a mean annual precipitation of about 550 mm and a mean annual discharge of 43 million cubic meters (MCM) [23].The Yarkon-Ayalon watershed consists of several soil types, as a factor of slopes and geological characteristics. Soft rocks (chalk and marl) and harder rocks (limestone and dolomite) are both abundant. The soils in this region are mostly terra rossa and brown rendzina along the eastern region of the watershed, while the plains along the central area are characterized by gromosols [32]. Most of the vegetation in the watershed area is composed of Mediterranean forests, woodlands, and shrubs, such as Palestine oaks, mastic shrubs, carob trees, terebinths,Poteriums.The research area contains a wide array of land-cover and land-use types, some of which are still natural while most are controlled and managed. These include forests, different types of residential land-uses (high/low density, built/unbuilt areas, large cities, and villages), industrial areas, commercial areas, roads, grasslands, agricultural lands (row crops, orchards), and quarries. All of these differ in magnitude, size, quality, and shape between Israel and the PA.The research was conducted on a subbasin scale. Therefore, the Yarkon-Ayalon watershed was partitioned into six subbasins according to the locations of the Israeli Hydrological Service hydrometric stations: Yarkon, Beit Dagan, Shilo, Natuf, Kana, and Lod (Figure2). Each subbasin was chosen according to hydrometric data availability and was studied independently. As demonstrated in Figure 2, the Yarkon subbasin contains both the Kana and Shilo subbasins, and the Beit Dagan subbasin comprises the Natuf and Lod subbasins. In 1957 a reservoir named Mishmar Ayalon was built at the eastern part of Lod subbasin, in order to minimize the peak discharges during severe flooding events, and was never breached. This means that the outflow volume runoff and discharge values that were measured at the hydrometric station represented only the runoff generated downstream from the reservoir. Therefore, the Mishmar Ayalon drainage area was excluded from the Lod and Beit Dagan subbasin areas (Figure 2).Figure 2
Yarkon-Ayalon watershed partitioned into six subbasins according to hydrometric stations. The number within each subbasin polygon represents its area.


### 2.2. Hydrological Model and Tools

The model used in the course of this research for simulating runoff within the Yarkon-Ayalon watershed was KINEROS2 [26, 27]. This model is designed for assessing runoff rates and volumes in ungauged watersheds in arid and semiarid regions. It is a physically based, distributed, event-oriented, and rainfall-runoff model, which simulates runoff response over the watershed. The watershed is divided into two types of spatially distributed model elements with a specified connectivity: (1) overland flow elements or planes (polygons) and (2) channel elements (lines). These elements can be oriented so that a one-dimensional flow can be assumed. Surface flow is simulated for all planes and channels using a finite difference solution of the one-dimensional kinematic wave equations [26, 33–35].KINEROS2 describes a Hortonian overland flow from complex watersheds and accounts for the spatial and temporal variability of rainfall input. It deals with processes of rainfall, interception, infiltration, overland flow, open channel flow, erosion, sediment transport, reservoir routing, and sedimentation [33, 34, 36]. The current research focuses on the rainfall-runoff process only. A detailed and thorough description of the model, its processes, and mechanisms can be found in Woolhiser et al. (1990) and Semmens et al. (2008).The choice of KINEROS2 was based on its capability to simulate runoff in semiarid regions, which is appropriate for the climate regime along the Yarkon-Ayalon watershed. Moreover, this study aimed to examine the spatial behavior and characteristics of one storm event. KINEROS2, being an event-based model that deals with specific storms, was suitable for the purpose of this study.The model is operated using a GIS interface named the Automated Geospatial Watershed Assessment (AGWA2) tool [27, 37, 38]. It was developed in order to implement two models: KINEROS2 and the Soil and Water Assessment Tool (SWAT) via a GIS interface. The data requirements for running the AGWA2 tool include elevation, classified land-cover map, soil map, and precipitation data. The model input variables are derived from these data using look-up tables (LUTs), provided by the tool. The model requires first generating the watershed outline followed by dividing it into model elements. The desired model is then chosen (KINEROS2 in this case), after which the soils and land-cover maps, along with LUT data, are entered into the system, and the storm event precipitation data is generated in order to complete the input dataset. The model can then be calibrated and run, and the simulation results can be displayed. The model output includes runoff, sediment yield, infiltration, peak runoff rate, and peak sediment discharge [27, 38]. The current study focused only on runoff simulations since gauged validation data for the other output variables were not available.


### 2.3. Model Input Data

In order to meet the model’s requirements, the Yarkon-Ayalon watershed was initially divided into six different subbasins according to available hydrometric data for validation of the predicted results (Figure2). Each of these subbasins was modeled separately, and model results were analyzed. Figure 3 shows the model’s flowchart that is specified hereinafter.Figure 3
The study flowchart. AGWA2 is the Automated Geospatial Watershed Assessment tool. VOL is the objective function used for calibration that calculates the total runoff volume ratio between predicted and observed measurements.A digital elevation model (DEM) was essential for providing elevation data. The DEM for this research was produced from optical stereo data acquired by the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) that operates on the NASA Terra platform and has a spatial resolution of 30 m [39].Soil data was represented through a digital soil map of the world provided by the Food and Agriculture Organization (FAO) of the United Nations (UN) and the UN Educational, Scientific and Cultural Organization (UNESCO). It was last updated in 2003 and has a spatial resolution of 1 : 5 million [40].The model used a specific extreme rainfall storm that occurred between December 29, 1991, and January 3, 1992, and had an average rainfall depth of 165 mm over the study watershed. Raw precipitation data from 24 stations, within and in vicinity of the watershed, were extracted from the Israeli Meteorological Service. The Thiessen polygon interpolation method [41] was conducted in order to receive estimations of precipitation values for each subbasin. Based on a Generalized Pareto Distribution method [42], this storm event is known to have had a precipitation return period of 50 years and led to the loss of 12 lives nationwide, damaged agricultural yields, and direct damage to buildings and infrastructure. Several main cities along the shoreline in central Israel were flooded, and the storm events during this winter were officially classified as a “natural disaster” [23]. Within the Yarkon-Ayalon watershed, some 153 million cubic meters (MCM) was discharged (compared to a mean annual discharge of 43 MCM) with a maximum discharge of 491 m3/s [23, 43]. Finally, a saturation index between 0.14 and 0.93 was entered into the model, representing the initial relative soil saturation [44].Land-cover data was extracted using satellite imagery. In order to understand the temporal effect of LCCs on rainfall-runoff relationships, a model simulation of the same storm event, given the land-cover data status of a time period 20 years later, was implemented. This required using two images, both from the Landsat Thematic Mapper (TM), with a 30 m spatial resolution; the earlier image, representing the land-cover status during the time of the selected storm event, was acquired on December 17, 1989, and the second image was acquired about 19 years later, on January 3, 2009. The two images underwent radiometric and atmospheric corrections, as well as a topographic correction, due to the steep topography of the research area. All of these processes were carried out using Atmospheric/Topographic Correction (ATCOR) software that is specifically detailed in Richter (1998, 2010) [45, 46].Six different land-cover classes were then decided upon in order to examine LCC through the years: (1) residential areas; (2) agricultural land; (3) orchards; (4) natural areas; (5) forests; and (6) bare exposed rocks, mines, and quarries. In order to extract land-cover patterns and distribution, a supervised maximum likelihood classification technique was conducted for both images [47], followed by an accuracy assessment procedure against high resolution orthophoto data, using 600 validation points. The methods chosen for classification accuracy assessment and evaluation were Kappa coefficient and overall accuracy as calculated from an error (confusion) matrix [48, 49].Since each of the land-cover classes has its own characteristics, they were expressed and linked to the land-cover classes through the LUT; each of the six land-cover classes received representative values for the following variables: estimated canopy cover, interception, Manning’s roughness coefficient, and imperviousness (percent paved area), specifically detailed by Miller et al. (2007) and Semmens et al. (2008). The LUTs are provided by the AGWA2 tool and can be modified by the user.


### 2.4. Model Calibration, Validation, and Output

The AGWA2 tool provides a number of output variables. In this work, the focus was on outflow runoff volume (m3) and maximum peak discharge (m3/s) (Figure 3). For each subbasin, these model output variables were simulated and examined. Thein situ measured data for these variables enabled validation of the model results. The output map given by the model is discretized and represented by a collection of plane elements (polygons) and channel elements (polylines). Each element has its own value of the variable being examined.The evaluation process was carried out using objective functions (OFs) as measures of the ratios and differences between the simulated values derived by the model, using the 1989 classification image as an input representing land-cover status during the storm, and the observed data measured during the actual rainstorm. In order to avoid bias, a number of OFs, rather than a single one, were chosen to evaluate the model accuracy [50, 51]. The calibration and validation OFs are summarized in Table 1. The objective was to keep the deviation between observed and simulated volume runoff values (VOL) as low as possible and in all cases under 1%, so that there was a minimum of 99% accuracy between measured and predicted total runoff volumes.Table 1
Calibration and validation objective functions (OFs). Vobs and Vsim are the observed and simulated runoff volume values, respectively; QMXobs and QMXsim are the observed and simulated peak discharge values, respectively; Qpred and Qsim are the simulated and observed discharge values, respectively.i denotes time increment.
Function
Description
VOL
=
Vobs
∗
100
Vsim
-
100
Total runoff volume ratio between simulated and observed measurements, VOL, in %.
PQR
=
QMXobs
*
100
QMXsim
-
100
Peak discharge ratio between simulated and observed measurements, PQR, in %.
ADPQ
=
QMXpred
-
QMXsim
Absolute deviation between peak discharges of  simulated and observed measurements, ADPQ, in m3/s.
MAQD
=
1
n
∑
Qpred
i
-
Qsim
i
Mean absolute discharge deviation between simulated and observed measurements, MAQD, in m3/s.
RMS
=
1
n
∑
Qpred
i
-
Qsim
i
2
Root mean square discharge deviation between simulated and observed measurements, RMS, in m3/s.
MXAQ
=
max
⁡
Qpred
i
-
Qsim
i
Maximum absolute discharge deviation between simulated and observed measurements, MXAQ, in m3/s.The calibration process included changing the model’s parameter multipliers, which are values used to adjust plane and channel parameters (e.g., percent cover, interception, roughness for planes, and width, depth, and field effective saturated hydraulic conductivity) [27]. This was conducted up to the point where VOL is under 1%.All of the OFs were calculated separately for the six subbasins. Each of the validation OF results was plotted against the subbasins’ sizes (in km2). A regression model, along with its significance level, was fitted in order to determine the relationship between the model accuracy and the subbasin size.Initially, the model was simulated using the 1989 classification map, representing land-cover status during the rainfall event of winter 1991-1992. In order to determine the changes in rainfall-runoff relations caused by LCCs during a period of 19 years, the 2009 land-cover map was entered into the model instead of the 1989 map, while all other input variables remained static (e.g., calibration parameters, soil saturation index, precipitation data, soil map, and DEM). This was performed in order to examine a scenario in which a rainfall event, with a fifty-year return period, occurs after 19 years of LCCs and to evaluate any changes and trends in future runoff volumes and peak discharges.


### 2.5. Analyses

Two major trends were analyzed from the model’s simulated results. The first was a temporal trend, determining whether LCCs throughout 19 years affected runoff volumes and peak discharges in each of the six Yarkon-Ayalon subbasins and to what extent. The second was a spatial trend evaluating the factors that best-explain the spatial variability of the runoff volume across the different subbasins.In order to conduct the temporal analysis, two factors were examined: (1) LCC between 1989 and 2009, using the land-cover maps; and (2) comparison of model results of the two time periods and how they were affected by LCCs. The first factor was analyzed by summing the total number of pixels for each land-cover in each subbasin, for both land-cover maps. The total number of pixels for each of the six land-cover classes was compared between the two images (1989 versus 2009) in order to analyze shifts. Then, the two land-cover maps were overlaid and a transition map was extracted, followed by summing up all of the transitions for every subbasin. For example, if a pixel with an agriculture land-cover type in the 1989 map was classified as a residential type in the 2009 map, it was defined as an agriculture-to-residential transition. For every subbasin, all of the pixels having this transition type were summed. This led to a total of 31 different transition phases (with “no-change” transition being the same for all types of land-cover undergoing no transition). The different transition phases for all subbasins were converted into a transition matrix and normalized to the area of each subbasin. Then, the transitions were analyzed using principal component analysis (PCA) ordination in order to determine which transitions were more pronounced in each subbasin.The second factor of the temporal analysis initially involved conducting, for every subbasin, a comparison between total simulated runoff volume and peak discharge values, between the real rainfall event and the 2009 rainfall event scenario, and between their land-cover representations (1989 and 2009) in order to examine the runoff trends. This step included spatially examining the LCCs. It was performed by zooming into every plane element (polygon) in the subbasin and exploring the changes within the subbasin’s elements, by calculating the percentage change between the two simulated runoff volume values (1989 and 2009) of every plane element. Only the most pronounced changes in runoff volumes were of interest, so only plane elements that had a value of above or below 1 standard deviation (SD) percent change were considered. The third step was to determine how LCCs affected the runoff changes in these plane elements. To do so, a zonal statistic tool for majority calculation of land-cover classes was used [52] for each plane element that was extracted in step 2. This means that each plane element (lower or higher than 1 SD of percent change) received one value of the land-cover class that had the largest representation within it. This was done for both time period images so that a qualitative comparison could be performed for each plane element that exhibited runoff changes satisfying the threshold of ±1 SD.The spatial analysis process was performed in order to gain a better understanding regarding the environmental factors spatially affecting the runoff trends and was conducted only for the 2009 land-cover map. For every subbasin, the simulated runoff volume values were considered as the dependent variable. Three explanatory variables were examined as predictors of the runoff volume: land-cover, represented by the land-cover classification map, elevation, represented by DEM, and vegetation cover, represented by a normalized difference vegetation index (NDVI) map [53], extracted from the 2009 image. In order to determine the relationship characteristics between these explanatory variables and the runoff volume, a geographically weighted regression (GWR) was conducted. This method is designed to consider spatial data and nonstationarity, where the variable estimates vary locally, and spatially varying relationships can be explored [54, 55]. The results were expected to point out the relations between the runoff volume and the three explanatory variables, represented by regression coefficient R
2, for each subbasin. The main goal was to understand which variable, or combination of variables, affects runoff volume variability the most for each subbasin.Seven GWR tests were built and considered for each subbasin. They were conducted for the runoff volume variable against each independent variable separately, followed by modeling paired independent variables against runoff volume and finally by running all three independent variables. The model with the highestR
2 for each subbasin implies which variable, or combination of variables, had the strongest effect on runoff.


## 2.1. Study Area

The research was conducted within the Yarkon-Ayalon watershed, located in central Israel (Figure1(a)), due to its hydrologic vulnerability. This watershed contains the most densely populated region of the country: Tel Aviv and some of its satellite cities [28]. The watershed is divided by the border between Israel and the West Bank, Palestinian Authority (PA), hence the limitations on conducting extensive field research and the inevitable need for a remote sensing tool and LCC monitoring methods. This area has experienced an extensive urbanization process during the past decades, at the expense of agricultural and natural areas [14, 29].Figure 1
Yarkon-Ayalon watershed: (a) watershed location within Israel and the West Bank (WB), Palestinian Authority (PA); (b) watershed elevation map.
(a)
(b)The Yarkon-Ayalon watershed has a total drainage area of 1,805 km2, originating in the Samaria Mountains, located in the PA in the east and draining west towards the Mediterranean Sea. There is a large variability in elevation along this watershed, with its headwaters defined as mountainous and reaching altitudes of about 800–900 meters above sea level (Figure 1(b)) [30, 31]. The watershed area is located in a semiarid climate zone with a mean annual precipitation of about 550 mm and a mean annual discharge of 43 million cubic meters (MCM) [23].The Yarkon-Ayalon watershed consists of several soil types, as a factor of slopes and geological characteristics. Soft rocks (chalk and marl) and harder rocks (limestone and dolomite) are both abundant. The soils in this region are mostly terra rossa and brown rendzina along the eastern region of the watershed, while the plains along the central area are characterized by gromosols [32]. Most of the vegetation in the watershed area is composed of Mediterranean forests, woodlands, and shrubs, such as Palestine oaks, mastic shrubs, carob trees, terebinths,Poteriums.The research area contains a wide array of land-cover and land-use types, some of which are still natural while most are controlled and managed. These include forests, different types of residential land-uses (high/low density, built/unbuilt areas, large cities, and villages), industrial areas, commercial areas, roads, grasslands, agricultural lands (row crops, orchards), and quarries. All of these differ in magnitude, size, quality, and shape between Israel and the PA.The research was conducted on a subbasin scale. Therefore, the Yarkon-Ayalon watershed was partitioned into six subbasins according to the locations of the Israeli Hydrological Service hydrometric stations: Yarkon, Beit Dagan, Shilo, Natuf, Kana, and Lod (Figure2). Each subbasin was chosen according to hydrometric data availability and was studied independently. As demonstrated in Figure 2, the Yarkon subbasin contains both the Kana and Shilo subbasins, and the Beit Dagan subbasin comprises the Natuf and Lod subbasins. In 1957 a reservoir named Mishmar Ayalon was built at the eastern part of Lod subbasin, in order to minimize the peak discharges during severe flooding events, and was never breached. This means that the outflow volume runoff and discharge values that were measured at the hydrometric station represented only the runoff generated downstream from the reservoir. Therefore, the Mishmar Ayalon drainage area was excluded from the Lod and Beit Dagan subbasin areas (Figure 2).Figure 2
Yarkon-Ayalon watershed partitioned into six subbasins according to hydrometric stations. The number within each subbasin polygon represents its area.


## 2.2. Hydrological Model and Tools

The model used in the course of this research for simulating runoff within the Yarkon-Ayalon watershed was KINEROS2 [26, 27]. This model is designed for assessing runoff rates and volumes in ungauged watersheds in arid and semiarid regions. It is a physically based, distributed, event-oriented, and rainfall-runoff model, which simulates runoff response over the watershed. The watershed is divided into two types of spatially distributed model elements with a specified connectivity: (1) overland flow elements or planes (polygons) and (2) channel elements (lines). These elements can be oriented so that a one-dimensional flow can be assumed. Surface flow is simulated for all planes and channels using a finite difference solution of the one-dimensional kinematic wave equations [26, 33–35].KINEROS2 describes a Hortonian overland flow from complex watersheds and accounts for the spatial and temporal variability of rainfall input. It deals with processes of rainfall, interception, infiltration, overland flow, open channel flow, erosion, sediment transport, reservoir routing, and sedimentation [33, 34, 36]. The current research focuses on the rainfall-runoff process only. A detailed and thorough description of the model, its processes, and mechanisms can be found in Woolhiser et al. (1990) and Semmens et al. (2008).The choice of KINEROS2 was based on its capability to simulate runoff in semiarid regions, which is appropriate for the climate regime along the Yarkon-Ayalon watershed. Moreover, this study aimed to examine the spatial behavior and characteristics of one storm event. KINEROS2, being an event-based model that deals with specific storms, was suitable for the purpose of this study.The model is operated using a GIS interface named the Automated Geospatial Watershed Assessment (AGWA2) tool [27, 37, 38]. It was developed in order to implement two models: KINEROS2 and the Soil and Water Assessment Tool (SWAT) via a GIS interface. The data requirements for running the AGWA2 tool include elevation, classified land-cover map, soil map, and precipitation data. The model input variables are derived from these data using look-up tables (LUTs), provided by the tool. The model requires first generating the watershed outline followed by dividing it into model elements. The desired model is then chosen (KINEROS2 in this case), after which the soils and land-cover maps, along with LUT data, are entered into the system, and the storm event precipitation data is generated in order to complete the input dataset. The model can then be calibrated and run, and the simulation results can be displayed. The model output includes runoff, sediment yield, infiltration, peak runoff rate, and peak sediment discharge [27, 38]. The current study focused only on runoff simulations since gauged validation data for the other output variables were not available.


## 2.3. Model Input Data

In order to meet the model’s requirements, the Yarkon-Ayalon watershed was initially divided into six different subbasins according to available hydrometric data for validation of the predicted results (Figure2). Each of these subbasins was modeled separately, and model results were analyzed. Figure 3 shows the model’s flowchart that is specified hereinafter.Figure 3
The study flowchart. AGWA2 is the Automated Geospatial Watershed Assessment tool. VOL is the objective function used for calibration that calculates the total runoff volume ratio between predicted and observed measurements.A digital elevation model (DEM) was essential for providing elevation data. The DEM for this research was produced from optical stereo data acquired by the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) that operates on the NASA Terra platform and has a spatial resolution of 30 m [39].Soil data was represented through a digital soil map of the world provided by the Food and Agriculture Organization (FAO) of the United Nations (UN) and the UN Educational, Scientific and Cultural Organization (UNESCO). It was last updated in 2003 and has a spatial resolution of 1 : 5 million [40].The model used a specific extreme rainfall storm that occurred between December 29, 1991, and January 3, 1992, and had an average rainfall depth of 165 mm over the study watershed. Raw precipitation data from 24 stations, within and in vicinity of the watershed, were extracted from the Israeli Meteorological Service. The Thiessen polygon interpolation method [41] was conducted in order to receive estimations of precipitation values for each subbasin. Based on a Generalized Pareto Distribution method [42], this storm event is known to have had a precipitation return period of 50 years and led to the loss of 12 lives nationwide, damaged agricultural yields, and direct damage to buildings and infrastructure. Several main cities along the shoreline in central Israel were flooded, and the storm events during this winter were officially classified as a “natural disaster” [23]. Within the Yarkon-Ayalon watershed, some 153 million cubic meters (MCM) was discharged (compared to a mean annual discharge of 43 MCM) with a maximum discharge of 491 m3/s [23, 43]. Finally, a saturation index between 0.14 and 0.93 was entered into the model, representing the initial relative soil saturation [44].Land-cover data was extracted using satellite imagery. In order to understand the temporal effect of LCCs on rainfall-runoff relationships, a model simulation of the same storm event, given the land-cover data status of a time period 20 years later, was implemented. This required using two images, both from the Landsat Thematic Mapper (TM), with a 30 m spatial resolution; the earlier image, representing the land-cover status during the time of the selected storm event, was acquired on December 17, 1989, and the second image was acquired about 19 years later, on January 3, 2009. The two images underwent radiometric and atmospheric corrections, as well as a topographic correction, due to the steep topography of the research area. All of these processes were carried out using Atmospheric/Topographic Correction (ATCOR) software that is specifically detailed in Richter (1998, 2010) [45, 46].Six different land-cover classes were then decided upon in order to examine LCC through the years: (1) residential areas; (2) agricultural land; (3) orchards; (4) natural areas; (5) forests; and (6) bare exposed rocks, mines, and quarries. In order to extract land-cover patterns and distribution, a supervised maximum likelihood classification technique was conducted for both images [47], followed by an accuracy assessment procedure against high resolution orthophoto data, using 600 validation points. The methods chosen for classification accuracy assessment and evaluation were Kappa coefficient and overall accuracy as calculated from an error (confusion) matrix [48, 49].Since each of the land-cover classes has its own characteristics, they were expressed and linked to the land-cover classes through the LUT; each of the six land-cover classes received representative values for the following variables: estimated canopy cover, interception, Manning’s roughness coefficient, and imperviousness (percent paved area), specifically detailed by Miller et al. (2007) and Semmens et al. (2008). The LUTs are provided by the AGWA2 tool and can be modified by the user.


## 2.4. Model Calibration, Validation, and Output

The AGWA2 tool provides a number of output variables. In this work, the focus was on outflow runoff volume (m3) and maximum peak discharge (m3/s) (Figure 3). For each subbasin, these model output variables were simulated and examined. Thein situ measured data for these variables enabled validation of the model results. The output map given by the model is discretized and represented by a collection of plane elements (polygons) and channel elements (polylines). Each element has its own value of the variable being examined.The evaluation process was carried out using objective functions (OFs) as measures of the ratios and differences between the simulated values derived by the model, using the 1989 classification image as an input representing land-cover status during the storm, and the observed data measured during the actual rainstorm. In order to avoid bias, a number of OFs, rather than a single one, were chosen to evaluate the model accuracy [50, 51]. The calibration and validation OFs are summarized in Table 1. The objective was to keep the deviation between observed and simulated volume runoff values (VOL) as low as possible and in all cases under 1%, so that there was a minimum of 99% accuracy between measured and predicted total runoff volumes.Table 1
Calibration and validation objective functions (OFs). Vobs and Vsim are the observed and simulated runoff volume values, respectively; QMXobs and QMXsim are the observed and simulated peak discharge values, respectively; Qpred and Qsim are the simulated and observed discharge values, respectively.i denotes time increment.
Function
Description
VOL
=
Vobs
∗
100
Vsim
-
100
Total runoff volume ratio between simulated and observed measurements, VOL, in %.
PQR
=
QMXobs
*
100
QMXsim
-
100
Peak discharge ratio between simulated and observed measurements, PQR, in %.
ADPQ
=
QMXpred
-
QMXsim
Absolute deviation between peak discharges of  simulated and observed measurements, ADPQ, in m3/s.
MAQD
=
1
n
∑
Qpred
i
-
Qsim
i
Mean absolute discharge deviation between simulated and observed measurements, MAQD, in m3/s.
RMS
=
1
n
∑
Qpred
i
-
Qsim
i
2
Root mean square discharge deviation between simulated and observed measurements, RMS, in m3/s.
MXAQ
=
max
⁡
Qpred
i
-
Qsim
i
Maximum absolute discharge deviation between simulated and observed measurements, MXAQ, in m3/s.The calibration process included changing the model’s parameter multipliers, which are values used to adjust plane and channel parameters (e.g., percent cover, interception, roughness for planes, and width, depth, and field effective saturated hydraulic conductivity) [27]. This was conducted up to the point where VOL is under 1%.All of the OFs were calculated separately for the six subbasins. Each of the validation OF results was plotted against the subbasins’ sizes (in km2). A regression model, along with its significance level, was fitted in order to determine the relationship between the model accuracy and the subbasin size.Initially, the model was simulated using the 1989 classification map, representing land-cover status during the rainfall event of winter 1991-1992. In order to determine the changes in rainfall-runoff relations caused by LCCs during a period of 19 years, the 2009 land-cover map was entered into the model instead of the 1989 map, while all other input variables remained static (e.g., calibration parameters, soil saturation index, precipitation data, soil map, and DEM). This was performed in order to examine a scenario in which a rainfall event, with a fifty-year return period, occurs after 19 years of LCCs and to evaluate any changes and trends in future runoff volumes and peak discharges.


## 2.5. Analyses

Two major trends were analyzed from the model’s simulated results. The first was a temporal trend, determining whether LCCs throughout 19 years affected runoff volumes and peak discharges in each of the six Yarkon-Ayalon subbasins and to what extent. The second was a spatial trend evaluating the factors that best-explain the spatial variability of the runoff volume across the different subbasins.In order to conduct the temporal analysis, two factors were examined: (1) LCC between 1989 and 2009, using the land-cover maps; and (2) comparison of model results of the two time periods and how they were affected by LCCs. The first factor was analyzed by summing the total number of pixels for each land-cover in each subbasin, for both land-cover maps. The total number of pixels for each of the six land-cover classes was compared between the two images (1989 versus 2009) in order to analyze shifts. Then, the two land-cover maps were overlaid and a transition map was extracted, followed by summing up all of the transitions for every subbasin. For example, if a pixel with an agriculture land-cover type in the 1989 map was classified as a residential type in the 2009 map, it was defined as an agriculture-to-residential transition. For every subbasin, all of the pixels having this transition type were summed. This led to a total of 31 different transition phases (with “no-change” transition being the same for all types of land-cover undergoing no transition). The different transition phases for all subbasins were converted into a transition matrix and normalized to the area of each subbasin. Then, the transitions were analyzed using principal component analysis (PCA) ordination in order to determine which transitions were more pronounced in each subbasin.The second factor of the temporal analysis initially involved conducting, for every subbasin, a comparison between total simulated runoff volume and peak discharge values, between the real rainfall event and the 2009 rainfall event scenario, and between their land-cover representations (1989 and 2009) in order to examine the runoff trends. This step included spatially examining the LCCs. It was performed by zooming into every plane element (polygon) in the subbasin and exploring the changes within the subbasin’s elements, by calculating the percentage change between the two simulated runoff volume values (1989 and 2009) of every plane element. Only the most pronounced changes in runoff volumes were of interest, so only plane elements that had a value of above or below 1 standard deviation (SD) percent change were considered. The third step was to determine how LCCs affected the runoff changes in these plane elements. To do so, a zonal statistic tool for majority calculation of land-cover classes was used [52] for each plane element that was extracted in step 2. This means that each plane element (lower or higher than 1 SD of percent change) received one value of the land-cover class that had the largest representation within it. This was done for both time period images so that a qualitative comparison could be performed for each plane element that exhibited runoff changes satisfying the threshold of ±1 SD.The spatial analysis process was performed in order to gain a better understanding regarding the environmental factors spatially affecting the runoff trends and was conducted only for the 2009 land-cover map. For every subbasin, the simulated runoff volume values were considered as the dependent variable. Three explanatory variables were examined as predictors of the runoff volume: land-cover, represented by the land-cover classification map, elevation, represented by DEM, and vegetation cover, represented by a normalized difference vegetation index (NDVI) map [53], extracted from the 2009 image. In order to determine the relationship characteristics between these explanatory variables and the runoff volume, a geographically weighted regression (GWR) was conducted. This method is designed to consider spatial data and nonstationarity, where the variable estimates vary locally, and spatially varying relationships can be explored [54, 55]. The results were expected to point out the relations between the runoff volume and the three explanatory variables, represented by regression coefficient R
2, for each subbasin. The main goal was to understand which variable, or combination of variables, affects runoff volume variability the most for each subbasin.Seven GWR tests were built and considered for each subbasin. They were conducted for the runoff volume variable against each independent variable separately, followed by modeling paired independent variables against runoff volume and finally by running all three independent variables. The model with the highestR
2 for each subbasin implies which variable, or combination of variables, had the strongest effect on runoff.


## 3. Results and Discussion

### 3.1. Model Input Data

Figures4(a) and 4(b) show the Yarkon-Ayalon watershed after the classification process for the two time periods investigated. The overall classification accuracy and Kappa statistics for the two classification maps are 81.3% and 0.75 for the 1989 map and 82.97% and 0.78 for the 2009 map, respectively. Using the Thiessen polygons method, followed by weighting each interpolated polygon relative to its portion within each subbasin, six rainfall depth (mm) values were extracted, one for every subbasin. The totals of these values range between 155.3 mm and 172.1 mm for the five-day event (December 29, 1991–January 3, 1992).Figure 4
Land-cover classification maps of (a) 1989 and (b) 2009.


### 3.2. Model Calibration, Validation, and Output

The 1989 model results are presented in Table2, describing the simulated and predicted runoff volume values. The 1989 model output values were introduced into the OF equations, along within situ measurements, in order to validate the model’s accuracy.Table 2
Observed total runoff volume (million cubic meters (MCM)) for the storm event that was used for model calibration. The right column shows the simulated runoff values from the 1989 model after calibration.
Subbasin
1991-1992 observed values (MCM)
1989 values (MCM) after calibration
Yarkon
42.93
42.77
Beit Dagan
30.63
30.35
Shilo
10.72
10.77
Natuf
16.74
16.65
Kana
8.37
8.39
Lod
8.81
8.88The OF results are shown in Table3, each calculated for all six subbasins. The relationships between the OF results and the subbasin sizes were studied using a regression model, supported by significance levels (P value). Figures 5(a)–5(e) display these relations, along with the supporting statistics. Mean absolute discharge deviation (MAQD), root mean square discharge deviation (RMS), and maximum absolute discharge deviation (MXAQ) all show linear, significant relationships to subbasin size, with R
2 values of 0.9, 0.9, and 0.95, respectively. The peak discharge ratio (PQR) and absolute peak deviation (ADPQ) OFs are both measures of peak discharges and exhibit a nonsignificant relationship to subbasin size. It can be concluded, from the OF relationships to the subbasin areas, that as the latter increases, so does the deviation between the model simulations andin situ measurements concerning discharge throughout the event (Figures 5(a)–5(c)). When dealing with peak discharge deviation (Figures 5(d) and 5(e)), however, there is no significance to the subbasin size, and the modeled results show large inaccuracies.Table 3
The results of the objective functions (OFs).
Subbasin
Size (km2)
Total runoff ratio (VOL)
Peak discharge ratio (PQR)
Absolute peak deviation (ADPQ)
Mean absolute discharge deviation (MAQD)
Root mean square discharge deviation (RMS)
Maximum absolute discharge deviation (MXAQ)
Yarkon
947
0.36
−14.08
80.48
138.71
182.73
356
Beit Dagan
557
0.93
−2.96
11.48
80.16
104.77
224.37
Shilo
357
−0.45
35.85
50.66
30.19
43.67
126.4
Natuf
242
0.53
−13.64
25.01
44.56
57.03
125.46
Kana
236
−0.22
−1.632
1.7
22
35.12
94.25
Lod
134
−0.77
10.81
11.77
20.24
29.28
83.62Figure 5
Evaluation of (a) peak discharge ratio (PQR); (b) absolute deviation between peak discharges (ADPQ); (c) mean absolute discharge deviation (MAQD); (d) root mean square discharge deviation (RMS); and (e) maximum absolute discharge deviation (MXAQ) as functions of the subbasin area by means of objective functions.
(a)
(b)
(c)
(d)
(e)The validation results demonstrate that a smaller subbasin is preferable to use with this model [56], as its simulations are most likely to perform more accurately.


### 3.3. Temporal LCC Analysis

The first factor of the temporal analysis includes examination of LCCs. Comparison results of the land-cover trends between 1989 and 2009 for each subbasin are shown in Table4. For all of the subbasins, the residential class went through the largest change. The number of pixels classified as residential was more than doubled during the 19 years examined. Lod shows the highest increase in the number of residential pixels by 2009, more than tripling the number of pixels in 1989. However, when considering the percent change relative to the entire subbasin area, a different trend emerges. The agriculture land-cover class in three of the subbasins went through the largest change, showing a decreasing trend between −6% and −13%. The other large trend belongs to the residential land-cover class, implying a shift from agricultural lands to residential areas. Shilo shows a quite different behavior, with a relatively large decrease in natural area (−15%) and an increase of bare rock (13.1%), due to urbanization processes that involve vegetation removal. Next, the transition phases for each pixel were derived by overlaying the land-cover classification maps. Figure 6 shows a biplot of PCA results. Components 1 and 2 explain a total of 85.2% of the data variance. The direction of each transition phase vector points at the subbasin it characterizes the most, while the length of the vectors expresses the portion of the transition occurrence relative to the other vectors. All of the subbasins are mostly characterized by no change at all, while Lod and Yarkon experienced the slightest change. The Yarkon subbasin experienced transitions from agricultural land, natural area, and orchards to residential area. The Shilo subbasin had greater transitions from natural areas, agricultural land, and orchards to bare rock than any other subbasin, along with transitions from natural area to orchards. The Natuf subbasin shows the highest transition portions from bare rock to orchards and natural areas. The main transition in the Kana subbasin was from agricultural land to orchards. The Lod and Beit Dagan subbasins did not experience any major transition, apart from agricultural land transformation to residential areas, as did the Yarkon subbasin.Table 4
Changes in land-cover classes between 1989 and 2009. Values in bold mark the largest changes.
Land-cover
Number of pixels in 1989
Number of pixels in 2009
Change of specific land-cover between years (%)
Change out of total subbasin (%)
Number of pixels in 1989
Number of pixels in 2009
Change of specific land-cover between years (%)
Change out of total subbasin (%)
Yarkon
Beit Dagan
Residential
45830
116316
153.80
6.97
30496
79637
161.14
8.27
Bare rock
92196
145011
57.29
5.23
88798
63719
−28.24
−4.22
Forest
5022
3106
−38.15
−0.19
7302
8132
11.37
0.14
Orchards
283741
315179
11.08
3.11
123224
121657
−1.27
−0.26
Natural area
290190
259873
−10.45
−3.00
159854
182431
14.12
3.80
Agricultural land
293927
171308
−41.72
−12.13
184744
138838
−24.85
−7.72
Shilo
Natuf
Residential
9208
19917
116.30
2.81
6955
18969
172.74
4.65
Bare rock
41795
91806
119.66
13.12
55905
49271
−11.87
−2.57
Forest
765
434
−43.27
−0.09
1343
3015
124.50
0.65
Orchards
115226
146357
27.02
8.17
72206
82477
14.22
3.97
Natural area
132613
75376
−43.16
−15.02
81506
80775
−0.90
−0.28
Agricultural land
81483
47199
−42.08
−9.00
40662
24070
−40.80
−6.42
Kana
Lod
Residential
6620
18149
174.15
4.57
7664
25545
233.31
12.50
Bare rock
22029
35127
59.46
5.19
10374
6229
−39.96
−2.90
Forest
839
760
−9.42
−0.03
4581
4048
−11.64
−0.37
Orchards
92137
103379
12.20
4.46
25067
12789
−48.98
−8.58
Natural area
80675
77806
−3.56
−1.14
33743
42701
26.55
6.26
Agricultural land
49997
16994
−66.01
−13.09
61650
51764
−16.04
−6.91Figure 6
Principal component analysis (PCA) biplot of the transition phases’ data. The land-cover classes represented are residential areas (res), agricultural land (crop), orchards (orch), natural areas (nat), and bare exposed rocks (rock).It is worth mentioning that the decrease in orchard land-cover occurred mostly in the Israeli side of the watershed, while increase in orchards was observed mainly along the PA part of the watershed. This can be explained by the extensive urbanization process in the western part of the watershed that led to a reduction in orchards in favor of urban development. Along the West Bank, however, orchards are an effective mean for Palestinian farmers to establish facts on the ground, thus the increase in this land-cover type along the Palestinian side of the watershed.


### 3.4. Temporal Runoff Change

The second factor of temporal analysis, the comparison between simulated runoff volume and peak discharge values for the two examined time periods, is presented in Table5, along with the percent change between the simulated values. All of the subbasins show a moderate increase in both runoff volumes and peak discharge, aside from the Natuf subbasin that presents an opposite trend. This trend corresponds to the findings of the Israeli Hydrological Service reports during the past 40 years [57].Table 5
Simulated runoff volume (million cubic meters (MCM)) and peak discharge values (m3/s) extracted from the 1989 and 2009 models.
Subbasin
1989 simulated value
2009 simulated value
Percent change between simulated values
Yarkon
Runoff volume (MCM)
42.77
43.74
2.26
Peak discharge (m3/s)
571.48
582
1.84
Beit Dagan
Runoff volume (MCM)
30.35
31.96
5.3
Peak discharge (m3/s)
388.48
405.99
4.51
Shilo
Runoff volume (MCM)
10.77
11.39
5.79
Peak discharge (m3/s)
141.34
147.245
4.18
Natuf
Runoff volume (MCM)
16.65
16.56
−0.57
Peak discharge (m3/s)
183.66
182.71
−0.52
Kana
Runoff volume (MCM)
8.39
8.78
4.64
Peak discharge (m3/s)
104.1
108.55
4.28
Lod
Runoff volume (MCM)
8.88
9.47
6.69
Peak discharge (m3/s)
108.93
114.68
5.28The percent changes in runoff volume above or below 1 SD within the planes are shown in Table6. It exhibits the portion of the area that had experienced a change of ±1 SD from the whole subbasin size. Figure 7 presents this data graphically. The Shilo subbasin went through the largest total change of ±1 SD, with 30% and 5% increase and decrease changes, respectively, followed by the Natuf subbasin, in which 20% of its total area experienced an increase in runoff values larger than 1 SD, as opposed to only a 10% decrease of those below 1 SD. At the same time, the Natuf subbasin had the largest decrease change value, compared to the other subbasins. The other cases show similar trends, with about a 9–18% increase change and a 5–8% decrease change of runoff volume above or below 1 SD.Table 6
Area size and percent change of ±1 SD increase, decrease, and total change in runoff volume values.
Subbasin
1 SD of % runoff change area
Subbasin area (km2)
Increase change area (km2)
Increase change portion
Decrease change area (km2)
Decrease change portion
Total change area (km2)
Total change
Yarkon
3.45
946.7
123.0
0.13
2.1
0.00
125.1
0.13
Beit Dagan
66.87
556.7
80.9
0.15
0.0
0.00
80.9
0.15
Shilo
7.5
356.9
110.4
0.31
18.0
0.05
128.4
0.36
Natuf
2.95
242.1
48.4
0.20
24.2
0.10
72.6
0.30
Kana
16.4
236.3
23.4
0.10
18.4
0.08
41.8
0.18
Lod
16.1
134
24.4
0.18
0.1
0.00
24.5
0.18Figure 7
Increase, decrease, and total runoff volume change of ±1 standard deviation (SD) out of the total subbasin area, within the planes.Examples of the temporal analysis results, for the two of the six subbasins examined, are presented in Figures8 and 9. Maps in Figures 8(a) and 9(a) present the positive and negative changes in runoff volume of ±1 SD. It should be noted that while Natuf exhibits a trend toward increasing runoff where these values are limited to ±1 SD, the general trend of runoff change (total increase and decrease percent change) is 47.3% increase and 52.7% decrease. This is complementary with Table 5 that implies a decreasing trend in runoff volume for this subbasin.Figure 8
Temporal land-cover changes (LCCs) in the Lod subbasin between the years 1989 and 2009: (a) change in simulated runoff volume of ±1 SD; (b) land-cover status in 1989, based on majority zonal statistics calculations; (c) land-cover status in 2009, based on majority zonal statistics calculations.
(a)
(b)
(c)Figure 9
Temporal land-cover changes (LCCs) in the Natuf subbasin between the years 1989 and 2009: (a) change in simulated runoff volume of ±1 SD; (b) land-cover status in 1989, based on majority zonal statistics calculations; (c) land-cover status in 2009, based on majority zonal statistics calculations.
(a)
(b)
(c)A land-cover comparison is shown in Figures8(b), 8(c), 9(b) and 9(c) through majority zonal statistics calculations (each polygon receives the value of the most common land-cover class within it) for the land-cover status of the polygons undergoing ±1 SD of runoff volume change (%), where maps in Figures 8(b) and 9(b) are the 1989 status and maps in Figures 8(c) and 9(c) represent 2009. This analysis was conducted for all six subbasins (Table 6), providing the following outcomes.(i)
Yarkon: most of the area that changed ±1 SD of the total percent runoff change has experienced a positive change in runoff volumes; this trend is owing to an increase in residential areas, mostly downstream towards the coastline. The shifts are usually from agricultural areas to built-up landscapes.
(ii)
Beit Dagan: there was no decrease in percent runoff change below 1 SD; rather, all of the changes are larger than 1 SD due to alterations from agricultural activities to residential land-cover.
(iii)
Shilo: most of the large changes are positive, located upstream along the headwaters. The land-cover shifts are of different characteristics. Some natural areas changed into orchards or bare exposed rock and so did the agricultural land. This greatly implies that the growing built-up trend along the headwaters is due to the activities of Palestinian villages and cities. The rocky nature of the shift is due to quarries, road paving, preparation of the area for urban development, and a general exposure of the ground.
(iv)
Lod (Figure8): this subbasin experienced mostly increase of larger than 1 SD runoff change, with an insignificant portion of the area going through decrease in lower than 1 SD runoff change (0.1
k
m
2). This is a result of a transformation of an agricultural landscape to an urbanized landscape.
(v)
Natuf (Figure9): similar to the Shilo subbasin, this area also experienced an increase in runoff change along the headwaters, shifting from natural and orchard land-cover to bare exposed rocks. However, the lower areas show a decrease in runoff change below 1 SD, due to a shift from bare exposed rock to natural covered area. The lowest, western part of this subbasin expresses a strong increase in runoff, mostly due to urbanization development at the expense of agricultural activities.
(vi)
Kana: in this subbasin, similar portions of the area went through both an increase and a decrease in the amounts of runoff (10% and 8%, resp.). The increase was mostly due to urbanization where some natural cover or agricultural activities used to take place, while the decreasing trend is a result of orchard-covered areas’ transformation to natural covered landscape.Most increase changes are explicitly due to a transformation from agricultural activities or natural land-cover to residential uses or bare rocks, mostly occurring along the rural-urban fringe and known as “urban-sprawl,” on the Israeli side of the watershed. This might be due to several reasons, including immigration and increase in population, natural growth, economic processes, and planning policies [58]. The Palestinian part of the watershed has also experienced a constant and rapid urban growth. The Ramallah governorate, which is mostly contained within the Yarkon-Ayalon watershed and affects the Natuf and Shilo subbasins, has a population size that constitutes about 12% of the total population in the West Bank, with an urban expansion reaching its limits in all directions of the municipal boundary. Between 1989 and 1994 the built-up area of Ramallah expanded by 16.1% and, between 1994 and 2000, it added another 24.5% to its area [29]. The inevitable outcome, in most cases of increased urban areas, is increases in runoff volumes, peak discharges, and flash floods [14, 59].


### 3.5. Geographically Weighted Regression (GWR) Statistical Tests

The spatial analysis results for the GWR statistical tests are shown in Table7, where the final adjusted R
2 values for each of the statistical models that were performed are presented. The GWR results indicate that, in all of the subbasins aside from Shilo, the vegetation cover variable (represented by NDVI values) alone was the best explanatory variable for the runoff variation along these subbasins (in Lod subbasin the coefficient of determination (R
2) value for NDVI was slightly higher than the value received from all of the explanatory variables combined). The NDVI variable was able to explain 89 to 93% of the runoff volume spatial variance. This corresponds to a large number of studies in this field, concluding that vegetation cover decrease has a great effect on runoff volumes, runoff peak discharge [10, 60], and the lag time between the rainfall event and runoff generation as well as on soil erosion and nutrient loss [61, 62]. For the Shilo subbasin, all of the explanatory variables put together produced the highest regression coefficient; while vegetation cover is the main contributor, the elevation variable provides an additional explanation to the runoff volume variance, owing to the steep topography along the headwaters of this subbasin.Table 7
Geographically weighted regression (GWR) results for runoff volume as the dependent variable and the digital elevation model (DEM), land-cover classes, and the normalized difference vegetation index (NDVI) as explanatory variables. The bold highlighted values are the highestR
2 results for each subbasin.
Spatial statistics and explanatory variables
Yarkon
Beit Dagan
Shilo
Natuf
Kana
Lod
GWR-classification
0.9
0.86
0.82
0.78
0.91
0.82
GWR-DEM
0.65
0.66
0.77
0.43
0.67
0.59
GWR-NDVI
0.93
0.92
0.94
0.89
0.92
0.93
GWR-DEM + classification
0.57
0.63
0.96
0.39
0.72
0.92
GWR-NDVI + classification
0.88
0.84
0.73
0.77
0.89
0.79
GWR-NDVI + DEM
0.59
0.63
0.96
0.37
0.72
0.92
GWR-all
0.5
0.60
0.97
0.87
0.66
0.93


## 3.1. Model Input Data

Figures4(a) and 4(b) show the Yarkon-Ayalon watershed after the classification process for the two time periods investigated. The overall classification accuracy and Kappa statistics for the two classification maps are 81.3% and 0.75 for the 1989 map and 82.97% and 0.78 for the 2009 map, respectively. Using the Thiessen polygons method, followed by weighting each interpolated polygon relative to its portion within each subbasin, six rainfall depth (mm) values were extracted, one for every subbasin. The totals of these values range between 155.3 mm and 172.1 mm for the five-day event (December 29, 1991–January 3, 1992).Figure 4
Land-cover classification maps of (a) 1989 and (b) 2009.


## 3.2. Model Calibration, Validation, and Output

The 1989 model results are presented in Table2, describing the simulated and predicted runoff volume values. The 1989 model output values were introduced into the OF equations, along within situ measurements, in order to validate the model’s accuracy.Table 2
Observed total runoff volume (million cubic meters (MCM)) for the storm event that was used for model calibration. The right column shows the simulated runoff values from the 1989 model after calibration.
Subbasin
1991-1992 observed values (MCM)
1989 values (MCM) after calibration
Yarkon
42.93
42.77
Beit Dagan
30.63
30.35
Shilo
10.72
10.77
Natuf
16.74
16.65
Kana
8.37
8.39
Lod
8.81
8.88The OF results are shown in Table3, each calculated for all six subbasins. The relationships between the OF results and the subbasin sizes were studied using a regression model, supported by significance levels (P value). Figures 5(a)–5(e) display these relations, along with the supporting statistics. Mean absolute discharge deviation (MAQD), root mean square discharge deviation (RMS), and maximum absolute discharge deviation (MXAQ) all show linear, significant relationships to subbasin size, with R
2 values of 0.9, 0.9, and 0.95, respectively. The peak discharge ratio (PQR) and absolute peak deviation (ADPQ) OFs are both measures of peak discharges and exhibit a nonsignificant relationship to subbasin size. It can be concluded, from the OF relationships to the subbasin areas, that as the latter increases, so does the deviation between the model simulations andin situ measurements concerning discharge throughout the event (Figures 5(a)–5(c)). When dealing with peak discharge deviation (Figures 5(d) and 5(e)), however, there is no significance to the subbasin size, and the modeled results show large inaccuracies.Table 3
The results of the objective functions (OFs).
Subbasin
Size (km2)
Total runoff ratio (VOL)
Peak discharge ratio (PQR)
Absolute peak deviation (ADPQ)
Mean absolute discharge deviation (MAQD)
Root mean square discharge deviation (RMS)
Maximum absolute discharge deviation (MXAQ)
Yarkon
947
0.36
−14.08
80.48
138.71
182.73
356
Beit Dagan
557
0.93
−2.96
11.48
80.16
104.77
224.37
Shilo
357
−0.45
35.85
50.66
30.19
43.67
126.4
Natuf
242
0.53
−13.64
25.01
44.56
57.03
125.46
Kana
236
−0.22
−1.632
1.7
22
35.12
94.25
Lod
134
−0.77
10.81
11.77
20.24
29.28
83.62Figure 5
Evaluation of (a) peak discharge ratio (PQR); (b) absolute deviation between peak discharges (ADPQ); (c) mean absolute discharge deviation (MAQD); (d) root mean square discharge deviation (RMS); and (e) maximum absolute discharge deviation (MXAQ) as functions of the subbasin area by means of objective functions.
(a)
(b)
(c)
(d)
(e)The validation results demonstrate that a smaller subbasin is preferable to use with this model [56], as its simulations are most likely to perform more accurately.


## 3.3. Temporal LCC Analysis

The first factor of the temporal analysis includes examination of LCCs. Comparison results of the land-cover trends between 1989 and 2009 for each subbasin are shown in Table4. For all of the subbasins, the residential class went through the largest change. The number of pixels classified as residential was more than doubled during the 19 years examined. Lod shows the highest increase in the number of residential pixels by 2009, more than tripling the number of pixels in 1989. However, when considering the percent change relative to the entire subbasin area, a different trend emerges. The agriculture land-cover class in three of the subbasins went through the largest change, showing a decreasing trend between −6% and −13%. The other large trend belongs to the residential land-cover class, implying a shift from agricultural lands to residential areas. Shilo shows a quite different behavior, with a relatively large decrease in natural area (−15%) and an increase of bare rock (13.1%), due to urbanization processes that involve vegetation removal. Next, the transition phases for each pixel were derived by overlaying the land-cover classification maps. Figure 6 shows a biplot of PCA results. Components 1 and 2 explain a total of 85.2% of the data variance. The direction of each transition phase vector points at the subbasin it characterizes the most, while the length of the vectors expresses the portion of the transition occurrence relative to the other vectors. All of the subbasins are mostly characterized by no change at all, while Lod and Yarkon experienced the slightest change. The Yarkon subbasin experienced transitions from agricultural land, natural area, and orchards to residential area. The Shilo subbasin had greater transitions from natural areas, agricultural land, and orchards to bare rock than any other subbasin, along with transitions from natural area to orchards. The Natuf subbasin shows the highest transition portions from bare rock to orchards and natural areas. The main transition in the Kana subbasin was from agricultural land to orchards. The Lod and Beit Dagan subbasins did not experience any major transition, apart from agricultural land transformation to residential areas, as did the Yarkon subbasin.Table 4
Changes in land-cover classes between 1989 and 2009. Values in bold mark the largest changes.
Land-cover
Number of pixels in 1989
Number of pixels in 2009
Change of specific land-cover between years (%)
Change out of total subbasin (%)
Number of pixels in 1989
Number of pixels in 2009
Change of specific land-cover between years (%)
Change out of total subbasin (%)
Yarkon
Beit Dagan
Residential
45830
116316
153.80
6.97
30496
79637
161.14
8.27
Bare rock
92196
145011
57.29
5.23
88798
63719
−28.24
−4.22
Forest
5022
3106
−38.15
−0.19
7302
8132
11.37
0.14
Orchards
283741
315179
11.08
3.11
123224
121657
−1.27
−0.26
Natural area
290190
259873
−10.45
−3.00
159854
182431
14.12
3.80
Agricultural land
293927
171308
−41.72
−12.13
184744
138838
−24.85
−7.72
Shilo
Natuf
Residential
9208
19917
116.30
2.81
6955
18969
172.74
4.65
Bare rock
41795
91806
119.66
13.12
55905
49271
−11.87
−2.57
Forest
765
434
−43.27
−0.09
1343
3015
124.50
0.65
Orchards
115226
146357
27.02
8.17
72206
82477
14.22
3.97
Natural area
132613
75376
−43.16
−15.02
81506
80775
−0.90
−0.28
Agricultural land
81483
47199
−42.08
−9.00
40662
24070
−40.80
−6.42
Kana
Lod
Residential
6620
18149
174.15
4.57
7664
25545
233.31
12.50
Bare rock
22029
35127
59.46
5.19
10374
6229
−39.96
−2.90
Forest
839
760
−9.42
−0.03
4581
4048
−11.64
−0.37
Orchards
92137
103379
12.20
4.46
25067
12789
−48.98
−8.58
Natural area
80675
77806
−3.56
−1.14
33743
42701
26.55
6.26
Agricultural land
49997
16994
−66.01
−13.09
61650
51764
−16.04
−6.91Figure 6
Principal component analysis (PCA) biplot of the transition phases’ data. The land-cover classes represented are residential areas (res), agricultural land (crop), orchards (orch), natural areas (nat), and bare exposed rocks (rock).It is worth mentioning that the decrease in orchard land-cover occurred mostly in the Israeli side of the watershed, while increase in orchards was observed mainly along the PA part of the watershed. This can be explained by the extensive urbanization process in the western part of the watershed that led to a reduction in orchards in favor of urban development. Along the West Bank, however, orchards are an effective mean for Palestinian farmers to establish facts on the ground, thus the increase in this land-cover type along the Palestinian side of the watershed.


## 3.4. Temporal Runoff Change

The second factor of temporal analysis, the comparison between simulated runoff volume and peak discharge values for the two examined time periods, is presented in Table5, along with the percent change between the simulated values. All of the subbasins show a moderate increase in both runoff volumes and peak discharge, aside from the Natuf subbasin that presents an opposite trend. This trend corresponds to the findings of the Israeli Hydrological Service reports during the past 40 years [57].Table 5
Simulated runoff volume (million cubic meters (MCM)) and peak discharge values (m3/s) extracted from the 1989 and 2009 models.
Subbasin
1989 simulated value
2009 simulated value
Percent change between simulated values
Yarkon
Runoff volume (MCM)
42.77
43.74
2.26
Peak discharge (m3/s)
571.48
582
1.84
Beit Dagan
Runoff volume (MCM)
30.35
31.96
5.3
Peak discharge (m3/s)
388.48
405.99
4.51
Shilo
Runoff volume (MCM)
10.77
11.39
5.79
Peak discharge (m3/s)
141.34
147.245
4.18
Natuf
Runoff volume (MCM)
16.65
16.56
−0.57
Peak discharge (m3/s)
183.66
182.71
−0.52
Kana
Runoff volume (MCM)
8.39
8.78
4.64
Peak discharge (m3/s)
104.1
108.55
4.28
Lod
Runoff volume (MCM)
8.88
9.47
6.69
Peak discharge (m3/s)
108.93
114.68
5.28The percent changes in runoff volume above or below 1 SD within the planes are shown in Table6. It exhibits the portion of the area that had experienced a change of ±1 SD from the whole subbasin size. Figure 7 presents this data graphically. The Shilo subbasin went through the largest total change of ±1 SD, with 30% and 5% increase and decrease changes, respectively, followed by the Natuf subbasin, in which 20% of its total area experienced an increase in runoff values larger than 1 SD, as opposed to only a 10% decrease of those below 1 SD. At the same time, the Natuf subbasin had the largest decrease change value, compared to the other subbasins. The other cases show similar trends, with about a 9–18% increase change and a 5–8% decrease change of runoff volume above or below 1 SD.Table 6
Area size and percent change of ±1 SD increase, decrease, and total change in runoff volume values.
Subbasin
1 SD of % runoff change area
Subbasin area (km2)
Increase change area (km2)
Increase change portion
Decrease change area (km2)
Decrease change portion
Total change area (km2)
Total change
Yarkon
3.45
946.7
123.0
0.13
2.1
0.00
125.1
0.13
Beit Dagan
66.87
556.7
80.9
0.15
0.0
0.00
80.9
0.15
Shilo
7.5
356.9
110.4
0.31
18.0
0.05
128.4
0.36
Natuf
2.95
242.1
48.4
0.20
24.2
0.10
72.6
0.30
Kana
16.4
236.3
23.4
0.10
18.4
0.08
41.8
0.18
Lod
16.1
134
24.4
0.18
0.1
0.00
24.5
0.18Figure 7
Increase, decrease, and total runoff volume change of ±1 standard deviation (SD) out of the total subbasin area, within the planes.Examples of the temporal analysis results, for the two of the six subbasins examined, are presented in Figures8 and 9. Maps in Figures 8(a) and 9(a) present the positive and negative changes in runoff volume of ±1 SD. It should be noted that while Natuf exhibits a trend toward increasing runoff where these values are limited to ±1 SD, the general trend of runoff change (total increase and decrease percent change) is 47.3% increase and 52.7% decrease. This is complementary with Table 5 that implies a decreasing trend in runoff volume for this subbasin.Figure 8
Temporal land-cover changes (LCCs) in the Lod subbasin between the years 1989 and 2009: (a) change in simulated runoff volume of ±1 SD; (b) land-cover status in 1989, based on majority zonal statistics calculations; (c) land-cover status in 2009, based on majority zonal statistics calculations.
(a)
(b)
(c)Figure 9
Temporal land-cover changes (LCCs) in the Natuf subbasin between the years 1989 and 2009: (a) change in simulated runoff volume of ±1 SD; (b) land-cover status in 1989, based on majority zonal statistics calculations; (c) land-cover status in 2009, based on majority zonal statistics calculations.
(a)
(b)
(c)A land-cover comparison is shown in Figures8(b), 8(c), 9(b) and 9(c) through majority zonal statistics calculations (each polygon receives the value of the most common land-cover class within it) for the land-cover status of the polygons undergoing ±1 SD of runoff volume change (%), where maps in Figures 8(b) and 9(b) are the 1989 status and maps in Figures 8(c) and 9(c) represent 2009. This analysis was conducted for all six subbasins (Table 6), providing the following outcomes.(i)
Yarkon: most of the area that changed ±1 SD of the total percent runoff change has experienced a positive change in runoff volumes; this trend is owing to an increase in residential areas, mostly downstream towards the coastline. The shifts are usually from agricultural areas to built-up landscapes.
(ii)
Beit Dagan: there was no decrease in percent runoff change below 1 SD; rather, all of the changes are larger than 1 SD due to alterations from agricultural activities to residential land-cover.
(iii)
Shilo: most of the large changes are positive, located upstream along the headwaters. The land-cover shifts are of different characteristics. Some natural areas changed into orchards or bare exposed rock and so did the agricultural land. This greatly implies that the growing built-up trend along the headwaters is due to the activities of Palestinian villages and cities. The rocky nature of the shift is due to quarries, road paving, preparation of the area for urban development, and a general exposure of the ground.
(iv)
Lod (Figure8): this subbasin experienced mostly increase of larger than 1 SD runoff change, with an insignificant portion of the area going through decrease in lower than 1 SD runoff change (0.1
k
m
2). This is a result of a transformation of an agricultural landscape to an urbanized landscape.
(v)
Natuf (Figure9): similar to the Shilo subbasin, this area also experienced an increase in runoff change along the headwaters, shifting from natural and orchard land-cover to bare exposed rocks. However, the lower areas show a decrease in runoff change below 1 SD, due to a shift from bare exposed rock to natural covered area. The lowest, western part of this subbasin expresses a strong increase in runoff, mostly due to urbanization development at the expense of agricultural activities.
(vi)
Kana: in this subbasin, similar portions of the area went through both an increase and a decrease in the amounts of runoff (10% and 8%, resp.). The increase was mostly due to urbanization where some natural cover or agricultural activities used to take place, while the decreasing trend is a result of orchard-covered areas’ transformation to natural covered landscape.Most increase changes are explicitly due to a transformation from agricultural activities or natural land-cover to residential uses or bare rocks, mostly occurring along the rural-urban fringe and known as “urban-sprawl,” on the Israeli side of the watershed. This might be due to several reasons, including immigration and increase in population, natural growth, economic processes, and planning policies [58]. The Palestinian part of the watershed has also experienced a constant and rapid urban growth. The Ramallah governorate, which is mostly contained within the Yarkon-Ayalon watershed and affects the Natuf and Shilo subbasins, has a population size that constitutes about 12% of the total population in the West Bank, with an urban expansion reaching its limits in all directions of the municipal boundary. Between 1989 and 1994 the built-up area of Ramallah expanded by 16.1% and, between 1994 and 2000, it added another 24.5% to its area [29]. The inevitable outcome, in most cases of increased urban areas, is increases in runoff volumes, peak discharges, and flash floods [14, 59].


## 3.5. Geographically Weighted Regression (GWR) Statistical Tests

The spatial analysis results for the GWR statistical tests are shown in Table7, where the final adjusted R
2 values for each of the statistical models that were performed are presented. The GWR results indicate that, in all of the subbasins aside from Shilo, the vegetation cover variable (represented by NDVI values) alone was the best explanatory variable for the runoff variation along these subbasins (in Lod subbasin the coefficient of determination (R
2) value for NDVI was slightly higher than the value received from all of the explanatory variables combined). The NDVI variable was able to explain 89 to 93% of the runoff volume spatial variance. This corresponds to a large number of studies in this field, concluding that vegetation cover decrease has a great effect on runoff volumes, runoff peak discharge [10, 60], and the lag time between the rainfall event and runoff generation as well as on soil erosion and nutrient loss [61, 62]. For the Shilo subbasin, all of the explanatory variables put together produced the highest regression coefficient; while vegetation cover is the main contributor, the elevation variable provides an additional explanation to the runoff volume variance, owing to the steep topography along the headwaters of this subbasin.Table 7
Geographically weighted regression (GWR) results for runoff volume as the dependent variable and the digital elevation model (DEM), land-cover classes, and the normalized difference vegetation index (NDVI) as explanatory variables. The bold highlighted values are the highestR
2 results for each subbasin.
Spatial statistics and explanatory variables
Yarkon
Beit Dagan
Shilo
Natuf
Kana
Lod
GWR-classification
0.9
0.86
0.82
0.78
0.91
0.82
GWR-DEM
0.65
0.66
0.77
0.43
0.67
0.59
GWR-NDVI
0.93
0.92
0.94
0.89
0.92
0.93
GWR-DEM + classification
0.57
0.63
0.96
0.39
0.72
0.92
GWR-NDVI + classification
0.88
0.84
0.73
0.77
0.89
0.79
GWR-NDVI + DEM
0.59
0.63
0.96
0.37
0.72
0.92
GWR-all
0.5
0.60
0.97
0.87
0.66
0.93


## 4. Conclusions

Between the years 1989 and 2009, some changes in land-cover occurred along the Yarkon-Ayalon watershed. Although most of the considered area experienced no change at all, the major trends that were observed at locations where change took place included a decrease in agricultural lands and an increase in residential cover or residential-related processes. Changes in land-cover, through this period of 19 years, were spotted via remote-sensing tools and coupled with hydrological simulation tools in order to find the relationship between LCC and rainfall-runoff relationships. LCCs were found to have some effect on the runoff volumes and peak discharges. In locations in which land-cover shifts towards natural covered areas occurred, a minor decrease in runoff was observed, thus pointing once more to the strong relationship between vegetative cover and runoff volumes and peak discharge values. The LCCs that had the largest impact on the runoff regime were those related to urbanization and vegetation removal.Spatially, this research demonstrates high correlations between vegetation cover and runoff volume values within the study area, implying that a continuing trend of vegetation removal will result in higher runoff volumes across the watershed during extreme rainfall events. Soil erosion trends were not examined due to a lack of validation data, although it can be expected that they would correspond to the runoff findings.An extreme rainfall event was used to examine the outcomes of LCCs on runoff volumes and peak discharges. While this might highlight different phenomena across the watershed, it would be interesting and important to examine the effects of frequent rainfall events on the runoff trends throughout periods of changing land-covers, since these events are far more frequent and have a greater impact on the watershed area.



---

*Source: 838070-2015-10-20.xml*