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.

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 [13]. 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 [58]. 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, 912].

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 [1417]. 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 (Figure 1(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].

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 (Figure 2). 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).

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, 3335].

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 (Figure 2). Each of these subbasins was modeled separately, and model results were analyzed. Figure 3 shows the model’s flowchart that is specified hereinafter.

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. The in 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.

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 , 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 highest 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

Figures 4(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).

3.2. Model Calibration, Validation, and Output

The 1989 model results are presented in Table 2, describing the simulated and predicted runoff volume values. The 1989 model output values were introduced into the OF equations, along with in situ measurements, in order to validate the model’s accuracy.

The OF results are shown in Table 3, 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 ( 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 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 and in 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.

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 Table 4. 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.

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 Table 5, 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].

The percent changes in runoff volume above or below 1 SD within the planes are shown in Table 6. 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.

Examples of the temporal analysis results, for the two of the six subbasins examined, are presented in Figures 8 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.

A land-cover comparison is shown in Figures 8(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 (Figure 8): 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 (). This is a result of a transformation of an agricultural landscape to an urbanized landscape.(v)Natuf (Figure 9): 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 Table 7, where the final adjusted 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 () 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.

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.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.


This study was conducted with partial financial support from the Ministry of Agriculture and the Jewish National Fund (JNF). The authors would like to thank The Israeli Water Authority for providing data, Shea Burns and David Goodrich (of the USDA/ARS) for their invaluable assistance, and the Jacob Blaustein Institutes for Desert Research Remote Sensing Laboratory staff and students for their most appreciated support and knowledge sharing.