GIS and remote sensing techniques for the assessment of land use change impact on flood hydrology : the case study of Yialias basin in Cyprus

Floods are one of the most common natural disasters worldwide, leading to economic losses and loss of human lives. This paper highlights the hydrological effects of multi-temporal land use changes in flood hazard within the Yialias catchment area, located in central Cyprus. A calibrated hydrological model was firstly developed to describe the hydrological processes and internal basin dynamics of the three major subbasins, in order to study the diachronic effects of land use changes. For the implementation of the hydrological model, land use, soil and hydrometeorological data were incorporated. The climatic and stream flow data were derived from rain and flow gauge stations located in the wider area of the watershed basin. In addition, the land use and soil data were extracted after the application of object-oriented nearest neighbor algorithms of ASTER satellite images. Subsequently, the cellular automata (CA)–Markov chain analysis was implemented to predict the 2020 land use/land cover (LULC) map and incorporate it to the hydrological impact assessment. The results denoted the increase of runoff in the catchment area due to the recorded extensive urban sprawl phenomenon of the last decade.


Introduction
Land use and floods are closely related; therefore, any changes in the land use, such as urbanization across the catchment's area, may trigger a sequence of flood occurrences (Hadjimitsis, 2010).The current and future development in water resources is very sensitive to land use and intensification of human activities.It is expected that flood risk will continue to rise, as a consequence of a combination of climate change (e.g., Kundzewicz et al., 2005;Tsanis et al., 2011;Grillakis et al., 2011) and an increase in exposure vulnerability (e.g., due to increasing flood plain occupancy), increase in endangered areas and changes in the terrestrial system (e.g., land cover changes and river regulation; see Elmer et al., 2012).Human transformation of the Earth's land surface seems to have multiple consequences for biophysical systems at all scales (Roosmalen et al., 2009).During the past decades, airborne and spaceborne remote sensing technologies along with geographical information systems (GISs) have been widely used for flood monitoring, including flash floods (Taubenbock et al., 2011).
Flash floods respond to the causative storms in a short period of time, with water levels in the drainage network reaching peak levels within a few minutes or hours, allowing for a very limited time window for warnings to be prepared and issued (Koutroulis and Tsanis, 2010;Grillakis et al., 2010).Modeling of floods has greatly improved in recent years, with the advent of GIS, satellite remote sensing imagery, highresolution digital elevation models (DEMs), distributed hydrologic models, and development of real time flood forecasting and delivery systems on the internet (Garrote and Bras, 1995;Bedient et al., 2003).Hydrological and hydraulic simulation models are essential tools to evaluate potential consequences of proposed strategies and to facilitate management decisions.Nowadays satellite remote sensing has the potential to provide extensive coverage of key variables such as precipitation and soil moisture as well as many of the parameters such as vegetation cover, vegetation change and imperviousness that are important inputs to modern hydrological models (De Fries and Eshelman, 2004).According to Mao and Cherkauer (2012), human activity is one of the major driving forces leading to changes in land cover characteristics and subsequently hydrologic processes.Land use influences the infiltration and soil water distribution process because saturated hydraulic conductivity is influenced by plant roots and pores resulting from the presence of soil fauna (Ragab and Cooper, 1993;Fohrer et al., 2000).A characteristic example is the influence of buildup areas and roads on overland flow, flood frequency and magnitude (Nejadhashemi et al., 2011).Therefore, land cover plays a key role in controlling the hydrologic regime of a catchment area through a number of different parameters such as leaf area index, evapotranspiration, soil moisture content and infiltration capacity, surface and subsurface flow regimes including base-flow contributions to streams and recharge, surface roughness, runoff as well as soil erosion through complex interactions among vegetation, soils, geology, terrain and climate processes.
Especially urban areas are prone to flooding due to the large proportion of impermeable surface cover such as concrete that increases the total volume of runoff and peak flows and shortens the time that the floodwaters take to arrive at peak runoff (Hall, 1984;Knebl, 2005).In various studies, historical and present land use/land cover patterns or extreme scenarios have been used as input in hydrologic models to determine hydrologic responses to different scenarios in a combined integrated approach (Moiwo et al., 2010;Hong et al., 2010).At different watershed scales, several researchers (Savary et al., 2009;Schilling et al., 2010;Turnbull et al., 2012) have developed various methods aiming at quantifying the hydrologic alterations in relation to land cover change.
This paper attempts to quantify the sensitivity of the distributed hydrological model to the land use and soil parameterizations, in order to simulate runoff processes in a catchment area in Cyprus, namely Yialias watershed.Specifically, the potential use of remote sensing in providing hydrological models with adequate, reliable and updated land use data is highlighted.The major flood event that occurred between 12 and 13 February 2003 was successively simulated with the use of multi-temporal land use data of the specific period (data of 2000) and data of 2010 (keeping the same meteorological parameters).
The aim of this approach was to assess the impact of land use changes (including conversions between different land use types and shifts in the geographic extent of those land use types) to the runoff processes and hydrologic response.In the following, a CA-Markov chain analysis was implemented to calculate and predict the area's land use/land cover (LULC) regime for 2020 and incorporate it to the hydrological model for assessing watershed's basin response.The hydrological model used is the HEC-HMS in distributed mode to utilize the modified Clark (Clark, 1945) method of transforming the excess rainfall to runoff.Moreover, the USDA Soil Conservation Service (SCS) curve number method (SCS, 1985) was used to account for the precipitation losses.The SCS curve number method is amongst the more widely used methods of assessing the effect of land use change in the hydrological response (Defries and Eshleman, 2003).

Study area and resources
Located in the central part of the island of Cyprus, the catchment area of the study is about 110 km 2 in size with an average slope value of 7.19 % (Fig. 1a).Specifically the study area is situated between longitudes 33 • 11 24.28 and 33 • 26 31.52 and latitudes 34 • 54 36.74 and 35 • 2 52.16 .In the past few years, the specific catchment area has been undergoing intensive land use change due to rapid economic growth and urbanization.The island of Cyprus is located in the northeasternmost corner of the Mediterranean Sea and, therefore, has a typical eastern Mediterranean climate: the combined temperature-rainfall regime is characterized by cool-to-mild wet winters and warm-to-hot dry summers (Michaelides et al., 2009).
For the purposes of the study, two ASTER images were utilized with 10 yr time interval in order to monitor the multitemporal urban spawl phenomenon.For this study, the first three spectral bands were used (visible and near-infrared (VNIR) and short-wavelength infrared (SWIR)) with spatial resolution of 15 m.The exact acquisition dates of the images were 12/05/2000 and 06/04/2010.
The meteorological data were provided from the Meteorological Service of Cyprus.More specifically, time-series rainfall data of six rain gauge stations for a period of 20 yr  were provided.Flow data from three stream gauges stations (Kotsiati, Nisou, Potamia) were provided from the Water Development Department of Cyprus (Table 1).From the available data, the most severe hydrological events were used to calibrate the hydrological model.The spatial distribution of rain and stream flow gauge stations in the vicinity of the catchment area is shown in Fig. 1b.

Methodology
Initially, pre-processing techniques such as geometric, radiometric and atmospheric corrections were applied to both satellite images.Sophisticated classification techniques, such as object-oriented analysis, were implemented and diachronic LULC maps of the study area were developed (for the time period of [2000][2001][2002][2003][2004][2005][2006][2007][2008][2009][2010].Using these LULC data sets, the CA-Markov algorithm was applied and the LULC map of 2020 time period was predicted.In addition, the area's soil map was developed in a GIS environment.In the following, LULC maps, soil map, DEM, meteorological and flow data for different time periods were incorporated in HEC-HMS (Hydrologic Engineering Center-Hydrologic Modeling System) hydrological software for implementing hydrological modeling in GIS environment.Firstly, the model was cali-brated for three precipitation events and then was validated for a major flood event that occurred in 2003.Finally, the same precipitation data of 2003 were once again incorporated to estimate the updated curve number map for 2010 and 2020 and assess the watershed's hydrological response under different land cover conditions.The overall proposed methodology is presented in the flowchart of Fig. 2.

Pre-processing techniques
Regarding the preprocessing of the images, geometric corrections were carried out using standard techniques with several ground control points (GCPs) and a second-order polynomial fit.For this purpose, detail topographical maps (scale 1 : 5000 and 1 : 2000) were used to track the position of GCPs in conjunction with the digital shoreline of Cyprus extracted from detail topographic maps (scale 1 : 5000).At a next step, radiometric corrections were applied to ASTER images.Radiometric corrections are essential for satellite images, since illumination's changes (e.g., Earth to Sun distance correction) and changes in viewing geometry (e.g., Sun elevation correction) should be minimized in multi-temporal analysis.Thus, the DN (digital number) values were converted to reflectance values.Atmospheric correction is considered to be one of the most difficult techniques since the distributions and intensities of these effects are often inadequately known.Despite the variety of techniques used to estimate the atmospheric effect, atmospheric correction remains a hard task in the preprocessing of image data.As it is shown by several studies (Hadjimitsis et al., 2004(Hadjimitsis et al., , 2010;;Agapiou et al., 2011), the darkest pixel (DP) atmospheric correction method can be easily and accurately applied either by using dark and nonvariant targets located in the image or by conducting in situ measurements.In the present study, water dams were used as dark targets, and the darkest pixel correction was applied to both images.

Object-oriented classification
According to Alexakis et al. (2012a), spectral mixing in satellite images between marl/chalk geological formations and urban areas was widely observed in Yialias catchment area and especially in its downwards part.This problem is clearly denoted in the spectral signature diagram derived from the use of the handheld GER 1500 spectroradiometer.The GER 1500 spectroradiometer can record electromagnetic radiation between 350 nm and 1050 nm.For the purposes of this study, different targets (Fig. 3) from the Yialias watershed basin were selected and their corresponding samples were collected (soil (Marl/chalk sample) A, B, C, -roof -tile).Laboratory spectroradiometric measurements were consecutively carried out for each different sample.A final mean measurement corresponding to ASTER bands was extracted from the 10 measurements, and for this transformation the relative spectral response (RSR) filters of ASTER satellite were used.RSR filters describe the relative sensitivity of the satellite sensor to radiance at various parts of the electromagnetic spectrum (Wu et al., 2010), and their values range from 0 to 1. Band-pass filters are used in the same way in spectroradiometers in order to transmit a certain wavelength band and block others.Therefore, the broadband reflectance from the spectroradiometer was calculated based on the wavelength of ASTER sensor and the RSR filter as follows: where R band is the reflectance at a range of wavelength (e.g., Band 1), R i the reflectance at a specific wavelength (e.g., i = 450 nm), and RSR i the relative response value at the specific wavelength.
According to the results (see Fig. 4), there is a spectral similarity between soil and urban signatures.This fact clearly depicted the unavoidable need for the application of alternative classification techniques such as object-oriented classification.
Object-based classification methodology begins with the construction of segmented objects at multiple levels of scales as major units for image analysis, instead of using a per pixel basis of single scale for image classification (Stow et al., 2007).Therefore, one of the main advantages of using objects in classification process is that, in addition to spectral information, objects have numerous geographical and geometrical features attributed to them, including shape, length and topological entities, such as adjacency (Baatz et al., 2004).A group of pixels having similar spectral and spatial properties is considered as an object in the object-based classification prototype.
Initially, the object-based approach involves the segmentation of image data into individual objects.According to Willhauck et al. (2000) and Alexakis et al. (2012b), the image segmentation is mainly influenced by the parameters of scale, color and form.The size of the image object is determined according to a scale parameter, which allows for more objects to be merged and fused as values become larger.The form parameter is a combination of the smoothness and the compactness of segment's borders.The weighting of these parameters establishes the homogeneity criterion for the object patterns (Whiteside et al., 2011).
In this study, appropriate values were assigned to three key parameters: shape, compactness and scale.The shape parameter, which adjusts spectral homogeneity compared to object's shape, was set to 0.1 in order to give less weight to shape since urban and marl/chalk classes did not have a specific shape.The compactness parameter balances compactness/smoothness and determines the object shape between compact edges and smooth boundaries.It was set to 0.5 in order to balance equally the compactness and smoothness of the objects.However, the most crucial factor of segmentation process is the adjustment of scale, which controls the object size.Thus, the higher the value of scale parameter, the larger the extracted segmented objects.Following the evaluation of several different scale parameters, a value of 10 was selected.Thus, the images were initially segmented (Fig. 4a) into object primitives or segments using the multi-resolution algorithm, which according to Baatz et al. (2003) follows the fractal net evolution algorithm.
The classification process identified and implemented seven major different classes (agriculture generic (general unidentified croplands), agriculture close grown (dense cultivated croplands (usually wheat)), herbaceous (mixture of grass, weeds, and low-growing brush, with brush the minor element), mixed forest, olive trees, urban fabric, water) by using the nearest neighbor classification algorithm (Fig. 4b,  c).The main advantage of the nearest neighbor classification algorithm is that it allows unlimited applicability of the classification process to other areas and requires only the additional selection of new training samples until a satisfactory result is obtained.
At the end, with the specific classification approach, the kappa coefficient values were increased from the initial val-ues of lower than 0.6 for both 2000 and 2010 images to 0.78 and 0.80, accordingly.

Soil map
The soil map was constructed in GIS environment according to local hydrogeological maps regime, local soil data and HEC-HMS soil classes.The final map was a three-class generalized soil map of the area (Fig. 5).Specifically, Vergennes is a very deep, moderately well drained soil of sandy loam composition concerning the specific area.Windsor consists of very deep, excessively drained soil, which for the specific area is of coarse sandy loam composition.Covington consists of very deep, poorly drained soil that is formed in calcareous glaciolacustrine and estuarine clays mainly found in the northeastern part of the basin.

Prediction of urban sprawl phenomenon
The stochastic Markov chain model was implemented to test whether urban expansion could be predicted for 2020 using the ASTER data of both 2000 and 2010.According to Ahmed and Ahmed (2012), this kind of predictive land cover change modeling is appropriate when the past trend of land cover is known.
Urban growth modeling has evolved over recent years to capture increasingly well the details of urban morphology and structure at a qualitative as well as a quantitative level (Rimal, 2005).Land use change transition probability in Markov analysis indicates the probability of making a transition from one land use class to another within two discrete time periods.The Markov chain equation was constructed using the land cover distributions at the beginning (Mt) and at the end (Mt + 1) of a discrete time period, as well as a transition matrix (MLc) representing the land cover changes that occurred during that period.In a Markov chain the probability of the next state is only dependent upon the current state.This is called Markov property as shown in the Eq. ( 2)  (Ahmed and Ahmed, 2012): where the probability Markov chain ξ 1, ξ 2, . . .can be calculated as (3) Under the assumption that the sample is representative of the region, these proportional changes become probabilities of land cover change over the entire sample area and form the transition matrices.However, the model is not spatially explicit and does not provide an explanation of the processes leading to changes and overlooks the spatial distribution of land cover in predicting land cover (Lambin et al., 1994;Adhikari et al., 2012).
The transition probability matrix explains the probability that each land cover category will change into another category.Specifically, it refers to the number of pixels that are expected to change from each land cover type to every other type over the specified number of time units (Kityuttachai et al., 2013).CA-Markov methodology underlies dynamics of the change events based on proximity concept so that the regions closer to existing areas of the same class are more probable to change to a different class (Memarian et al., 2012).A combined Markov and cellular automata (CA-Markov) was used to predict the area's land cover regime for the year 2020.The CA-Markov analysis was run to test a pair of land cover images (2000 and 2010) and output the transition probability matrix (Table 2).As it is indicated in Table 2, the mixed forest and olive tree classes have significant possibility to change to urban land cover in the near future.
After the implementation of CA-Markov model, the land use area statistics were thoroughly examined (Fig. 6a).The results indicated a steady increase of urban land cover within the catchment area, which is expected to range in a percentage of around 100 % until 2020 as well as a respective increase of agricultural generic and olive tree classes.In addition, significant decrease of agricultural close grown land cover is recorded for 2010 and is predicted for 2020.Those characteristic changes are also presented in the relative differences (%) diagram (Fig. 6b).

The hydrological model HEC-HMS
The Hydrologic Modeling System (HEC-HMS) is designed to simulate the precipitation-runoff processes of dendritic watershed systems.It is developed to be applicable in a wide range of geographic areas for solving the widest possible range of problems.This includes large river basin water supply and flood hydrology and small urban or natural watershed runoff (HEC-HMS User's Manual, 2001).
The basic rainfall runoff processes that need to be simulated in HEC-HMS for flood flow estimation using rainfall data as input, are the rainfall losses and the transformation of excess rainfall to runoff.For calculating rainfall losses, the SCS curve number method was used and for the transformation of excess rainfall to runoff the ModClark method was used both applied in GIS environment.

The ModClark method
The modified Clark (ModClark) model in HEC-HMS is a distributed parameter model in which spatial variability of characteristics and processes are considered explicitly (Kull and Feldman, 1998;Peters and Easton, 1996).This model accounts explicitly for variations in travel time to the watershed outlet from all regions of a watershed.The ModClark algorithm is a version of the Clark unit hydrograph transformation modified to accommodate spatially distributed precipitation (Clark, 1945).Runoff computations with the ModClark model explicitly account for translation and storage.Storage is accounted for within the same linear reservoir model incorporated in the Clark model.Translation is accounted for within a grid-based travel-time model t cell = t c × (d cell / d max ) (HEC, 2000), where t c is the time of concentration for the subwatershed and is a function of basin's length and slope, t cell is the travel distance from the cell to the outlet, and d max is the travel distance from the cell furthest from the outlet.The method requires an input coefficient for storage, R, where R accounts for both translation and attenuation of excess precipitation as it moves over the basin toward the outlet.Storage coefficient R is estimated as the discharge at the inflection point on the recession limb of the hydrograph divided by the slope at the inflection point.The translation hydrograph is routed using the equation: where Q(t) is the outflow from storage at time t, t is the time increment, R is the storage coefficient, I (t) is the average inflow to storage at time t and Q(t − 1) is the outflow from storage at previous time (t − 1).

SCS curve number loss method
The SCS (Soil Conservation Service) curve number loss method is a simple, widely used and efficient method for computing excess rainfall (direct runoff) from a rainfall event in a particular area.The curve number is based on the area's hydrologic soil group, land use, treatment and hydrologic condition, with the first two having the greatest importance.The SCS runoff curve number (CN) method is described in detail in National Engineering Handbook,  , 1985).The SCS runoff equation is where Q is the runoff volume, P the precipitation volume, IA is the initial abstraction and S capacity.
A linear relationship between IA and S was suggested by SCS (1985), as shown in Eq. ( 6).
where λ = initial abstraction ratio.With λ = 0.2 in Eq. ( 3), Eq. ( 2) is transformed into the following equation: For convenience in practical applications, S is mapped into a dimensionless parameter CN (i.e., the curve number), which varies in the more appealing range between 0 and 100.The chosen mapping equation is presented as follows, for SI units:

Nash-Sutcliffe efficiency E
The efficiency E proposed by Nash and Sutcliffe (1970) is defined as one minus the sum of the absolute squared differences between the predicted and observed values normalized by the variance of the observed values during the period under investigation.It is estimated by equation where O indicates observed and P predicted values; bars indicate mean values.The normalization of the variance of the observation series results in relatively higher values of E in catchments with higher dynamics and lower values of E in catchments with lower dynamics.To obtain comparable values of E in a catchment with lower dynamics, the prediction has to be better than in a basin with high dynamics.The range of E lies between 1.0 (perfect fit) and −∞.A result lower than zero indicates that the mean value of the observed time series would have been a better predictor than the model.

Phase error (PE)
Phase error is defined as the difference in hours between the peak of the observed and the simulated flow.

Peak discharge error (PD err )
Peak discharge error is defined as the percent difference between the observed and the simulated peak discharges: 6 Case study

Model setup
The HEC-HMS model was set up in distributed mode, enabling the utilization of the spatial information of the land use via the curve number coefficient.The rainfall losses component was based solely on the SCS curve number method (USDA SCS, 1972).This method assumes an initial abstraction before ponding that is related to curve number.Curve numbers in this study were determined from USDA National Engineering Handbook (USDA-SCS, 1972).The curve number method in HEC-HMS relates runoff to soil type, land use management and antecedent soil moisture conditions.The transformation method used was the modified Clark that considers the spatial variability of characteristics and processes explicitly.The curve number was estimated using a 10 m resolution digital elevation model, land use classification for 2000 and soil classification of the area.
Yialias basin was modeled using a three subbasin setup following the available flow gauge locations within the basin.The outlets of the subbasins were set at Kotsiatis (75.15 km 2 ), Nisou (21.71 km 2 ) and Potamia (16.29 Km 2 ).

Hydrological data
The HEC-HMS model was calibrated using three available rainfall-runoff events (2000,2001,2004), while was validated using a recorded flood event.A list of the rainfall runoff events is given in Table 3.Four precipitation events were selected for the calibration-validation of HEC-HMS hydrological model.The calibration events were the most intense that could be found in the recorded data.Three events were selected, to calibrate the model for the flood of 2003, which served as validation event.The specific event occurred in the watershed's urban area (downstream) between 12 and 13 (peak time) of February of 2003.During this event, a driver of a school bus was killed and much damage was caused all over the catchment area.
The hydrological characteristic of each event is presented in Table 3.The three events (dated 2000, 2001 and 2004) served for the calibration of the hydrological model.The calibrated model was then evaluated for its performance on the fourth event of 2003, which was a major flood event of the basin.The total precipitation depths (as estimated by the areal interpolation of the available rain gauge data for the entire period of the rainfall events) are also given in Table 3, along with the total duration of the event.It can be observed that the flood event of 2003 had the greatest precipitation height compared to the rest of the calibration events.To identify the driving forces of the flood event, the return period of each maximum hourly rainfall rate was estimated for each rainfall station and event (Table 3).It can be observed that the flood event distinguishes from the rest of the rainfall-runoff events mainly due to the relatively high return period that occurred simultaneously at two stations (Leukara and Analiontas), compared to the rest of the events.

Result and discussion
In this study the multi-temporal land use regime of Yialias watershed in Cyprus was thoroughly searched with the use of object-oriented classification technique and application of CA-Markov model.The specific model appears to have certain advantages as well as specific disadvantages in its application.Initially, it does not require deep insight into the mechanisms of dynamic change, but it can help to indicate areas where such insight would be valuable and hence act as both a guide and stimulant to further research.On the other hand, Markov analysis ignores the forces and processes that produced the initial land use patterns, and also it assumes that changes will continue to do so in the future by sometimes ignoring social, human and economic dynamics.However, in order to give a spatial dimension to the Markov model, we applied the CA-Markov model.Through the 2000-2010 decade's analysis, results denote an increase in agricultural generic, olive tree cultivation and herbaceous areas, putting stress onto the close growth agricultural land, which is the main decreasing land use category.The forested is shown to occupy roughly the same land portion.The same tendency seems to be for the next decade affecting the potential hydrological response of the basin.Specifically, the simultaneous increase of residential areas and the decrease of agricultural close grown cover throughout the basin is expected to enhance the potential devastating surface runoff processes.
Regarding the hydrological modeling, the calibration of the model was performed using the Nash-Sutcliffe estimator (E), with respect to the correct representation of the peak discharge and the correct timing of it.The calibration and validation results are shown in Table 4.
The calibration and validation hydrographs are presented in Fig. 7.The results of the calibration show that the distributed setup of HEC-HMS model adequately describes the timing and the peak discharge of Yialias basin.E ranged from 0.9 to 0.46 between the calibration events and the three subbasins.For the validation event, the E ranged between 0.45 and 0.62.The phase error ranged between 0 and 1 h, except for Event 2 and 3 (Table 3) simulation at Potamia in which the phase error was 2 h.Finally, the peak discharge error was kept under 17.6 %, in all subbasins and calibration events, while, for the validation event, it ranged between 0 % and 6 %.
Having calibrated and validated the HEC-HMS model, the land use map of 2010 and the projected 2020 were used to estimate the changed curve number map for 2010 and 2020, respectively (Fig. 8).
The changes between the 2000 curve number and the ones of 2010 and 2020 are also demonstrated in Fig. 9.The changes indicate that, between 2000 and 2010, 2020, the area-weighted curve number for all the land use categories except the urban areas retains a relatively constant value around 53 (from CN = 52.9 for year 2000, to 52.7 and 53.3 for 2010s and 2020s respectively).In contrast, the areal weighted CNs for all the land use categories retain a more robust increasing trend from 53.8 to 55 and then to 56.2.It is shown here that the increase in the urban land use in 2010 (from 1.85 to 5 % of area) outweighs the slight decrease in the CN in the rest of the basins' land use classes.Accordingly, the 2020 projected land use shows that the CN is projected to increase from 55.0 in 2010 to 56.2.This increase by 1.2 units is both attributed to the change to non-urban land uses and the further urbanization of the basin (from 5 to 6.5 %).
All the events that were used to calibrate and validate the hydrological model were then run under changed land use/curve number conditions.The results are shown in Table 5.The results show an increase in the peak discharge.The magnitude of the increase in peak flow is different for the four simulated events and for each subbasin in the catchment.
Results for the validation event (which consisted a flood event in 2003) indicate an increase in the runoff response under the changed land use conditions of 2010.The changes were estimated to be 10.2, 7 and 11.1 % (Table 4), for the three subbasins Kotsiatis, Nisou and Potamia, respectively, in comparison to those of 2000.The outcome indicates that the runoff dynamics of the basin are changing due to the land use transition among different categories.Next, the CA-Markov chain predicted 2020 land use was used to simulate the 2003 event under future land use conditions.The results show a noteworthy increase in the peak discharge that reached 22.4, 14.6, and 19.9 % compared to the 2000 land use runs for the above three subbasins, respectively.The simulated changes in the runoff are presented in Fig. 7.The changes in the simulated peak discharges can be explained by the overall increase in the curve number of the basin for both 2010 and 2020 simulations.Moreover, the pattern of the CN increase between the urban and non-urban land use classes can stand as positive proof that the change in 2010 peak discharge is wholly attributed to the urban area increase, while the 2020 further increase is merely attributed to urban area increase as well as to the trade-off of non-urban land uses.It has to be noted that the above rationale explains in general the mechanism of the land use change effect on the peak runoff, but it accounts neither for the spatial distribution of the land use changes nor the distribution of the precipitation.

Conclusions
This study presented an integrated methodology for searching and forecasting a catchment's area hydrologic response with the use of HEC HMS model and satellite remote sensing techniques.The preliminary results denoted the crucial role of urban sprawl phenomenon as well as the significant change of land cover regime in the increase of runoff rate within the spatial limits of a catchment area and highlighted the importance of searching land use regime with the use of satellite remote sensing imageries.It was proved that the incorporation of multi-temporal remote sensing data in hydrological models can effectively support decision making in the areas of risk and vulnerability assessment, sustainable development and general management before and after flood events.In addition, the implementation of CA-Markov provided indication of the potential impact of land use change on flood vulnerability in the near future.
The comparison of observed flow results concerning the flood event of 2003 with the simulated flow results (with the use of different land use data concerning 2000, 2010 and 2020 land use regimes) proved that, in the case of "2010" and "2020" model, the runoff rates were steadily higher due to the expanded urban area cover that increased the phenomenon of surface runoff.This tendency was verified after incorporating the land use data for the 2020 time period.Knowing from past events that the area between the Nisou and Potamia is highly prone to flooding, the already increased dynamics of the surface runoff indicate higher flooding hazard for the area.Moreover, the projected changes in land use, which is simulated to increase the peak discharge by 14.6 and 19.9 % by 2020 for the Nisou and Potamia, dictate actions have to be taken to mitigate the flood hazard.
The results of this study can be used as a road map for taking specific actions in land use management changes to achieve sustainable water resources goals in the near future.The research team will continue to study the hydrological response of the catchment area with more updated meteorological and stream data as well as satellite images of higher spatial resolution.

Fig. 1 .
Fig. 1.(a) Study Area as indicated in the RGB-321 of ASTER image, (a) location of rain and stream gauges stations, (b) a digital elevation model (DEM) of 10 m pixel created with the use of orthorectified stereo pairs of aerial photos (scale 1 : 5000) covering the study area.

Fig. 4 .
Fig. 4. (a) Image segmentation of ASTER 2010 image, (b) LULC map of the study area for 2000, and (c) LULC map of the study area for 2010.

Fig. 5 .
Fig. 5. Soil map of the study area.

Table 1 .
Characteristics of the study area's rain and flow gauge stations.

Table 2 .
Transition probability matrix for each land cover class using the Markov chain equation.

Table 3 .
List of calibration and validation events.The precipitation represents the average precipitation of the entire watershed.The return period for each rain gauge for each event is also provided.

Table 4 .
Calibration and validation results of HEC-HMS.The Nash-Sutcliffe (E), phase error (PE) and peak discharge error (PDE) are presented.Negative values of phase error (PE) indicate simulated peak before the observed event.

Table 5 .
Changes in Yialias simulated peak discharge due to land cover change in 2010 and 2020 compared to the 2000 land use.