Abstract

Classic rainfall–runoff models usually use historical data to estimate model parameters and mean values of parameters are considered for predictions. However, due to climate changes and human effects, model parameters change temporally. To overcome this problem, normalized difference vegetation index (NDVI) derived from remotely sensed data was used in this study to investigate the effect of land cover variations on hydrological response of watersheds using a conceptual rainfall–runoff model. The study area consists of two sub-watersheds (Hervi and Lighvan) with varied land cover conditions. Obtained results show that the one-parameter model generates runoff forecasts with acceptable level of the considered criteria. Remote sensing data were employed to relate land cover properties of the watershed to the model parameter. While a power form of the regression equation could be best fitted to the parameter values using available images of Hervi sub-watershed, for the Lighvan sub-watershed the fitted equation shows somewhat lower correlation due to higher fluctuations of the model parameter. The average values of the Nash–Sutcliffe efficiency criterion of the model were obtained as 0.87 and 0.55, respectively, for Hervi and Lighvan sub-watersheds. Applying this methodology, the model's parameters might be determined using temporal NDVI values.

INTRODUCTION

The impact of land use and land cover change on runoff generated from a river basin is of major interest as the tendency towards increased change continues. Studies have shown that land use change can have a substantial effect on the hydrological response of a watershed to rainfall, resulting in quicker response times (Huang et al. 2008), greater river flow volumes (Hawley & Bledsoe 2011), and higher recurrence of floods (Braud et al. 2013).

Conceptual unit hydrograph (UH)-based models, derived from linear systems theory, have been used to investigate the effects of land use change on rainfall–runoff process (e.g., see Kang et al. 1998; Cheng & Wang 2002; Huang et al. 2008). Cheng et al. (2010) explored the effects of population density and impervious areas on the Nash model (Nash 1957) parameters for the urbanized Wu-Tu watershed. Subsequently, Huang et al. (2012) derived a power law between the Nash model parameters and population and impervious areas in the Wu-Tu watershed.

Unlike lumped models, semi-distributed geomorphologic models require large amounts of spatial data in which a combination of remote sensing, geographic information systems (GIS), and hydrological concepts can be used to represent spatial hydro-geomorphological data and information at landscape to global scales (Pickup 1995). Since geomorphologic models use the physical properties of the watershed to inform specification of values for the model's parameters, they are less dependent on the availability of adequate historical event data and can therefore be applied in ungauged watersheds with only limited rainfall–runoff records.

In support of geomorphological modeling, remotely sensed data about land cover type are widely used to derive input variables for a variety of hydrologic-response environmental models (Miller et al. 2007). Several studies assessed the effects of land use changes urbanization on rainfall–runoff response using a modeling strategy supported by remote sensing (e.g., see Verbeiren et al. 2013; Miller et al. 2014; Nourani et al. 2015).

As stated by Yang et al. (2011), during the past two decades remote sensing indices have been proposed to represent land cover types of vegetation (Soil-Adjusted Vegetation Index), bare soil (Normalized Multi-band Drought Index), impervious (Normalized Difference Built-up Index), and water area (Normalized Difference Water Index). However, well-known normalized difference vegetation index (NDVI) and leaf area index (LAI) are the two most commonly used indices to monitor and study vegetation responses to climatic changes at different scales (e.g., Wang et al. 2003; Beck et al. 2006; Eckert et al. 2015; Jiang et al. 2015). NDVI is related to the LAI and vegetation coverage (Baret & Guyot 1991; Turner et al. 1999; Goswami et al. 2015); the higher value of NDVI indicates the larger value of the LAI, and thus the higher amount of vegetation coverage. Li et al. (2009) incorporated the remotely sensed (MODIS-LAI) data into Xinanjiang rainfall–runoff model and assessed the model performance on 210 catchments in south-east Australia. The results showed that the inclusion of LAI data improved both the model calibration results as well as the daily runoff prediction in ungauged catchments. Ali et al. (2011) combined an empirical land use change model and an event scale rainfall–runoff model to quantify the impacts of potential land use change on storm–runoff generation in the Lai Nullah Basin. They calibrated the hydrologic engineering center-hydrologic modeling system (HEC-HMS) rainfall–runoff model and validated for five storm events in the study area. Their results showed a good consistency between the simulated and measured hydrographs at the outlet of the basin. Lakshmi et al. (2011) studied the influence of land surface on hydrometeorology and ecology, addressing various relational investigations among vegetation properties such as NDVI, LAI, surface temperature, and vegetation water content derived from satellite sensors. Also, they examined the effects of vegetation and its relationship with soil moisture on the simulated land–atmospheric interactions through the weather research and forecasting model. Li et al. (2012) used the Xinanjiang-ET and SIMplified HYDrology model (SIMHYD-ET) rainfall–runoff models (Zhang et al. 2009) to investigate impacts of vegetation change and climate variability on streamflow in a Southern Australian catchment. The models incorporated remotely sensed LAI data. The results of their study suggested that increase in plantations can reduce streamflow substantially, even more than climate variability. Liu et al. (2012) analyzed a correlation, based on different vegetation types, between time series of monthly NDVI and Palmer drought severity index (PDSI) during the growing season from April to October in Laohahe catchment. Their results showed that NDVI and PDSI had a good correlation, especially for shrub and grass. The highest value of correlation coefficient was in June when the vegetation is growing and the lower correlation occurred at the end of the growing season for all vegetation types. Zhou et al. (2013) investigated a modified version of Xinanjiang rainfall–runoff model incorporating dynamic remote sensing data for south-east Australian catchments. The results showed that using remote sensing LAI and albedo data in the modified Xinanjiang model led to model performance improvement in wetter bushfire impacted catchments. However, use of vegetation dynamics did not improve runoff time series simulation in a dry catchment. Törnros & Menzel (2014) derived the LAI using the NDVI data for the years 1982–2004. They also made a correlation analysis between precipitation and LAI for several land uses and each month of the year. Based on their results, LAI could be simulated as a function of precipitation. Tesemma et al. (2015a) assessed the effect of observed monthly LAI on hydrological model performance and the simulation of runoff using the variable infiltration capacity (VIC) hydrological model in the Goulburn–Broken catchment of Australia. VIC was calibrated with both observed monthly LAI and long-term mean monthly LAI, which were derived from the global land surface satellite (GLASS) LAI dataset covering the period from 1982 to 2012. Tesemma et al. (2015b) combined a nonlinear model for estimating changes in LAI due to climatic fluctuations with the VIC hydrological model to improve catchment streamflow prediction under a changing climate. The combined model was applied to 13 gauged sub-catchments with different land cover types (crop, pasture, and tree) in the Goulburn–Broken catchment, Australia. Incorporating climate-induced changes in LAI in the VIC model reduced the projected declines in streamflow and confirmed the importance of including the effects of changes in LAI in future projections of streamflow. Zhiqiang et al. (2016) classified the vegetation communities in Poyang Lake wetland using NDVI and linked it with the water regimes by a Gaussian regression model. Tian et al. (2017) examined the effects of revegetation on soil moisture using a biophysically based ecohydrological model. Vegetation covers were investigated based on four types including trees, shrub, grass, and crop. Results showed that trees consume more water than shrub and grass and therefore would result in more soil desiccation.

One of the major effects of urbanization is land cover degradation, therefore, NDVI as a remote sensing index can also be effectively used to monitor urbanization and land cover changes, especially for the arid condition of Iran, where natural watersheds are gradually changing to urbanized areas.

Whereas it is common to consider changes in storm runoff volume due to increases in impervious surfaces (Boyd et al. 1993; Shuster et al. 2005), our investigation of the literature indicates that the relations between temporal changes of NDVI and the conceptual parameters of a rainfall–runoff model are only rarely investigated. In the aforementioned research (e.g., Cheng & Wang 2002; Huang et al. 2008, 2012; Cheng et al. 2010), the effects of urbanization on parameters of the Nash model were investigated using population growth and impervious area data as indices of urbanization. This ignores other causes of increased urbanization such as industrial development. One could mention that impervious surface cover mapping is not straightforward due to heterogeneity of urban areas and possible confusion with bare soils, therefore, vegetation mapping and monitoring is easier. We examine the proposed GUHCR (geomorphological unit hydrograph based on cascade of linear reservoirs) model where the geomorphological properties of the watershed can compensate for a lack of observed event data.

Moreover, while most classical rainfall–runoff models use the historical event data to estimate the model parameters considering the average value for predictions, the effect of climate and anthropogenic changes during the time, which would lead to significant changes in the model parameters, are taken into account for detecting the temporal changes of NDVI using remotely sensed data.

The main aim of this study was to explore the link between the temporal change of NDVI and that of the parameters of the GUHCR model by looking into the consequent impact of land cover changes on flood behavior of the Hervi and Lighvan sub-watersheds in East Azerbaijan, Iran. In this regard, we tried to: (1) compute the NDVIs and determine the land cover changes from 1998 to 2012 using Landsat satellite imagery; (2) propose the rainfall–runoff model and calibrate the model parameter using available storm events within the same period; and (3) relate the change of the model parameter to NDVI values during the study period. These objectives provide the structural sub-headings used in the following Methodology, Results and Discussion sections.

STUDY AREA AND DATA

Study area

The study area comprises a 77-km2 sub-watershed of Lighvan and an approximately 60-km2 sub-watershed of Hervi located in north west Iran. Elevation ranges from 1,920 m to 3,453 m above sea level in Hervi and from 2,180 m to 3,453 m in Lighvan. The annual average rainfall in these sub-watersheds is 461 mm, mainly occurring in May to July. The average of mean daily temperature is 12.6 °C and the annual pan evaporation is equal to 1,866 mm. All mentioned values are the long-term mean annual values for the study area. Vegetation in the study area consists of grassland, farmland, and gardens. Although both sub-watersheds contain residential areas, the Hervi sub-watershed is more highly developed area (14.5% versus 10.1% for Lighvan). For the study area, Figure 1 shows an aerial photograph, the location map, and the shuttle radar topography mission (SRTM) digital elevation model (DEM) retrieved via http://data.cgiar-csi.org/srtm/tiles/.

Figure 1

The aerial photograph, the location, and the SRTM DEM map of the study area showing hydrometric (Hervi and Lighvan) and rain gauge (Lighvan) stations.

Figure 1

The aerial photograph, the location, and the SRTM DEM map of the study area showing hydrometric (Hervi and Lighvan) and rain gauge (Lighvan) stations.

Due to the degree of land use change, this area was selected to address the primary objective of investigating the relationship between NDVI change and model parameters. Table 1 shows the land cover properties of Hervi and Lighvan sub-watersheds in 1976 and 2010. As can be seen in Table 1, dense pasture is the dominant land cover in both Hervi and Lighvan sub-watersheds. The percentages of different land cover have experienced different changes. On the one hand, the percentage of dense pasture, characterized by high NDVI values, showed a considerable decrease from 1976 to 2010. The percentage of sparse pasture, on the other hand, increased considerably in both watersheds.

Table 1

The percentage of different land covers at Hervi and Lighvan sub-watersheds during the study period

Land cover Hervi
 
Lighvan
 
1976 2010 1976 2010 
Irrigate farming 18.40 8.43 19.01 9.47 
Dry farming 24.42 36.65 31.6 48.26 
Sparse pasture 3.35 27.10 2.94 21.06 
Dense pasture 44.49 3.10 37.53 4.74 
Garden 4.96 10.35 3.00 6.34 
Residential 4.38 14.46 5.92 10.13 
Total 100 100 100 100 
Land cover Hervi
 
Lighvan
 
1976 2010 1976 2010 
Irrigate farming 18.40 8.43 19.01 9.47 
Dry farming 24.42 36.65 31.6 48.26 
Sparse pasture 3.35 27.10 2.94 21.06 
Dense pasture 44.49 3.10 37.53 4.74 
Garden 4.96 10.35 3.00 6.34 
Residential 4.38 14.46 5.92 10.13 
Total 100 100 100 100 

Discharge measurements at 1-hr intervals at the Lighvan and Hervi gauges are available for 13 events from 1998 until 2012. Rainfall data were collected at Lighvan station by means of a weighing recording rain gauge.

The rainfall of this rain gauge was considered as the input rainfall for both sub-watersheds. The rainfall was assumed to be distributed uniformly over both sub-watersheds because of the small size of the sub-watersheds. The direct runoff hydrograph which can be separated by straight line, fixed based, and variable slope methods, was specified by the fixed gradient method as a simple and commonly used technique (Chow et al. 1988) for each event. It should be noted that the probable bias in this simple method would be offset in the calibration process. Thereafter, the Φ-index method, which applies the continuity equation, was selected among the other methods such as curve number (CN) or SCS to extract the excess hyetograph of each event from observed hyetographs (Chow et al. 1988). In this way, any rainfall prior to the beginning of direct runoff was taken as initial abstraction. Since the study area is small, the rainfall distribution was assumed uniform over the whole watershed. In this study, 13 storm events occurring simultaneously in both sub-watersheds (during 1998–2012) were selected to examine the proposed methodology. Assuming the linear system theory for instantaneous unit hydrograph (IUH) model, the runoff hydrographs for the two sub-watersheds were separated using the concept of linearity. Ten storm events were used for calibrating the model and three events used for verification. Characteristics of the available events are presented in Table 2.

Table 2

Hydrological characteristics of selected events for both calibration and verification datasets

  Event date Rainfall depth (mm) Runoff depth (mm)
 
Hervi Lighvan 
Calibration events 31/Jul/1998 22.72 2.10 2.32 
04/May/1999 5.66 1.72 2.40 
08/Apr/2003 26.02 1.84 1.20 
11/Jun/2003 19.55 1.82 1.47 
22/Apr/2006 5.03 0.49 0.66 
28/Jun/2006 6.78 2.51 0.45 
07/Apr/2010 7.26 0.55 0.33 
18/Jun/2010 7.36 0.81 0.39 
22/Apr/2011 13.36 0.71 0.96 
18/Jul/2012 5.31 1.68 1.72 
Verification events 27/May/2003 26.99 0.95 0.71 
22/May/2006 9.6 1.31 1.47 
28/Jun/2012 11.1 0.99 0.22 
  Event date Rainfall depth (mm) Runoff depth (mm)
 
Hervi Lighvan 
Calibration events 31/Jul/1998 22.72 2.10 2.32 
04/May/1999 5.66 1.72 2.40 
08/Apr/2003 26.02 1.84 1.20 
11/Jun/2003 19.55 1.82 1.47 
22/Apr/2006 5.03 0.49 0.66 
28/Jun/2006 6.78 2.51 0.45 
07/Apr/2010 7.26 0.55 0.33 
18/Jun/2010 7.36 0.81 0.39 
22/Apr/2011 13.36 0.71 0.96 
18/Jul/2012 5.31 1.68 1.72 
Verification events 27/May/2003 26.99 0.95 0.71 
22/May/2006 9.6 1.31 1.47 
28/Jun/2012 11.1 0.99 0.22 

Satellite imagery

A well-developed global archive of Landsat images with 30 m spatial resolution is available, and is widely used to detect and monitor land cover change (Kepner et al. 2008). Satellite images for the study area (located on path 168, row 34) for the rainy season months May (3), June (2), and July (3) in the period 1998–2012 were downloaded from the USGS website (https://earthexplorer.usgs.gov/login/). Since land cover typically does not change rapidly, one image was assigned to each rainfall–runoff event based on the closest cloud-free image to the date of event. The weather around a storm event's day is usually cloudy, also in some cases there is no image on the exact day of event. Thus we tried to select the cloud-free image of day which is closest to the event date to calculate the NDVI, representing the land cover of the study area for the certain event. Accordingly, seven images were selected for the period 1998–2012, and an extra image from July 1985 was selected to represent land cover prior to the study period. Characteristics of the selected images are presented in Table 3.

Table 3

Characteristics of remote sensing sensors and mean values of NDVI for the selected events

  Image date Landsat sensor Mean NDVI
 
Hervi Lighvan 
Calibration images 29/Jul/1985 5 TM 0.231 0.286 
01/Jul/1998 5 TM 0.219 0.283 
12/Jul/1999 7 ETM + 0.212 0.262 
04/May/2003 7 ETM + 0.156 0.123 
23/Jul/2006 5 TM 0.186 0.263 
23/May/2010 7 ETM + 0.175 0.213 
03/Jun/2011 5 TM 0.112 0.142 
29/Jun/2012 7 ETM + 0.172 0.210 
Verification images 04/May/2003 7 ETM + 0.156 0.123 
23/Jul/2006 5 TM 0.186 0.263 
17/Sep/2012 7 ETM + 0.172 0.210 
  Image date Landsat sensor Mean NDVI
 
Hervi Lighvan 
Calibration images 29/Jul/1985 5 TM 0.231 0.286 
01/Jul/1998 5 TM 0.219 0.283 
12/Jul/1999 7 ETM + 0.212 0.262 
04/May/2003 7 ETM + 0.156 0.123 
23/Jul/2006 5 TM 0.186 0.263 
23/May/2010 7 ETM + 0.175 0.213 
03/Jun/2011 5 TM 0.112 0.142 
29/Jun/2012 7 ETM + 0.172 0.210 
Verification images 04/May/2003 7 ETM + 0.156 0.123 
23/Jul/2006 5 TM 0.186 0.263 
17/Sep/2012 7 ETM + 0.172 0.210 

METHODOLOGY

Remote sensing analysis

The Landsat images were processed to compute the NDVI for each image pixel. Available cloud-free satellite images were used to detect the land cover changes of the study area.

Image processing

Since the Landsat images are from different times and seasons, the position of the sun, the angle of the terrain, and the atmospheric effects are different. Also, different sensors might not have the same image data. In order to make comparisons between the images obtained by different sensors in different times and seasons, the Landsat images have to be pre-processed by the common methods. The image pre-processing carried out in this study includes radiometric, atmospheric, and topographic corrections. The radiometric correction decreases signal variations uncorrelated to the brightness of the image surface. The purpose of atmospheric correction is to eliminate the atmospheric effects to specify true surface reflectance values. This correction is done by the FLAASH method in ENVI (ENVI 2009). For a DEM, the topographic correction removes the effects resulting from the differences in illumination due to the position of the sun and the angle of the terrain (Lu et al. 2008). This topographic correction has been already applied by the provider (see, http://data.cgiar-csi.org/srtm/tiles/) to the SRTM DEM used in this study. The Landsat images were all geo-referenced to the same co-ordinate system and there was no need to apply geometric correction as this is done by the U.S. Geological Survey using globally available digital topographic data.

Bands 3 (red visible) and 4 (near-infrared) were corrected and the NDVI was calculated for each image. The watershed boundary was then used to subset the NDVI layer via region of interests (ROIs) operation (ENVI 2009).

NDVI

The NDVI is historically one of the first vegetation indices. The processed RED (band 3) and NIR (band 4) images were used to calculate NDVI values for the study area, using the formulation introduced by Rouse et al. (1974):  
formula
(1)

The main concept behind the NDVI is that for vegetated surface, red and near-infrared wavelengths are characterized by high and low absorptions, respectively (Chen et al. 2003). When sunlight strikes objects, certain wavelengths of this spectrum are absorbed and other wavelengths are reflected. The pigment in plant leaves, chlorophyll, strongly absorbs visible light (from 0.4 to 0.7 μm) for use in photosynthesis. The cell structure of the leaves, on the other hand, strongly reflects near-infrared light (from 0.7 to 1.1 μm). The more leaves a plant has, the more these wavelengths of light are affected, respectively. In general, if there is much more reflected radiation in near-infrared wavelengths than in visible wavelengths, then the vegetation in that pixel is likely to be dense and may contain some type of forest. If there is very little difference in the intensity of visible and near-infrared wavelengths reflected, then the vegetation is probably sparse and may consist of grassland, tundra, or desert. Healthy vegetation absorbs most of the visible light that hits it, and reflects a large proportion of the near-infrared light. Unhealthy or sparse vegetation reflects more visible light and less near-infrared light. Nearly all satellite vegetation indices employ this difference formula to quantify the density of plant growth on the Earth. The NDVI value for a specific pixel always varies from −1 to +1. Water bodies are specified by extreme negative values, surfaces with no vegetation cover result in near to zero NDVIs, and the highest density of green leaves is indicated by NDVI value close to +1 (0.8–0.9). NDVI indicates the vegetation cover level and also acts as a beneficial index for monitoring vegetation variations and land cover changes.

The aforementioned steps were applied to obtain the NDVI maps of the study area for all images. The ‘mean NDVI’ value (NDVIw) indicates the average of NDVIs in all pixels of land cover image of watershed. It should be noted that the negative values of the pixels (water bodies) reduce the mean values of NDVI. As these negative values exist in all images, they would affect all images during the whole study period and so the mean values would be comparable.

GUHCR model and calibration performance

In this study, the geomorphologic rainfall–runoff model of GUHCR with a single model parameter is investigated as a semi-distributed conceptual model. There are some other conceptual models such as Zoch and Nash models (Singh 1988), which could be used here, but this model was selected as a semi-distributed model to analyze the relation between the model parameter – which relates to the geomorphologic properties of the watershed – and NDVI.

In the GUHCR model, the watershed is divided into sub-watersheds, and each sub-watershed is represented by a linear reservoir. In this way, each watershed is represented by a sequence of reservoirs distributed according to the watershed geomorphology, resulting in an unequal linear reservoir cascade model with distributed rainfall input. The rainfall input is attributed proportionally between sub-basins on the basis of their areas. The excess rainfall (I(t)) and the outflow (Q(t)) of the prior reservoir for each reservoir is the input and output, respectively. Considering distributed excess rainfall in proportion to sub-watershed areas and applying linear reservoir storage and continuity equations together (Chow et al. 1988) to all the reservoirs, the system equations for N reservoirs are (Singh 1988):  
formula
(2)
 
formula
(3)
in which D is differential operation , ki is the storage coefficient of the ith reservoir, Ci is the ith sub-watershed area and A is the watershed area. Representing the inflow rainfall excess as a Dirac delta function , system equations for the instantaneous unit precipitation are:  
formula
(4)
Equation (4) gives the IUH of each sub-watershed (hi(t)) from which the related outflow hydrograph can be computed. Using the semi-distributed GUHCR model, it is possible to incorporate hydrological conditions of the interior watershed parts.
In the classic model of a cascade of unequal linear reservoirs (Singh 1988), the storage parameter ki (corresponding to lag time) of any sub-watershed, represented by a reservoir, is determined by calibration. However, in the GUHCR model, ki can be related to the geomorphological properties of the sub-watershed and just one unknown parameter as:  
formula
(5)
in which (Saeidifarzad et al. 2014):  
formula
(6)
where is the model parameter with dimension [TL−1/2]. S0i is the average overland slope, Ci is the ith sub-watershed area and L, with dimension [L], is the longest flow path in the drainage network of that sub-watershed. According to Singh (1988), there are several methods for determining the model parameters, such as maximum likelihood, linear and nonlinear programming techniques, but since the moment method uses the characteristics of statistical distribution of the data, like mean and standard deviation values, this method was applied in this study as:  
formula
(7)
where M1 is the first moment of the quantities within parentheses. Equation (7) shows that the model parameter, , is explicitly linked to geomorphological properties of the sub-watersheds. In modeling via GUHCR, only one parameter (i.e., ki) is computed by the empirical equation (Equation (5)) which is related to geomorphological properties (via Ki) of the sub-watersheds.

The root mean square error (RMSE) measure and a related normalization called the Nash–Sutcliffe efficiency (NSE, Nash & Sutcliffe 1970), are the two model performance criteria most widely used for calibration and evaluation of the hydrological models by comparison with observed data (Gupta et al. 2009). Legates & McCabe (1999) indicated that a hydrological model can be sufficiently evaluated by NSE (Equation (8)) and RMSE (Equation (9)), but owing to the importance of peak discharges in flood control and the volume of runoff in water resources management, three other important criteria, the ratio of error of peak flow (EP), the ratio of error of time to peak (tP), and the ratio of error of the hydrograph's volume (Ev) were also used in this study to evaluate the proposed methodology as shown in the appendix (available with the online version of this paper). There are some other criteria that can be used in order to assess the performance of the method. The used criteria are the most widely used measures for calibration and evaluation of the hydrological models by comparison with observed data.

The model parameter was calibrated based on minimizing (maximizing) the RMSE (NSE) and then the performance evaluation was done using the other criteria.

There are several methods for calibration of model parameters such as PEST (Doherty 2007), and meta-heuristic methods like the genetic algorithm (Nourani 2008). In this study, the calibration of the parameter was performed using a semi-automatic trial and error method because there is only one parameter in the model.

ANALYSES AND RESULTS

Using the calculated values of NDVI and the calibrated parameter of the GUHCR model for each event, a regression equation could be fitted to obtain a relationship between NDVI and the parameter of the model for the study area. A schematic diagram of the proposed methodology is shown in Figure 2.

Figure 2

A schematic diagram of the presented methodology showing different steps of modeling.

Figure 2

A schematic diagram of the presented methodology showing different steps of modeling.

At the first step, the hydro-climatologic (rainfall–runoff events and geomorphological characteristics of watersheds) data were gathered for both watersheds and the contemporarily cloud-free land cover images were processed to calculate related NDVI values. Then, hydro-geomorphologic data were applied to calibrate and verify the conceptual rainfall–runoff model (GUHCR). Finally, statistical relationship was created between model parameter and NDVI values.

To develop relationships between the model parameter and temporal trends of NDVI due to the land cover change, the 13 storm events occurring simultaneously over both Lighvan and Hervi sub-watersheds were selected and used along with the corresponding Landsat images. All available 13 events data were used for calibration and verification of the hydrological model to reduce the uncertainty of the calibrated model parameter. However, for detecting the land cover changes (via NDVI), since there are not significant changes in land cover in a short period of each year (within 2, 3 months), only one image (most appropriate image) for each year (totalling seven images, see Table 3) was used to find long-term relationship between NDVI of image and the calibrated parameter of the hydrologic model.

Computed values of NDVI

Figure 3 shows the NDVI maps, derived from the respective Landsat images. NDVI values smaller than 0.1 correspond to no vegetation or the urbanized part of the watersheds, while the higher NDVI values correspond to higher vegetation coverage. In order to investigate seasonal NDVI dynamics, a representative dynamic and stable pixel for each sub-watershed was selected.

Figure 3

NDVI classification map of Landsat images for different dates showing temporal land cover changes.

Figure 3

NDVI classification map of Landsat images for different dates showing temporal land cover changes.

Figure 4 shows the location of pixels, also changes of NDVI in selected pixels for both sub-watersheds. As can be seen in Figure 4, one of the group pixels in both sub-watersheds (Pixel 1) does not change significantly during the seasons while the other pixels which contain dense vegetation cover clearly show a dynamic NDVI behavior. This indicates the seasonal change of land cover. Constant NDVI value for stable pixels demonstrates there are no/limited inconsistencies in derived NDVI maps. Considering both kinds of these pixels, the values of NDVI may have been determined with some degrees of uncertainty, but overall, a decreasing trend can be detected from the NDVI values (see Figure 3 and Table 3). The mean NDVI values of the start and end dates of the study period (1998 and 2012) were respectively computed as 0.283, 0.210 for the Lighvan sub-watershed, and 0.219, 0.172 for the Hervi sub-watershed.

Figure 4

The location of the selected pixels for NDVI analysis and their time series graphs: (a) Hervi and (b) Lighvan sub-watersheds (H and L denote Hervi and Lighvan, respectively).

Figure 4

The location of the selected pixels for NDVI analysis and their time series graphs: (a) Hervi and (b) Lighvan sub-watersheds (H and L denote Hervi and Lighvan, respectively).

The average of mean NDVI values for all images in Hervi and Lighvan sub-watersheds are respectively 0.17 and 0.21. As can be seen clearly from the images (Figure 3), the NDVI values generally decreased during the 14-year study period (as expected) for both sub-watersheds. However, because of inter-annual changes inconsistent behavior is observed. The changes in the Hervi sub-watershed are more considerable, mostly due to land cover changes, with the northern part experiencing a considerable loss of vegetation between 1998 and 2012 (see Figure 3). In general, the residential portion of the Lighvan sub-watershed is smaller and so the NDVI values are larger, with changes being mainly due to a land-use shift from grassland to farmland. The following sub-sections present and discuss the relationships between the model parameter and temporal evolution of NDVI.

Calibration of model parameters

The geomorphological parameters of the GUHCR model (Ci, S0i, Li) were first extracted using GIS tools (Table 4). Then, the parameter was determined via Equations (5)–(7), using the sub-watersheds' physical parameters (Table 4) and data of 10 calibration storm events.

Table 4

Geomorphological parameters of both study sub-watersheds extracted via GIS

Sub-watershed Ci (km2S0i (%) Li (m) Ki (m1/2
Hervi 59.56 14.43 47,071 131.05 
Lighvan 77.15 24.95 17,220 132.89 
Sub-watershed Ci (km2S0i (%) Li (m) Ki (m1/2
Hervi 59.56 14.43 47,071 131.05 
Lighvan 77.15 24.95 17,220 132.89 

Note:Ci, S0i, Li, and Ki are the area, the average overland slope, the longest flow path in the drainage network, and geomorphological parameter of sub-watersheds, respectively.

Table 5 displays the parameters of the GUHCR models during 1998–2012 calibrated for each rainfall–runoff event. The GUHCR model parameter is related to the watershed's geomorphology and the event's characteristics. As it can be seen from Table 5, parameter decreases from 1998 to 2012 in both sub-watersheds, although it has more fluctuations in Lighvan which has smaller urbanized areas (see Figure 3).

Table 5

The calibrated parameter of GUHCR model for Hervi and Lighvan sub-watersheds for different calibration events

Event date Hervi Lighvan 
31/Jul/1998 0.085 0.073 
04/May/1999 0.048 0.055 
08/Apr/2003 0.030 0.030 
11/Jun/2003 0.035 0.035 
22/Apr/2006 0.036 0.032 
28/Jun/2006 0.025 0.025 
07/Apr/2010 0.030 0.041 
18/Jun/2010 0.022 0.021 
22/Apr/2011 0.016 0.032 
18/Jul/2012 0.019 0.024 
Event date Hervi Lighvan 
31/Jul/1998 0.085 0.073 
04/May/1999 0.048 0.055 
08/Apr/2003 0.030 0.030 
11/Jun/2003 0.035 0.035 
22/Apr/2006 0.036 0.032 
28/Jun/2006 0.025 0.025 
07/Apr/2010 0.030 0.041 
18/Jun/2010 0.022 0.021 
22/Apr/2011 0.016 0.032 
18/Jul/2012 0.019 0.024 

Table 6 summarizes the model performance criteria (NSE, RMSE, Ep, tP and Ev) for the calibration events which indicate acceptable performance of the GUHCR model. Based on the Ep, tP and Ev values in Table 6, the GUHCR model simulates the peak flow, time to peak, and the volume of runoff accurately. The results of Table 6 show that the variation of parameter before 2000 is greater than those in the other years for the study period. This may be associated with more rapid increases in urbanization and gardening, and decrease in grassland and farmland. Beginning around the year 2000, the study area was designated as an experimental watershed by East Azerbaijan local water organization, and human activities were restricted.

Table 6

Calculated efficiency criteria in calibration stage using GUHCR model for both study sub-watersheds

Sub-watershed Parameter Date
 
31/Jul/1998 04/May/1999 08/Apr/2003 11/Jun/2003 22/Apr/2006 28/Jun/2006 07/Apr/2010 18/Jun/2010 22/Apr/2011 18/Jun/2012 Average 
Hervi  0.085 0.048 0.030 0.035 0.036 0.025 0.030 0.022 0.016 0.019 0.035 
NSE 0.74 0.86 0.86 0.83 0.60 0.88 0.86 0.74 0.87 0.96 0.82 
RMSE (mm) 0.045 0.026 0.030 0.037 0.019 0.045 0.008 0.026 0.021 0.029 0.03 
Ep (%) −13.98 −6.72 31.08 −10.56 −16.81 −8.73 21.14 18.80 10.98 −8.71 1.65 
tp (%) −6.67 −25.00 −20.00 −20.00 −33.33 0.00 −16.67 −25.00 −25.00 −25.00 −19.67 
EV (%) −20.79 −10.15 2.47 −14.43 −26.15 −16.58 −3.50 −15.01 −16.31 −11.70 −13.22 
Lighvan  0.073 0.055 0.030 0.035 0.032 0.025 0.041 0.021 0.032 0.024 0.037 
NSE 0.71 0.72 0.88 0.83 0.69 0.77 0.85 0.95 0.80 0.75 0.80 
RMSE (mm) 0.05 0.04 0.03 0.03 0.02 0.02 0.01 0.01 0.03 0.07 0.03 
Ep (%) −22.01 −21.35 −20.13 −14.92 −7.71 −25.04 −8.73 6.46 −28.65 −14.74 −15.68 
tp (%) −6.67 −33.33 0.00 −25.00 −28.57 −25.00 −22.22 0.00 −12.50 −16.67 −17.00 
EV (%) −28.17 −24.24 −16.51 −24.61 −14.90 −3.66 −5.79 9.16 −17.48 −28.22 −15.44 
Sub-watershed Parameter Date
 
31/Jul/1998 04/May/1999 08/Apr/2003 11/Jun/2003 22/Apr/2006 28/Jun/2006 07/Apr/2010 18/Jun/2010 22/Apr/2011 18/Jun/2012 Average 
Hervi  0.085 0.048 0.030 0.035 0.036 0.025 0.030 0.022 0.016 0.019 0.035 
NSE 0.74 0.86 0.86 0.83 0.60 0.88 0.86 0.74 0.87 0.96 0.82 
RMSE (mm) 0.045 0.026 0.030 0.037 0.019 0.045 0.008 0.026 0.021 0.029 0.03 
Ep (%) −13.98 −6.72 31.08 −10.56 −16.81 −8.73 21.14 18.80 10.98 −8.71 1.65 
tp (%) −6.67 −25.00 −20.00 −20.00 −33.33 0.00 −16.67 −25.00 −25.00 −25.00 −19.67 
EV (%) −20.79 −10.15 2.47 −14.43 −26.15 −16.58 −3.50 −15.01 −16.31 −11.70 −13.22 
Lighvan  0.073 0.055 0.030 0.035 0.032 0.025 0.041 0.021 0.032 0.024 0.037 
NSE 0.71 0.72 0.88 0.83 0.69 0.77 0.85 0.95 0.80 0.75 0.80 
RMSE (mm) 0.05 0.04 0.03 0.03 0.02 0.02 0.01 0.01 0.03 0.07 0.03 
Ep (%) −22.01 −21.35 −20.13 −14.92 −7.71 −25.04 −8.73 6.46 −28.65 −14.74 −15.68 
tp (%) −6.67 −33.33 0.00 −25.00 −28.57 −25.00 −22.22 0.00 −12.50 −16.67 −17.00 
EV (%) −28.17 −24.24 −16.51 −24.61 −14.90 −3.66 −5.79 9.16 −17.48 −28.22 −15.44 

Relationship between NDVI and IUH parameter

The previously shown results for the study sub-watersheds demonstrate the effect of land cover changes on the results of the rainfall–runoff model. As mentioned before, the changes of the model parameter are somewhat similar to NDVI changes in the study period.

As can be seen in Tables 3 and 5, both and NDVI values show a decreasing trend from 1998 to 2012. NDVI change in time mainly reflects the impact of human activities, urbanization, and conversion from sparse to dense pasture, as well as seasonal changes. Parameter of the model for both sub-watersheds can therefore be related to the NDVI values through a regression analysis.

Figure 5 illustrates the relationship between the parameter and the decreasing percent of average NDVI values (Indvi) during 1998–2012 for Hervi and Lighvan sub-watersheds, respectively. It should be mentioned that the parameter is dependent on each rainfall–runoff event; also, the impact of climate condition is implicitly involved in both rainfall and runoff behaviors and data.

Figure 5

The variation of parameter under different levels of NDVI change in (a) Hervi and (b) Lighvan sub-watersheds; fitted regression and fitting efficiency are also shown for each study sub-watershed.

Figure 5

The variation of parameter under different levels of NDVI change in (a) Hervi and (b) Lighvan sub-watersheds; fitted regression and fitting efficiency are also shown for each study sub-watershed.

The NDVI values of year 1985 for each sub-watershed (Figure 3, NDVI values are 0.231 and 0.286 for Hervi and Lighvan sub-watersheds, respectively) was considered as the reference NDVI and the percent of NDVI decreasing (Indvi) during 1998–2012 was calculated with regard to this reference value. While a power form of the regression equation could be best fitted to the data on parameter using available images of Hervi sub-watershed, for the Lighvan sub-watershed the fitted equation shows somewhat lower correlation between the parameters due to more fluctuations of the parameter. This might be due to the more limited land cover changes in the Lighvan sub-watershed during 1998–2012. One of the most important features of land cover change can be urbanization. This means that the sub-watershed has been mostly covered by vegetation cover with highly temporal fluctuations (from each season and one year to the next) during the study period.

To evaluate the performance of the fitted equations to estimate the model parameter, the three rainfall–runoff events which were not used in the calibration stage were examined to verify the rainfall–runoff modeling for both sub-watersheds. For the GUHCR model, which consists of only one parameter of , the parameter was calculated for each sub-watershed using the NDVI value of the verification events and the fitted power equation (see Figure 5) as:  
formula
(8)
 
formula
(9)
for Hervi (Equation (8)) and Lighvan (Equation (9)) sub-watersheds.

Table 7 presents the verification results of GUHCR model in the Hervi sub-watershed. As is clear in Table 7, the verification results of the GUHCR model show acceptable accuracy in predicting the peak flow, time to peak, and the volume of the hydrograph.

Table 7

Calculated efficiency criteria for verification events using the GUHCR model for both study sub-watersheds

Sub-watershed Parameter Date
 
Average 
27/May/2003 22/May/2006 28/Jun/2012 
 NDVI 0.16 0.19 0.17  
Hervi  0.023 0.032 0.027 0.027 
NSE 0.82 0.88 0.90 0.87 
RMSE (mm) 0.030 0.022 0.019 0.024 
Ep (%) −3.21 23.44 −12.46 2.59 
tp (%) −25 0.00 −20.00 −15.00 
EV (%) −20.33 −5.88 −11.74 −12.65 
 NDVI 0.12 0.26 0.21  
Lighvan  0.027 0.039 0.031 0.033 
NSE 0.69 0.58 0.38 0.55 
RMSE (mm) 0.025 0.052 0.017 0.031 
Ep (%) −11.30 −40.50 −52.14 −34.65 
tp (%) −50.00 −33.33 −75.00 −52.78 
EV (%) −14.86 −28.48 −21.05 −21.46 
Sub-watershed Parameter Date
 
Average 
27/May/2003 22/May/2006 28/Jun/2012 
 NDVI 0.16 0.19 0.17  
Hervi  0.023 0.032 0.027 0.027 
NSE 0.82 0.88 0.90 0.87 
RMSE (mm) 0.030 0.022 0.019 0.024 
Ep (%) −3.21 23.44 −12.46 2.59 
tp (%) −25 0.00 −20.00 −15.00 
EV (%) −20.33 −5.88 −11.74 −12.65 
 NDVI 0.12 0.26 0.21  
Lighvan  0.027 0.039 0.031 0.033 
NSE 0.69 0.58 0.38 0.55 
RMSE (mm) 0.025 0.052 0.017 0.031 
Ep (%) −11.30 −40.50 −52.14 −34.65 
tp (%) −50.00 −33.33 −75.00 −52.78 
EV (%) −14.86 −28.48 −21.05 −21.46 

Figure 6 compares the observed and simulated hydrographs of the verification events using the GUHCR model for both sub-watersheds. As can be seen from Table 7 and Figure 6, although the output hydrographs of Hervi sub-watershed are accurately predicted, hydrographs of Lighvan sub-watershed are not well simulated. This may be due to temporal oscillation of the NDVI values of the Lighvan sub-watershed which leads to a weak correlation between the model's parameter and the NDVI values. Hervi sub-watershed shows a clear trend during the study period; however, Lighvan does not show a certain trend in vegetation cover and thus the parameter.

Figure 6

Comparison of simulated and observed hydrographs in Hervi and Lighvan sub-watersheds for verification events.

Figure 6

Comparison of simulated and observed hydrographs in Hervi and Lighvan sub-watersheds for verification events.

Based on the calibration and verification results, it can be seen that the values of the criteria are reasonable for Hervi sub-watershed, however not so satisfying for the Lighvan sub-watershed. The average values for NSE, RMSE, Ep, tP, and Ev at the Hervi sub-watershed using the GUHCR model are 0.87, 0.024, 2.59, −15.00, and −12.65, respectively. In contrast, the averages of the aforementioned criteria for Lighvan sub-watershed are 0.55, 0.031, −34.65, −52.78, and −21.46, respectively.

DISCUSSION

Computed values of NDVI

In the current study, it was attempted to examine a proposed geomorphology-based model and remote sensing data to relate the land cover properties of the watershed to the model parameter. This would lead to a more reliable result due to the precise detection of the land use/cover conditions, considering seasonality and permanent land cover changes. Applying this methodology, the model's parameter is determined using the temporal NDVI values. A clear trend of NDVI can be detected in Hervi sub-watershed in the study period. However, in Lighvan sub-watershed, the NDVI values are fluctuate significantly because of a smaller urbanized area and larger natural land cover, and could be more related to climate parameters (seasonality) than anthropogenic effects.

Calibration of model parameters

Figure 6 shows a good agreement between simulated and observed hydrographs in the Hervi sub-watershed. It is also worth noting that the GUHCR model is able to simulate both the rising limbs of hydrographs (which are mostly related to storm properties) and the recession limbs (which are usually related to the watershed morphology) appropriately. As indicated by Ep, tP, and Ev in Table 7, the error of peak flow, time to peak, and the volume of hydrographs show acceptable accuracy for the Hervi sub-watershed. However, the model performance is acceptable in calibration events in both watersheds; the performance criteria are not good in the Lighvan sub-watershed due to the use of regression equation which is poor in Lighvan.

Note that the GUHCR model has an acceptable performance in both sub-watersheds, which is attributable to the use of geomorphologic properties for calculating the model parameter. This model is dependent on physical properties of the watersheds which have less uncertainty in comparison to the historical rainfall–runoff event data (Nourani et al. 2015). Such geomorphological models, which use the physical properties of the watershed to estimate the model parameters, can be employed in ungauged watersheds lacking sufficient data due to less dependency of the model on rainfall–runoff event data. It should be noted that most rainfall–runoff events over the watershed are occurring only in a few months of the year which reduces the uncertainty of the calibrated parameter.

In the GUHCR model, there is only one parameter which is calculated using the geomorphologic properties of the watershed. The model parameter is related to area, slope, and length of the watershed; however, the other geomorphologic parameters can be taken into account in the model. This calibrated parameter is related to land cover variation via a power equation using remote sensing data and can be estimated using the NDVI values. A single power equation has been fitted to create a relation between the NDVI and model parameter; however, having more data for fitting would lead to a more precise model or even to developing an artificial intelligence-based model.

Relationship between NDVI and IUH parameter

Based on the result of the verification step, higher performance criteria are clearly due to the higher correlation between the NDVI and model parameter in Hervi sub-watershed. The high fluctuation values of NDVI in Lighvan sub-watershed might be more related to climate parameters than human effects. Also, this might be due to less uncertainty in watersheds with urbanized parts where runoff to rainfall ratio is larger than in the natural watershed. For small runoff to rainfall ratios in the natural sub-watershed, the uncertainties associated with rainfall measurement errors and spatial variability can form a dominant part of the overall model prediction uncertainty (high noise to signal ratio). In addition, the parameters' values (i.e., physical and land cover properties) in the urbanized part of a watershed are approximately constant while these values are more variable over the season in natural watersheds. This characteristic would show less uncertainty in the urbanized watershed which leads to a more precise result at Hervi sub-watershed than that of the Lighvan sub-watershed. Here, the purpose was to find an overall trend for the model parameter; however, because of fluctuations in the model parameter around the mean value, the fitted equation could not show an acceptable correlation.

In former studies, the parameters of the Nash model were related to urbanization using only the factors of the impervious area and population density via a power equation (Cheng & Wang 2002; Huang et al. 2008, 2012; Cheng et al. 2010). In the mentioned studies, the population was considered as one of the urbanization factors. In countries with population control programs the urbanization development is not related to the population and therefore this factor would not be an indicator of urbanization in such countries. Table 8 presents a summary of the relevant literature on the effect of land cover changes in rainfall–runoff modeling showing the present study in the last row. As can be seen in Table 8, the effects of land cover changes in conceptual models have been calibrated only based on the event data. The invention of GIS also considering the geomorphologic properties of the watersheds has led to the use of semi-distributed rainfall–runoff models. Although many studies have studied the effect of land cover properties on rainfall–runoff models, there are few studies that have investigated the temporal land cover properties of a watershed using the NDVI and the parameter of a semi-distributed model. As mentioned before, in the present study we tried to relate the changes of NDVI to the changes of the model parameter to address the gap of previous studies (see Table 8). Using five criteria to evaluate the investigated model, it could be demonstrated that the obtained regression equation provides a promising approach for estimating lumped model parameters using remote sensing data for one of the study sub-watersheds.

Table 8

Summary of relevant rainfall–runoff modeling studies considering land cover changes (the current paper is added for completeness)

Study Data/Rainfall–runoff model Location Key results 
Kang et al. (1998)  Gauge data/linear reservoir model (Nash model) On-Cheon watershed, Pusan, Korea (1) The peak and time to peak values increased and decreased respectively due to urbanization. (2) The lag time is decreased due to urbanization. (3) N/A 
Cheng & Wang (2002)  Gauge data, impervious area percentage/linear reservoir model (Nash model) Wu-Tu watershed, Taiwan (1) Three decades of urbanization had increased the peak flow by 27%, and the time to peak was decreased by 4 hr. (2) The parameter n (number of reservoirs) decreased with the increase in impervious area. (3) N/A 
Cheng et al. (2010)  Gauge data, population intensity, impervious area percentage/linear reservoir model (Nash model) Wu-Tu watershed, Taiwan (1) N/A. (2) Parameter n responds more sensitively than parameter k to increasing impervious areas and population densities. Additionally, parameter n responds more strongly to imperviousness than to population. (3) N/A 
Ali et al. (2011)  Gauge data/HEC-HMS model Lai Nullah basin, Islamabad, Pakistan (1) Using future land use scenario and calibrating the HEC-HMS rainfall–runoff model, their results showed a good consistency between the simulated and measured hydrographs at the outlet of the basin. (2) N/A. (3) N/A 
Li et al. (2012)  Gauge data, AVHRR-LAI data/Xinanjiang-ET and SIMHYD-ET rainfall–runoff models Crawford River catchment, Victoria, Australia (1) Increase in plantations can reduce streamflow substantially, even more than climate variability. (2) N/A. (3) N/A 
Huang et al. (2012)  Gauge data, population intensity, impervious area percentage/linear reservoir model (Nash model) Wu-Tu watershed, Taiwan (1) The time to peak is diminished as a result of increasing impervious coverings while the peak flow increases. (2) The parameter n (number of reservoirs) decreased with the increase in impervious area. (3) N/A 
Verbeiren et al. (2013)  Gauge data, Landsat and SPOT images/physically based rainfall–runoff model (WetSpa) Tolka River basin, Dublin, Ireland (1) The steady urban growth had a considerable impact on peak discharges. The hydrological response is quicker as a result of urbanization. (2) N/A. (3) N/A 
Miller et al. (2014)  Gauge data, topographic map/catchment hydrological cycle assessment tool (CAT) Haydon Wick brook and Rodbourne catchments, Swindon, UK (1) Reduction in the Muskingum routing parameter k, and increase in peak flow are observed. Reduction in flood duration and response time of a catchment is greatest at low levels of urbanization. (2) N/A. (3) N/A 
Nourani et al. (2015)  Gauge data, Landsat images/geomorphological rainfall–runoff model (GCUR) La Terraza, watershed, Arizona, USA (1) The interior hydrographs in different land use sub-watersheds can be determined by using the index related to land cover. (2) The model parameter has a significant change in different land use sub-watersheds. (3) N/A 
Tesemma et al. (2015b)  Gauge data, GLASS images/VIC hydrological model Goulburn–Broken catchment, Australia (1) Incorporating climate-induced changes in LAI in the VIC model reduced the projected declines in streamflow and confirmed the importance of including the effects of changes in LAI in future projections of streamflow. (2) N/A. (3) N/A 
This study Gauge data, Landsat images/geomorphological rainfall–runoff model (GUHCR) Hervi and Lighvan watersheds, Tabriz, Iran (1) The precise detection of the land use/cover conditions by remotely sensed data would lead to more reliable results. (2) The change of the model parameter is somehow similar to NDVI changes in the study period. (3) A power form of the regression equation could be fitted to the data on parameters of the model and NDVI values using available images 
Study Data/Rainfall–runoff model Location Key results 
Kang et al. (1998)  Gauge data/linear reservoir model (Nash model) On-Cheon watershed, Pusan, Korea (1) The peak and time to peak values increased and decreased respectively due to urbanization. (2) The lag time is decreased due to urbanization. (3) N/A 
Cheng & Wang (2002)  Gauge data, impervious area percentage/linear reservoir model (Nash model) Wu-Tu watershed, Taiwan (1) Three decades of urbanization had increased the peak flow by 27%, and the time to peak was decreased by 4 hr. (2) The parameter n (number of reservoirs) decreased with the increase in impervious area. (3) N/A 
Cheng et al. (2010)  Gauge data, population intensity, impervious area percentage/linear reservoir model (Nash model) Wu-Tu watershed, Taiwan (1) N/A. (2) Parameter n responds more sensitively than parameter k to increasing impervious areas and population densities. Additionally, parameter n responds more strongly to imperviousness than to population. (3) N/A 
Ali et al. (2011)  Gauge data/HEC-HMS model Lai Nullah basin, Islamabad, Pakistan (1) Using future land use scenario and calibrating the HEC-HMS rainfall–runoff model, their results showed a good consistency between the simulated and measured hydrographs at the outlet of the basin. (2) N/A. (3) N/A 
Li et al. (2012)  Gauge data, AVHRR-LAI data/Xinanjiang-ET and SIMHYD-ET rainfall–runoff models Crawford River catchment, Victoria, Australia (1) Increase in plantations can reduce streamflow substantially, even more than climate variability. (2) N/A. (3) N/A 
Huang et al. (2012)  Gauge data, population intensity, impervious area percentage/linear reservoir model (Nash model) Wu-Tu watershed, Taiwan (1) The time to peak is diminished as a result of increasing impervious coverings while the peak flow increases. (2) The parameter n (number of reservoirs) decreased with the increase in impervious area. (3) N/A 
Verbeiren et al. (2013)  Gauge data, Landsat and SPOT images/physically based rainfall–runoff model (WetSpa) Tolka River basin, Dublin, Ireland (1) The steady urban growth had a considerable impact on peak discharges. The hydrological response is quicker as a result of urbanization. (2) N/A. (3) N/A 
Miller et al. (2014)  Gauge data, topographic map/catchment hydrological cycle assessment tool (CAT) Haydon Wick brook and Rodbourne catchments, Swindon, UK (1) Reduction in the Muskingum routing parameter k, and increase in peak flow are observed. Reduction in flood duration and response time of a catchment is greatest at low levels of urbanization. (2) N/A. (3) N/A 
Nourani et al. (2015)  Gauge data, Landsat images/geomorphological rainfall–runoff model (GCUR) La Terraza, watershed, Arizona, USA (1) The interior hydrographs in different land use sub-watersheds can be determined by using the index related to land cover. (2) The model parameter has a significant change in different land use sub-watersheds. (3) N/A 
Tesemma et al. (2015b)  Gauge data, GLASS images/VIC hydrological model Goulburn–Broken catchment, Australia (1) Incorporating climate-induced changes in LAI in the VIC model reduced the projected declines in streamflow and confirmed the importance of including the effects of changes in LAI in future projections of streamflow. (2) N/A. (3) N/A 
This study Gauge data, Landsat images/geomorphological rainfall–runoff model (GUHCR) Hervi and Lighvan watersheds, Tabriz, Iran (1) The precise detection of the land use/cover conditions by remotely sensed data would lead to more reliable results. (2) The change of the model parameter is somehow similar to NDVI changes in the study period. (3) A power form of the regression equation could be fitted to the data on parameters of the model and NDVI values using available images 

Not all papers performed the effect of land cover properties on the model parameters and the temporal changes of the parameters. These are denoted with N/A representing ‘not applicable’ in the relevant part of the Key results column. In the Key results column, the above-mentioned three components are identified by the code: (1) effect of land cover changes on runoff; (2) effect of land cover changes on the model parameters; and (3) investigation of temporal vegetation cover index.

CONCLUSION

In this study, the remotely sensed data were used to analyze the NDVI values during 1998–2012. Considering a reference value for the NDVI, the percent of NDVI decreasing (Indvi) during the study period was calculated, which was used to determine a relationship between the model parameter and NDVI. Regression relationships were established between temporal variations in remotely sensed NDVI and the parameter of the geomorphologic IUH model. Thereby, a mechanism is provided to monitor and detect the temporal variation of watershed land cover and its consequent impacts on the rainfall–runoff response of the watershed. A single power equation was fitted on the NDVI and the model parameter values. However, having more data would contribute to a more precise model.

The obtained results indicate that decreasing NDVI percentages (e.g., due to urbanization or other land cover/use changes, such as seasonal variation) may show a considerable impact on values of the model parameter. While a power form of the regression equation could be best fitted to the data on parameter using available images of Hervi sub-watershed, the fitted equation for the Lighvan sub-watershed implies slightly lower correlation between the parameters due to more fluctuations of the parameter. As a result, despite the accurate prediction of the output hydrographs of Hervi sub-watershed, hydrographs of Lighvan sub-watershed are not well simulated. This might be the consequence of temporal oscillation of the NDVI values in the Lighvan sub-watershed which leads to a weak correlation between the model's parameter and the NDVI values. To evaluate the performance of the fitted equations, three rainfall–runoff events were examined to verify the rainfall–runoff modeling for both sub-watersheds. The verification results of the GUHCR model show acceptable accuracy in predicting the peak flow, time to peak, and the volume of the hydrograph. The average NSE values for the GUHCR model in the presented study are 0.87 and 0.55 in Hervi and Lighvan sub-watersheds, respectively.

The NDVI analyses of the ‘stable pixels’ of both sub-watersheds show that combination of seasonal and inter-annual changes of NDVI would lead to an overall downwarding NDVI which, in turn, can affect the rainfall–runoff model parameters.

Our results confirmed that the single parameter of the GUHCR model is related to the geomorphologic properties of the watershed, which has a physical meaning. Employing the physical properties of the watershed to calculate the model parameters, such geomorphological models can be efficiently applied in ungauged watersheds lacking sufficient data due to less dependency of the model on rainfall–runoff event data.

As a future research plan, it can be suggested to employ other conceptual IUHs and geomorphologic IUHs or fully distributed models using the proposed methodology. Also, the effect of temporal and seasonal changes of NDVI on models parameters can be analyzed more carefully considering the seasonal variations in type and phenology of vegetation cover. It can be suggested to obtain a combined description by different linear equations which would have better performance. Implementation of other remote sensing indices with different characteristics and abilities (such as LAI, instead of NDVI) to examine the land cover change effects on the performance of the hydrologic model can be considered as a research plan for future studies. The results of different climatic watersheds can be investigated and compared with each other to obtain a global relationship.

The applied modeling and method have some limitations just like other models. The modeling was performed assuming the watershed as a linear system which is nonlinear in reality. This assumption might result in some errors in the model's output. Furthermore, the number of rainfall–runoff events was limited to some seasons, which might be a limitation in determining the parameters' trend; however, the model can be updated to have access to more data in the future.

ACKNOWLEDGEMENTS

This paper was supported by a research grant of the University of Tabriz. We thank East Azerbaijan local water organization in Tabriz, Iran, that provided the events data related to the study area (to obtain the data contact pubrel@azarwater.ir). The satellite images of the studied area were downloaded from the USGS website (https://earthexplorer.usgs.gov/login/).

REFERENCES

REFERENCES
Ali
,
M.
,
Khan
,
S. J.
,
Aslam
,
I.
&
Khan
,
Z.
2011
Simulation of the impacts of land-use change on surface runoff of Lai Nullah Basin in Islamabad, Pakistan
.
Landscape and Urban Planning
102
(
4
),
271
279
.
Baret
,
F.
&
Guyot
,
G.
1991
Potentials and limits of vegetation indices for LAI and APAR assessment
.
Remote Sensing of Environment
35
(
2–3
),
161
173
.
Beck
,
P. S. A.
,
Atzberger
,
C.
&
Høgda
,
K. A.
2006
Improved monitoring of vegetation dynamics at very high latitudes, a new method using MODIS NDVI
.
Remote Sensing of Environment
100
,
321
336
.
Boyd
,
M. J.
,
Bufill
,
M. C.
&
Knee
,
R. M.
1993
Pervious and impervious runoff in urban catchments
.
Hydrological Sciences Journal
38
(
6
),
463
478
.
Braud
,
I.
,
Breil
,
P.
,
Thollet
,
F.
,
Lagouy
,
M.
,
Branger
,
F.
,
Jacqueminet
,
C.
,
Kermadi
,
S.
&
Michel
,
K.
2013
Evidence of the impact of urbanization on the hydrological regime of a medium-sized periurban catchment in France
.
Journal of Hydrology
485
,
5
23
.
Chen
,
P. Y.
,
Srinivasan
,
R.
,
Fedosejevs
,
G.
&
Kiniry
,
J. R.
2003
Evaluating different NDVI composite techniques using NOAA-14 AVHRR data
.
International Journal of Remote Sensing
24
(
17
),
3403
3412
.
Cheng
,
S. J.
,
Lee
,
C. F.
&
Lee
,
J. H.
2010
Effects of urbanization factors on model parameters
.
Water Resources Management
24
,
775
794
.
Chow
,
V. T.
,
Maidment
,
D. R.
&
Mays
,
L. W.
1988
Applied Hydrology
.
McGraw-Hill
,
New York
,
USA
.
Doherty
,
J.,
2007
FORTRAN 90 Modules for Implementation of Parallelised, Model-Independent, Model-Based Processing
.
Watermark Numerical Computing
,
Australia
.
Eckert
,
S.
,
Hüsler
,
F.
,
Liniger
,
H.
&
Hodel
,
E.
2015
Trend analysis of MODIS NDVI time series for detecting land degradation and regeneration in Mongolia
.
Journal of Arid Environments
113
,
16
28
.
ENVI
2009
ENVI User's Guide, Version 4.7
.
ITT Visual Information Solutions
.
Goswami
,
S.
,
Gamon
,
J.
,
Vargas
,
S.
&
Tweedie
,
C.
2015
Relationships of NDVI, biomass, and Leaf Area Index (LAI) for six key plant species in Barrow, Alaska
.
PeerJ PrePrints
3
,
e1127
.
https://doi.org/10.7287/peerj.preprints.913v1
.
Gupta
,
H. V.
,
Kling
,
H.
,
Yilmaz
,
K. K.
&
Martinez
,
G. F.
2009
Decomposition of the mean squared error and NSE performance criteria: implications for improving hydrological modeling
.
Journal of Hydrology
377
,
80
91
.
Huang
,
H. J.
,
Cheng
,
S. J.
,
Wen
,
J. C.
&
Lee
,
J. H.
2008
Effect of growing watershed imperviousness on hydrograph parameters and peak discharge
.
Hydrological Processes
22
,
2075
2085
.
Huang
,
S. Y.
,
Cheng
,
S. J.
,
Wen
,
J. C.
&
Lee
,
J. H.
2012
Identifying hydrograph parameters and their relationships to urbanization variables
.
Hydrological Sciences Journal
57
(
1
),
144
161
.
Jiang
,
W.
,
Yuan
,
L.
,
Wang
,
W.
,
Cao
,
R.
,
Zhang
,
Y.
&
Shen
,
W.
2015
Spatio-temporal analysis of vegetation variation in the Yellow River basin
.
Ecological Indicators
51
,
117
126
.
Kang
,
I. S.
,
Park
,
J. I.
&
Singh
,
V. P.
1998
Effect of urbanization on runoff characteristics of the On-Cheon Stream watershed in Pusan, Korea
.
Hydrological Processes
12
,
351
363
.
Kepner
,
W. G.
,
Semmens
,
D. J.
,
Hernandez
,
M.
&
Goodrich
,
D. C.
2008
Evaluating hydrological response to forecasted land use change
.
Chapter 15
. In:
The North America Land Cover Summit
(J. C. Campbell, K. Bruce Jones, J. H. Smith, M. K. Knoppe, eds)
.
Association of American Geographers
,
Washington, DC
,
USA
, pp.
275
292
.
Legates
,
D. R.
&
McCabe
,
G. J.
Jr
1999
Evaluating the use of ‘goodness-of-fit’ measures in hydrologic and hydro-climatic model validation
.
Water Resources Research
35
.
doi:10.1029/1998WR900018
.
Li
,
H. X.
,
Zhang
,
Y. Q.
,
Chiew
,
F. H. S.
&
Xu
,
S. G.
2009
Predicting runoff in ungauged catchments by using Xinanjiang model with MODIS leaf area index
.
Journal of Hydrology
370
(
1–4
),
155
162
.
Liu
,
X.
,
Ren
,
L.
,
Yuan
,
F.
,
Xu
,
J.
&
Liu
,
W.
2012
Assessing vegetation response to drought in the Laohahe catchment North China
.
Hydrology Research
43
(
1–2
),
91
101
.
Lu
,
D.
,
Ge
,
H.
,
He
,
S.
,
Xu
,
A.
,
Zhou
,
G.
&
Du
,
H.
2008
Pixel-based Minnaert correction method for reducing topographic effects on a Landsat 7 ETM+ image
.
Photographic Engineering and Remote Sensing
74
(
11
),
1343
1350
.
Miller
,
S. N.
,
Semmens
,
D.
,
Goodrich
,
D.
,
Hernandez
,
M.
,
Miller
,
R.
,
Kepner
,
W.
&
Guertin
,
D. P.
2007
The automated geospatial watershed assessment tool
.
Environmental Modelling & Software
22
,
365
377
.
Miller
,
J. D.
,
Kim
,
H.
,
Kjeldsen
,
T. R.
,
Packman
,
J.
,
Grebby
,
S.
&
Dearden
,
R.
2014
Assessing the impact of urbanization on storm runoff in a peri-urban catchment using historical change in impervious cover
.
Journal of Hydrology
515
,
59
70
.
Nash
,
J. E.
1957
The form of the instantaneous unit hydrograph
.
Hydrological Sciences Bulletin
3
,
114
121
.
Nash
,
J. E.
&
Sutcliffe
,
J. V.
1970
River flow forecasting through conceptual models I: a discussion of principles
.
Journal of Hydrology
10
(
3
),
282
290
.
Nourani
,
V.
,
Fakheri Fard
,
A.
,
Niazi
,
F.
,
Gupta
,
H. V.
,
Goodrich
,
D. C.
&
Valizadeh Kamran
,
K.
2015
Implication of remotely sensed data to incorporate land cover effect into a linear reservoir-based rainfall-runoff model
.
Journal of Hydrology
529
,
94
105
.
Rouse
,
J. W.
,
Haas
,
R. H.
,
Schell
,
J. A.
,
Deering
,
D. W.
&
Harlan
,
J. C.
1974
Monitoring the vernal advancements and retrogradation (green wave effect) of nature vegetation
.
NASA/GSFC Final Report
,
Greenbelt, MD
,
USA
.
Saeidifarzad
,
B.
,
Nourani
,
V.
,
Aalami
,
M. T.
&
Chau
,
K.-W.
2014
Multi-site calibration of linear reservoir based geomorphologic rainfall-runoff models
.
Water
6
,
2690
2716
.
Shuster
,
W.
,
Bonta
,
J.
,
Thurston
,
H.
,
Warnemuende
,
E.
&
Smith
,
D.
2005
Impacts of impervious surface on watershed hydrology: a review
.
Urban Water Journal
2
(
4
),
263
275
.
Singh
,
V. P.
1988
Hydrologic Systems. Rainfall–Runoff Modeling
,
Vol. I
.
Prentice-Hall
,
Englewood Cliffs, NJ
,
USA
.
Tian
,
F.
,
Feng
,
X.
,
Zhang
,
L.
,
Fu
,
B.
,
Wang
,
S.
,
Lv
,
Y.
&
Wang
,
P.
2017
Effects of revegetation on soil moisture under different precipitation gradients in the Loess Plateau, China
.
Hydrology Research
48
(
5
),
1378
1390
.
doi:10.2166/nh.2016.022
.
Törnros
,
T.
&
Menzel
,
L.
2014
Leaf area index as a function of precipitation within a hydrological model
.
Hydrology Research
45
(
4–5
),
660
672
.
Turner
,
D. P.
,
Cohen
,
W. B.
,
Kennedy
,
R. E.
,
Fassnacht
,
K. S.
&
Briggs
,
J. M.
1999
Relationships between Leaf Area Index and landsat TM spectral vegetation indices across three temperate zone sites
.
Remote Sensing of Environment
70
,
52
68
.
Verbeiren
,
B.
,
Van De Voorde
,
T.
,
Canters
,
F.
,
Binard
,
M.
,
Cornet
,
Y.
&
Batelaan
,
O.
2013
Assessing urbanisation effects on rainfall-runoff using a remote sensing supported modelling strategy
.
International Journal of Applied Earth Observation and Geoinformation
21
,
92
102
.
Wang
,
J.
,
Rich
,
P. M.
&
Price
,
K. P.
2003
Temporal responses of NDVI to precipitation and temperature in the central Great Plains, USA
.
International Journal of Remote Sensing
24
,
2345
2364
.
Zhang
,
Y. Q.
,
Chiew
,
F. H. S.
,
Zhang
,
L.
&
Li
,
H. X.
2009
Use of remotely sensed actual evapotranspiration to improve rainfall-runoff modeling in Southeast Australia
.
Journal of Hydrometeorology
10
(
4
),
969
980
.
Zhiqiang
,
T.
,
Qi
,
Z.
,
Mengfan
,
L.
,
Yunliang
,
L.
,
Xiuli
,
X.
&
Jiahu
,
J.
2016
A study of the relationship between wetland vegetation communities and water regimes using a combined remote sensing and hydraulic modeling approach
.
Hydrology Research
47
(
S1
),
278
292
.
doi:10.2166/nh.2016.216
.
Zhou
,
Y. C.
,
Zhang
,
Y. Q.
,
Vaze
,
J.
,
Lane
,
P.
&
Xu
,
S. G.
2013
Improving runoff estimates using remote sensing vegetation data for bushfire impacted catchments
.
Agricultural and Forest Meteorology
182
,
332
341
.

Supplementary data