Using SCS-CN and Earth Observation for the Comparative Assessment of the Hydrological Effect of Gradual and Abrupt Spatiotemporal Land Cover Changes

Abstract

In this study a comparative assessment of the impacts of urbanization and of forest fires as well as their combined effect on runoff response is investigated using earth observation and the Soil Conservation Service Curve Number (SCS-CN) direct runoff estimation method in a Mediterranean peri-urban watershed in Attica, Greece. The study area underwent a significant population increase and a rapid increase of urban land uses, especially from the 1980s to the early 2000s. The urbanization process in the studied watershed caused a considerable increase of direct runoff response. A key observation of this study is that the impact of forest fires is much more prominent in rural watersheds than in urbanized watersheds. However, the increments of runoff response are important during the postfire conditions in all cases. Generally, runoff increments due to urbanization seem to be higher than runoff increments due to forest fires affecting the associated hydrological risks. It should also be considered that the effect of urbanization is lasting, and therefore, the possibility of an intense storm to take place is higher than in the case of forest fires that have an abrupt but temporal impact on runoff response. It should be noted though that the combined effect of urbanization and forest fires results in even higher runoff responses. The SCS-CN method, proved to be a valuable tool in this study, allowing the determination of the direct runoff response for each soil, land cover and land management complex in a simple but efficient way. The analysis of the evolution of the urbanization process and the runoff response in the studied watershed may provide a better insight for the design and implementation of flood risk management plans.

1. Introduction

Since past decades, impact assessment of land use/cover changes (LUCC) on runoff response attracts a constantly growing attention by the scientific community, due to their immense effect on flood risk increase, soil erosion evolution, hydrological regime and aquatic ecosystems functioning, transfer of pollutants, etc. [1,2,3,4,5,6,7]. These impairments classify LUCC as a key factor in global environmental change. In most studies, the terms land use (LU) and land cover (LC) are used collectively (LULC), while in fewer occasions they are interchangeably met. However, they display fundamental differences. Land use refers to the purpose for which the land is used, emphasizing on its functional role on economic activities, without describing the actual cover of the earth’s surface [8]. On the other hand, land cover refers to the cover of the earth’s surface, i.e., vegetation (natural or agricultural), bare soil or rock, artificial structures and buildings, water, ice, etc. Information on both LU and LC of a given area is a fundamental component of planning and management actions; however, focusing either on LU or LC depends on the purpose of each study. In those related to the LUCC effect on runoff response, like the current one, focus is attracted on LC, considering moreover the land management conditions that depend on LU.

Up till now, numerous studies have underlined the marked effect of LUCC on runoff response [7,8,9,10,11,12,13,14,15,16]. Decline of forested areas due to cultivations, urbanization, deforestation, and forest deterioration result in increased direct runoff, total discharge, and peak flows [13,14,15,16], which may lead to larger flood risk [7]. Contrary, expanding forest cover may increase evapotranspiration and reduce runoff [9,13]. Extended agricultural areas can increase evapotranspiration and surface flows, leading to a significant decrease in percolation and lateral flow [10,11]. Water abstractions related to agriculture and urban uses may also reduce total discharges [11,12].

Urbanization is one of the most pervasive anthropogenic LU changes, with a range of physical and biochemical impacts on hydrological systems and processes [17]. The global urban population has proliferated from 751 million in 1950 to 4.2 billion in 2018 [18]. The proportion of urban area inhabitants is expected to increase from 55% in 2018 to 68% by 2050 [18]. Thus, understanding the impacts related to urbanization is crucial for the implementation of the “2030 Agenda for Sustainable Development,” and the successful management of urban growth. Furthermore, a significant proportion of the recent urban expansion is directed towards the peri-urban environment, i.e., the space at the cities’ surroundings that merges into the rural landscape [19].

Several studies have investigated the effects of urbanization-related LUCC on hydrological processes and more importantly on runoff generation [4,17,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34]. Land development usually entails increased impervious surfaces, considered as the main cause of increased urban surface runoff [20,27], higher peak discharges [28,29], faster runoff response [21], and alteration of the hydrological regimes and the water balance in general [32,33,34]. Specifically, Hu et al. [24] reported that changes in surface runoff are most strongly correlated with changes on impervious surfaces when compared with the changes of runoff correlated with other types of land use. Furthermore, Miller et al. [21] observed that increase in impervious cover could have a much greater impact on peak flows for rural catchments than for an existing urban area. However, as it is explained by McGrane [23], the heterogeneous urban landscape at contrasting scales is complicating our capacity to quantify infiltration and evapotranspiration dynamics. Additionally, while impervious surfaces exacerbate runoff processes, runoff from pervious areas remains uncertain owing to variable infiltration dynamics.

In addition to the gradual LUCC (e.g., urbanization) effect on the hydrological cycle, several studies have also highlighted the impact of forest fires on the latter, including reduced infiltration and evapotranspiration rates, and increased overland flow [13,35,36,37,38,39]. Forest fires cause abrupt and widespread LUCC, that may have a critical effect on a watershed’s hydrological response and associated hydrological risks [35,36,38,40,41]. However, these impacts are (mostly) temporary, and the hydrological processes usually recover in a short period (1 or 2 years) after the event. Under certain conditions, such as harsh meteorological and hydrological regimes, steep relief, highly erodible soils, grazing etc., longer (5–10 years) recovery periods may be required [36,42,43,44].

One of the most widely used methods for the assessment of LUCC impact on hydrological response is the Soil Conservation Service Curve Number (SCS-CN) method, also known as the US Natural Resources Conservation Service Curve Number (NRCS-CN) method [8,17,20,24,45,46,47]. SCS-CN is one of the simplest and well-documented conceptual methods for predicting runoff and is used in many hydrologic applications, such as flood design and water balance calculation models. It features easy to obtain and well-documented environmental inputs, and it accounts for many of the factors affecting runoff generation such as soil, land use, and land management, incorporating them in a single CN parameter [45,46,48,49,50,51,52,53].

Up till now, many research efforts have utilized SCS-CN to quantify the effect of LUCC dynamics on runoff response [8,20,47,54,55]. SCS-CN is routinely used in studies investigating the impact of urbanization on runoff response [17,24,56,57]. Furthermore, the SCS-CN method is generally used in studies estimating the effect of forest fires on the hydrological response of the vegetated areas and the associated hydrological risks [36,40,41,58,59,60,61,62,63,64].

Earth Observation (EO) satellites provide high spatial resolution, systematic revisit times and synoptic view of the Earth’s surface; hence they are well suited to support several kinds of studies of local, regional, and national interest [65,66,67]. The multitemporal, high-resolution airborne and satellite data are an accurate and rapid tool to map the spatial distribution of land use/cover and provide a remarkable diachronic historical record for the examination of LULC temporal changes [68,69,70]. Therefore, they can be utilized as ideal inputs to LULC analyses and hydrological studies [69,71].

Remote sensing techniques can deliver spatial and temporal data of the land surface, appropriate for classification and analysis of watershed areas and runoff models. On the one hand, they can quantify surface factors (such as LULC type, surface roughness, etc.) and their temporal changes, and on the other hand, they can identify hydrologically critical areal phenomena for spatial model output confirmation [67,72]. LUCC in a watershed area with urbanization and deforestation phenomena (mainly due to wildfires), affects the soil’s surface characteristics (such as impervious cover and vegetation), that can influence the nature, extent, and speed of surface runoff [70,73]. Various studies have been conducted worldwide utilizing EO information with different methods and datasets, in order to classify and map LULC, urban expansion and wildfire events, and their subsequent impact on the watershed’s hydrological characteristics (runoff change, flood frequency, and severity) [72,73,74,75,76,77,78,79].

Starting from the early 1970s, the National Aeronautics and Space Administration (NASA) and the Department of Interior United States Geological Survey (USGS) mission of Landsat satellites, offered the longest continuous space-based record of Earth’s surface [80]. Moreover, the European Space Agency (ESA) Copernicus Programme, with the optical satellite Sentinel-2 (S2), made available a new land monitoring and classification mission with 13 spectral bands, having 10 m spatial analysis. The high-resolution multispectral data of Landsat and Sentinel along with very high- resolution aerial imagery, and the CORINE Land Cover (CLC) dataset, constitute a powerful toolbox for LUCC mapping, urban expansion record, wildfire severity, and postfire effects monitoring [67,73,81,82,83,84,85,86].

To the best of our knowledge, up to now, very few (if any) studies have investigated the combined effect of gradual (e.g., urbanization) and abrupt (e.g., forest fires) LUCC. The purpose of this study is to make a comparative assessment of the impacts of urbanization, causing gradual but lasting LC changes, and of forest fires, causing abrupt but temporary LC changes, as well as their combined effect, on runoff response. To that end, earth observation (to describe the spatiotemporal land cover changes between 1945 and 2018) and the well-known SCS-CN direct runoff estimation method were used. The study is performed at a Mediterranean peri-urban watershed located in the east side of Penteli Mountain, Attica, Greece. The watershed was subjected to several natural and anthropogenic disturbances such as extensive urbanization and forest fires. Furthermore, the existence of an experimental watershed (namely, Lykorema) in the upper part of the study area, which was also affected by a major forest fire, provided the opportunity to use specifically calibrated CN values based on real data before and after the fire event effect [36].

2. Materials and Methods

2.1. Study Area

The Rafina watershed is located at the eastern coast of Attica Prefecture, Greece (coordinates: UL 23°45′ E–38°05′ N, LR 24°00′ E–37°56′ N) (Figure 1) approximately 25 km northeast of the Athens city center. It is a peri-urban area in the greater southeast Mesogeia region [87], occupying a surface of 135.32 km². The Hymettus and Penteli mountains are the basin’s dominant topographic elements, delimited by shallow catchments of northwest to southeast and with northeast to southwest orientation [88]. Its northern parts are characterized by strong relief (i.e., steep slopes), while the eastern parts and the lowlands towards its outlet to Evoikos Gulf by milder topography. Mean elevation is equal to 228 m, ranging from 0 to 940 m a.s.l. The drainage network forms a quite dense dendritic pattern, developing more complex branches at the eastern and northern regions of the basin. In accordance to the climatic conditions of the broader eastern Attica, local climate is of Mediterranean type (Csa), featuring mild, rainy winters and hot, dry summers with frequent thunderstorms [89]. The bedrock comprises four major lithological units, namely, Neogene formations (Upper Pliocene–Lower Pleistocene), alluvial deposits (Quaternary; Pleistocene), marbles, and limestones [90,91]. Regarding its soils, the upper part of the study area is mostly covered by coarse-textured and very permeable soils, while the lower part is mostly covered by medium- and moderately fine-textured soils and only small parts are covered by fine-textured soils with bad drainage.

The mean annual precipitation of the study area is approximately equal to 552 mm (2003–2018) ranging from 453 to 665 mm, with the highest depth being met at the northern part of the basin, where the experimental watershed is located, and the lowest at its southern part near the Spata settlement [91]. Snowfall events in the studied area are rare. The dry period extends from May to October. Mean annual temperature fluctuates between 11 and 27 °C.

At the Northern part of the study area; the east side of Penteli Mountain in particular, an experimental watershed is established, operated by the Agricultural University of Athens (AUA), Greece. The latter is equipped with a dense hydrometeorological network, including one hydrometric station at the outlet of the watershed, and five weather stations (Lykorema stream experimental watershed, coordinates: UL 23°53′33″ E–38°04′13″ N, LR 23°56′00″ E–38°02′28″ N) (Figure 1). The watershed is operational from September 2004 up to this day, providing the opportunity to determine the CN values for various land cover types of the study area (especially for natural vegetation and impermeable surface) for the prevalent prefire and postfire conditions, based on real data [35,36,92]. The experimental watershed is characterized by increased permeability with a runoff coefficient for the natural conditions (without the effect of forest fires) of about 0.05 [49]. The Base Flow Index (BFI) for the natural conditions was estimated equal to 0.79. In the period after the fire, the estimated BFI was decreased to 0.54 [35]. The maximum total direct runoff depth observed for the postfire conditions was (54 mm), which is much higher than the maximum total direct runoff depth observed for the natural conditions (7 mm). Finally, the maximum flow rate observed during the postfire period was 81.25 m³ sec⁻¹, which is considerably higher than the maximum flow rate observed before the fire (6.9 m³ sec⁻¹) [36].

Details about the characteristics of the experimental watershed can be found in [35,36,49,92,93].

2.2. Background

2.2.1. Urbanization in the Region

Rafina catchment comprises six municipalities of the eastern Attica region, based on the Kallikratis scheme (established in 2011), namely: (a) Rafina—Pikermi (almost entirely, 83.3%), (b) Pallini (almost entirely, 98.4%), (c) Penteli (partly, 56.8%), (d) Spata—Artemida (partly, 31.6%), (e) Paiania (partly, 26.6%), and (f) Agia Paraskevi (partly, 18.3%) [94]. According to the demographic data for the region (1940, 1951, 1981, 1991, 2001, and 2011), a significant increase in the municipalities’ population has been recorded at the last 70 years, which in some cases reached values greater than 1000% (Figure 2) [95].

Furthermore, public construction projects of noteworthy scale have been made in the region, especially throughout the last two decades. The most significant of them, involve the development of the new national airport of Athens in Spata, along with an extensive supportive road network, including the Attiki Odos motorway, which connects the Attica Prefecture’s center with the southeastern part of Attica, as well as the Rafina port. In addition, various sports facilities were manufactured for the 2004 Olympic Games of Athens, and several accompanying construction projects were established in the area (like shopping malls and leisure facilities) [87].

2.2.2. Historical Fires in the Region

The basin’s northern part is forested with flammable material [59], making it exceptionally prone to fire outbreaks. Over the last decades, a notable number of wildfires has occurred in the broader study area. The fire of July 21, 1995 started at Penteli Mountain, Attica. The event caused the incineration of approximately 251 km² of pine forest, and severe damages on almost 105 buildings, yet without loss of human lives [96]. On August 2, 1998, a fire broke out at the northeast suburbs of Athens, reaching the limits of the previously (1995) burned area, and entering “deeper” into the adjacent settlements. Recorded impairments include the reincineration of already regenerated vegetation, forest burn down at the Penteli village outskirts, damages on numerous private and public buildings, and the loss of one human life [97]. The two front fires starting at Rafina, Attica on July 28, 2005 ended up burning 4 km² of land and damaging several buildings [97]. The wildfire of August 21, 2009 is considered one of the greatest events, the Attica Prefecture has ever experienced. The fire broke out in the eastern part of the Prefecture, expanding shortly over the entire northeastern Attica, heavily affecting local communities. It burned approximately 850 km² of land (pine forest predominantly), damaging 72 houses moreover. The event affected a large part of the experimental watershed of the AUA, providing the opportunity to study the effect of such events on the hydrological response of a basin and to estimate the derivative variation of CN values. Finally, on July 23, 2018, one of the most devastating fire events in the history of modern Greece took place. The fire broke out at the forested outskirts of Penteli mountain and descended its way down towards the settlement of Mati, Attica (eastern Attica coast), having affected numerous communities in its path. The extreme weather conditions and specific characteristics of the area (i.e., morphology, microclimate, town configuration, and fuel material) enhanced its magnitude, rapid spread, and overall catastrophic action. The fire caused 101 human casualties, and severe consequences on the local environment (incineration of hundreds of square kilometer area covered by pine forest), personal assets (over 1500 buildings were damaged, and several vehicles were destroyed), infrastructure networks (water supply, telecommunication, and electricity), and regional economic development [91].

2.2.3. Historic Flood Events

The study area drainage network (paired by the public hydraulic works) is considered inadequate to convey flash floods [98]. This fact, in conjunction with the rapid (ongoing) urbanization, is responsible for flood outbreaks, even during precipitation events of less intensity [87]. In this context, four significant flood events of the recent past are presented. The first one occurred on November 2, 1977, classified as the second more important event in the Attica District history, causing severe damages throughout (37 casualties, 1924 houses damaged, 172 cars carried away, and 15% of the road network destroyed) [99]. Although the phenomenon displayed intense manifestation at the study area, no impairments were registered [87]. The event on February 25, 1988 affected the local settlements of Paiania and Marathon, flooding approximately 20 km² of land. The flood of March 26, 1998 resulted in the destruction of two bridges [100]. Finally, on February 22, 2013 a flood outbroke caused by a 9-h precipitation event with recorder height of 95.6 and 13.6 mm at the Penteli and Pikermi settlements, respectively [101]. Despite the low depth of the Pikermi event, the induced flood caused several damages (on bridges, houses, and road network) at the regions of Spata and Rafina [87].

The historic floods database of the Greek Ministry of Environment, Energy and Climate Change (MEECC) contains records of additional worth-mentioning flood events in the region [102]. These occurred at the settlements of Pallini (November 8, 1991) and Nea Makri (January 24, 2003), having affected areas very close to the basin’s northeastern (November 8,1991, Paiania settlement) and southwestern (January 24, 2003, Kokkino Limanaki settlement) boarders. Unfortunately, no additional information is available regarding the magnitude and derivative effects of the events.

2.3. Materials

2.3.1. LULC Earth and Fire Mapping Observation Data

The CLC records and their LULC classification for the years 1990, 2000, 2012, and 2018 were utilized as a base to create the respective maps of the study area [103]. For the delineation of the multi-temporal LULC maps (from 1945 to 2018), several orthophotos datasets were used from 1945 (at the scale of 1:20.000) to 1994, 2007, and 2015 (at the scale of 1:1.000). Landsat-5 (L5, 30 m resolution) and Sentinel-2 (10 m resolution) satellite images of 1986, 1988, 1993, 2000, 2004, and 2018 were also utilized (

Table 1. Earth Observation datasets.

EO System	Instrument	Acquisition Date	Usage
Orthophotos Aerial Imagery	Camera	1945	LULC mapping Scale 1:20.000 (1945) and 1:1.000 (1994, 2007, and 2015)
		1994
		2007
		2015
Landsat 5	Thematic Mapper (TM)	May 22, 1986	LULC mapping (pixel size 30 m)
		August 15, 1988
		June 10, 1993
		March 25, 2000
		September 17, 2000
		February 01, 2004
		July 26, 2004
Landsat 5	Thematic Mapper (TM)	September 13, 1995	Historical Burn scars delineation (pixel size 30 m)
		September 05, 1998
		September 08, 2005
		September 03, 2009
Sentinel-2	Multispectral Instrument (MSI)	July 05, 2018	Burn scar of 2018 wildfire delineation and LULC mapping (pixel size 10 m)
		August 04, 2018

).

Moreover, four Landsat 5 Thematic Mapper images were acquired, to record the most severe historical fires that occurred in the area (Table 1). In addition, two S2 images (prefire and postfire) were acquired for the 2018 fire event and the burn scar delineation. The free and open-access Landsat and Sentinel data were obtained via the EarthExplorer system of the United States Geological Survey (https://earthexplorer.usgs.gov/) and through the European Space Agency open access hub (https://scihub.copernicus.eu/), respectively. The L5 available data were in processing level 1T, having undergone geometrical and standard terrain correction, while the S2 products were in processing level 2A, having been geometrically and atmospherically corrected (Table 1). All satellite images in this study were selected to represent clear atmospheric conditions, without cloud cover. The digital image processing and spatial analysis were accomplished utilizing ENVI (5.5, Harris Geospatial Solutions, USA), ArcGIS (10.5.1, Environmental Systems Research Institute, Redlands, CA, USA), and SNAP (6.0, European Space Agency) software.

2.3.2. Soil Data

Due to the lack of a detailed soil map for the entire study area, data from various sources were integrated and interpreted in order to create the required Hydrologic Soil Groups (HSG) map (Figure 3) according to the method documented [104]. The basic information was the soil map of Greece (scale 1:30,000) [105], which, however, covers mostly the cultivated areas. For areas that are not included in the map, various other soil data sources such as the European Soil Database (European Soil Database v2.0, scale 1:1,000,000) [106], topsoil physical properties for Europe based on LUCAS topsoil data [106,107], and information on 14 soil profiles provided by the Greek National Agricultural Research Foundation (NAGREF) were utilized. Regarding the experimental watershed, the detailed soil survey data obtained during the establishment and monitoring of the experimental watershed were utilized [49].

2.3.3. Hydrometeorological Data

Meteorological data were obtained from 13 rain gauge stations located inside or near the studied watershed (Figure 1). Five of them belong to the hydrometeorological network of the experimental watershed, operational since September 2004 [35,36,49]. Eight additional rain gauges, scattered throughout the basin, belong to the METEONET network [101], operational since 2005. In all cases, the data are recorded in a time step of 10 min. Historical precipitation data in a daily time step were also obtained from the Hellenic National Meteorological Service [108].

The hydrological data utilized for the estimation of CN values were obtained from the hydrometric station of the experimental watershed. More information about the hydrological data and the corresponding methods are provided in [35,36,109].

2.3.4. Topographic Data

The Hellenic Military Geographical Service (HMGS) provided the required topographic maps (“Kifissia,” “Rafina,” and “Athinai-Koropi” sheets) in a 1:50,000 scale. The latter were used to extract background information data layers (e.g., geomorphology, rivers, roads, infrastructures, etc.). A detailed DEM, with a spatial resolution of 5 m, 0.5 m geolocation accuracy that was developed based on the 1:5000 scale, and 4 m contour interval topographical maps of HMGS, was also used in this study.

2.3.5. Historical Data of Fires

The diachronic inventory database of forest fires since 1984 (BSM, National Observatory of Athens, http://ocean.space.noa.gr/diachronic_bsm/) and the European Forest Fire Information System (EFFIS, http://effis.jrc.ec.europa.eu/) were used for the cross-checking of significant fire events of the last 40 years in the area.

2.4. Methodology

2.4.1. Wildfires Burn Scar Mapping

The delineation of the wildfire burn scars was divided into two parts: (i) the historical fires that occurred in the area, for which L5 images were utilized and (ii) the recent wildfire of 2018 in Mati [91], for which S2 images were utilized.

Concerning the historical events, L5 images of 1995, 1998, 2005, and 2009 false color composite (FCC) images and the Normalized Difference Vegetation Index were created using the appropriate short-wave (SWIR, band 7), near infrared (NIR, band 5) and red (band 4) spectral bands [FCC: 7,4,2 as RGB and NDVI = (B4 − B3)/(B4 + B3)]. Thereinafter, the burn scars were digitized using visual interpretation procedure.

Regarding the prefire and postfire images of S2, the Normalized Burn Ratio (NBR) index was used to mark burned areas. NBR reveals fire characteristics and evaluates fire impact on vegetation. Healthy vegetation displays very high reflectance in the NIR band, and low reflectance in the shortwave infrared (SWIR) one, and on the other hand, burned areas display the opposite characteristics. High NBR values indicate healthy vegetation in general, whereas low NBR values represent bare ground and recently burned areas [91,110,111]. The NBR formula utilizes the NIR (band 8) and SWIR (band 12) spectral bands (Equation (1)).

$$\mathrm{N}\mathrm{B}\mathrm{R}=\frac{\mathrm{N}\mathrm{I}{\mathrm{R}}_{\mathrm{B}8}-\mathrm{S}\mathrm{W}\mathrm{I}{\mathrm{R}}_{\mathrm{B}12}}{\mathrm{N}\mathrm{I}{\mathrm{R}}_{\mathrm{B}8}+\mathrm{S}\mathrm{W}\mathrm{I}{\mathrm{R}}_{\mathrm{B}12}}$$

(1)

Fire frequency is the major controlling factor of vegetation existence and recovery. The information on the historical fire frequency of an area is essential in LUCC and also to explore its impact on the hydrological processes. In sight of the above, a fire frequency map was created by overlaying the polygons of calculated historical burn scars.

2.4.2. Remote Sensing LUCC

LULC was mapped at five different years, based on the CLC records, aerial orthophotos, and satellite images. These data have different spectral and spatial characteristics; therefore, a uniform classification and mapping scale was set prior to the analysis, using the CLC maps. All were resampled to the Greek Grid projection system. CLC is a well-classified continental-scale land cover map, making use of a Minimum Mapping Unit (MMU) of 25 ha. Taking into consideration its limitations (25 ha MMU, higher resolution parcels may be missed, particularly, in urban surrounding regions), a critical “update” was made in entire classes, by additionally introducing a new one named “Urban Forest—UF,” in order to distinguish such structures (mainly present at the eastern part of the basin). For that purpose, the very high-resolution aerial orthophotos and satellite images were utilized.

More specifically, the UF class characterizes the Wildland-Urban Interface (WUI) areas, i.e., “the transition zones between cities and wildland, where configurations and other anthropogenic constructions are mixed without clear boundaries with natural wildland or vegetative fuels” (FAO, 2002). According to several researchers, WUI areas constitute a hazardous and associated high-risk ecosystem, in which severe wildfires can provoke vastly disastrous impacts, especially in the Mediterranean region where special attention should be given [91,112,113].

The LULC map of 1945 was created based on the aerial orthophotos of the same year, while the maps of 1990, 2000, 2012, and 2018 by utilizing the corresponding CLC maps, which were updated using mainly aerial orthophotos and in a few cases, supplemental satellite images. The digitization was made by implementing visual photointerpretation technique, with the simultaneous use of additional resources (thematic layers, population censuses, reports, etc.) and information from the burn scars (as described in Section 2.4.1.). As a result, a database of five detailed LULC maps were developed, covering a time period of 73 years (Figure 4).

2.4.3. Runoff Response Estimation

In this study, the Soil Conservation Service Curve Number (SCS-CN) method, which is one of the most widely used and well-documented conceptual methods for the estimation of runoff response of storm events, was used [8,17,20,24,45,46,47]. This method is a standardized approach in studies investigating the effect of LUCC on runoff response, as well as flood risk assessment and flood design studies [8,17,20,24,36,40,41,47,54,55,56,57,58,59,60,61,62,63,64], because it directly accounts for many of the factors affecting runoff generation such as soil, land use, and land management, incorporating them in a single CN parameter [45,46,48,49,50,51,52,53].

In the SCS-CN method, the direct runoff depth is calculated based on the following equation:

$$\begin{array}{c}Q=\frac{{\left(P-\lambda S\right)}^2}{P+\left(1-\lambda \right)S}\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.33em}for\hspace{1em}P>\lambda S\\ Q=0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}for\hspace{1em}P\le \lambda S\end{array}$$

(2)

where, P is the total rainfall depth, λ is the initial abstraction ratio, Q is the estimated direct runoff depth, and S is the potential maximum retention depth. When S is expressed in millimeter, it is calculated from the dimensionless curve number (CN) parameter using the following relationship:

$$S=\frac{25400}{CN}-254$$

(3)

The CN values corresponding to any given soil, land cover, and land management conditions are typically selected from the various tables included in the method’s documentation, and/or the international literature [104]. However, it is always advantageous to estimate CN values from available rainfall-runoff data from the studied watershed or nearby ones of similar characteristics.

According to the method’s documentation, the determination of the CN values assumes that the initial abstraction rate is a constant value, λ = 0.2, for CN to be the only unknown parameter of the method. Some studies, having analyzed rainfall-runoff data from several watersheds, suggested that the use of a lower λ value, equal to 0.05, may lead to a better performance [114,115,116]. In this study, the initial abstraction ratio is set to the typical value (λ = 0.2) in order to maintain the compatibility with the method documentation.

In order to account for the effect of the soil’s antecedent moisture conditions (AMC) on runoff response, the estimated CN values can be adjusted according to the 5-day rainfall preceding the studied event, using the equation described in the method’s documentation [106]. The CN values corresponding to the typical conditions (AMC-II) are termed CNII, while the ones corresponding to dry (AMC-I) and wet (AMC-III) conditions are termed CNI and CNIII, respectively.

Even if AMC was initially considered as the primary source of intrastorm variability, this assumption is of questionable origin, and its use is not directly recommended anymore [117,118]. In the latest version of the documentation, the description of AMC was revised as follows: “Variability is incorporated by considering the CN as a random variable and the AMC-I and AMC-III categories as bounds of the distribution.” Other researchers also observed that in many cases, the adaptation of CN value according to AMC category does not improve the runoff prediction [46,48,49,119]. Accordingly, the CNII values were used for the current analysis.

More details on the SCS-CN method can be found in [36,49,53,92,104,109].

As it is explained in detail by Soulis and Valiantzas [109], Grove et al. [120], and Lantz and Hawkins [121], using the SCS-CN method in lumped (i.e., single composite area-weighted CN values for the watershed) and not spatially distributed (i.e., area-weighted runoff estimates from spatially distributed CNs) form can result in significant errors in runoff estimates. The main reason is the nonlinearity of the SCS-CN equation. In this study, the estimation of runoff for each scenario is made in an entirely spatially distributed form using the ArcGIS software package (ESRI, ArcGIS 10.5.1).

To this end, the areal rainfall depth for each examined event was estimated based on the data obtained from the 13 rain gauges (Figure 1), using the elevation regression method. This method was selected considering its simplicity and efficiency, the small size of the studied watershed, and the notably uneven density of the rain gauges in the study area. More specifically, many gauges reside at the upper part of the watershed, where the experimental watershed is located, but only a few at the lower part. Subsequently, the soil/LC complexes’ spatial distribution was estimated by overlaying the HSG map (Figure 3) to the respective LC, one for each scenario. The CN values that correspond to each soil–LC complex were estimated using the lookup tables included in the SCS-CN documentation [104] as well as CN values obtained from recorded rainfall-runoff data in the experimental watershed, which is part of the study area.

In most relevant studies, CN values for the various LUCC scenarios are estimated by selecting from the lookup tables, the appropriate LC and management types that correspond to each soil/LC complex [8,17,20,24,47,54,55,56,57]. However, in the case of forest fires, CN values for the postfire conditions are typically estimated either by substituting the prefire LC with other LC types that seem more representative of such conditions during the assignment of the CN values from the lookup tables [63,64] or by increasing the CN values of the affected areas [60,122,123] based on general guidelines for postfire CN determination [58,64,123,124,125,126,127].

The main problem is that the CN values corresponding to postfire conditions are not well known. Very few studies have investigated the CN variation following forest fires based on real rainfall-runoff datasets, allowing the comparison between prefire and postfire conditions [125,128], mainly due to the very limited availability of such data. Springer and Hawkins [125] concluded that “more observations and analyses following fires are needed to support definition of CNs for post-fire response and mitigation efforts.” Accordingly, this can be a significant source of error. Optimally, CN values should be estimated through the analysis of real hydrological observations comparing prefire and postfire conditions.

In this study, the results of Soulis [36] on the postfire CN variation (using detailed prefire and postfire rainfall-runoff datasets in the upper part of the study area (Lykorema experimental watershed)) are utilized in order to increase the accuracy of the yielded results. Specifically, the CN values for the soil-LC complexes existing in the experimental watershed for the prefire and postfire conditions were estimated using the “two-CN model” proposed by Soulis and Valiantzas (2012) [109] as it was generalized for heterogeneous watersheds by Soulis and Valiantzas (2013) [92]. Detailed measured rainfall-runoff data for 29 storm events producing significant direct runoff that took place from September 2004 to August 2009, corresponding to prefire conditions, and 60 events that took place from August 2009 to May 2015, corresponding to postfire conditions, were utilized. The approach followed exploits the observed correlation between calculated CN values and measured rainfall depths in order to identify the spatial distribution of CN values of watersheds taking into account their specific spatial characteristics. Specific details on the CN values determination can be found in [36,92,109].

2.4.4. Assessment Steps

In this section, the various steps followed for the comparative assessment of the impacts of urbanization and forest fires on runoff response are presented. The latter comprise:

• Burned area delineation during the study period: The area burned at each forest fire event during the last 40 years (1995, 1998, 2005, 2009, and 2018) was delineated using earth observation data.
• LUCC delineation from 1945 to 2018: To analyze the LUCC, and especially the urbanization and deforestation in the study area, the LC at five points for the years 1945, 1990, 2000, 2012, and 2018 was delineated using earth observation data.
• Design of the investigated scenarios: The scenarios investigated in this study were designed according to the LC and forest fires’ time points. The following scenarios were examined: (1) LC-1945, (2) LC-1990, (3) LC-2000, (4) LC-2012, (5) LC-2018, (6) postfire-1995, (7) postfire-1998, (8) postfire-1995+1998, (9) postfire-2005, (10) postfire-2009, (11) postfire-2018, (12) virtual fire-1945, and (13) virtual fire-2018. The postfire scenarios were based on the last available LC before each fire and the burned area of the corresponding event. Regarding the 1995 and 1998 wildfires, the case of combining the burned areas of both events was investigated, since the 1998 fire took place before the complete vegetation recovery. Finally, the two virtual scenarios examine the effect of a virtual forest fire, affecting the area burned by 50% or more of the forest fire events, based on the (a) 1945 LC and (b) 2018 LC. To this end, a fire frequency map was created considering the burned area maps of the five events that occurred during the study period.
• Selection of storm events: Towards selecting a representative set of storm events for the current analysis, all events that took place during the period the Lykorema experimental watershed was monitored were considered. In total, 89 storm events producing significant direct runoff took place from September 2004 to May 2015. The selection of these events was based on the analysis of the hydrograph at the outlet of the experimental watershed. A threshold peak flow rate was identified (0.25 m³ s⁻¹) based on visual inspection of the latter. The criterion for the separation of the storm events was set to a 3-h interval with rainfall intensity lower than 0.4 mm h⁻¹, considering the concentration time of the watershed. Subsequently, a set of 5 characteristic storm events was selected in order to simplify the analysis. To this end, the events were sorted according to their total rainfall depth, and the greater storm of each year was selected. The resulted rainfall depths were classified into five groups, and the greater event of each group was selected, in order to eliminate those with similar rainfall depths. The resulted areal average total rainfall depths were 160.7, 113.8, 103.9, 91.3, and 72.6 mm. These depths correspond to return periods ranging from 500 to 5 years according to the “Flood Risk Management Plan of the Attica Water District” [129].
• Estimation of runoff: The CN spatial distribution for each scenario is calculated, and the runoff response for each event and scenario is estimated using the SCS-CN method, as described in detail in the previous section. In this manner, the runoff response for each LULC for different rainfall depths is estimated, as well as the increase of runoff after each forest fire.
• Statistical analysis of LUCC and runoff variation: The last step includes the analyses of the LUCC over the studied period (1945–2018), the urbanization and its trends, the affected area by forest fires, and finally the effect of all these factors on runoff response. A comparison between the impact of the gradual (but lasting) LC changes caused by urbanization and of the abrupt (but temporary) LC changes caused by forest fires on runoff response was also performed. Finally, the combined effect of all LUCC is investigated and the trends in runoff response are analyzed.

3. Results

3.1. LUCC-Burned Areas in the Study Area—Deforestation

The statistical analysis of the fire events in the area (performed from 1990 until today) demonstrated the occurrence of one extreme event every 6 years. During the period of interest, the most severe of them, manifested in 1995, 1998, 2005, 2009, and 2018 (one of the deadliest and disastrous wildfires in Greek history). The fires of 1995, 1998, and 2009 affected mainly the northern part of the basin, with the burned area covering 37.1% (45.72 km²), 27.3% (33.7 km²), and 27.48% (33.88 km²) of the basin, respectively. In 1995, the 30.01 km² (24.3%) of the burned vegetation was forest and seminatural land cover, while in 1998 and 2009, 22.12 km² (19%) and 8.75 km² (7.1%) of the same type were affected, respectively. The decrease of the burned forest cover in 1998 and 2009, even if the affected areas were considerably broad, is probably due to the continuous deterioration of vegetation from forest (312, 323, and 324 codes) to sparsely vegetated areas (code 333) type. The events of 2005 and 2018 affected smaller portions of the basin, at its eastern part, covering an area of 7.8 km² (6.3%) and 5.9 km² (4.8%), respectively (Figure 5). In these two cases, the forest areas burned were small (less than 1.5 km²); contrary, the affected peri-urban (code 112) and forest urban areas (code 119), in the 2018 wildfire (

Table 2. LULC coverage (Km²) for the years 1945, 1990, 2000, 2012, and 2018.

Code	Land Use/Cover	1945	1990	2000	2012	2018
111	Continuous urban fabric	0.00	12.42	13.24	13.89	14.87
112	Discontinuous urban fabric	2.17	22.81	27.11	30.03	30.24
119	Urban forest—WUI	0.00	2.22	2.02	2.29	3.10
121	Industrial or commercial units	0.00	0.00	0.00	0.30	0.54
122	Roads and rails	0.00	0.00	0.00	1.32	1.32
133	Construction sites	0.00	2.87	0.58	0.00	0.00
142	Sport and leisure facilities	0.00	0.00	0.00	0.21	0.21
211	Nonirrigated arable land	3.46	0.00	0.00	0.00	0.00
221	Vineyards	24.23	18.07	17.35	15.69	15.16
223	Olive groves	1.01	0.61	0.61	0.61	0.61
242	Complex cultivation patterns	48.39	27.44	24.95	21.55	20.17
243	Land principally occupied by agriculture, with significant areas of natural vegetation	1.61	0.52	0.52	0.25	0.25
312	Coniferous forest	2.50	5.70	5.69	5.79	4.28
313	Mixed forest	0.00	0.00	0.00	0.00	0.22
321	Natural grassland	0.15	0.00	0.00	0.00	0.00
322	Moors and heathlands	2.39	0.00	0.00	0.22	0.00
323	Sclerophyllous vegetation	12.10	6.78	3.66	3.66	3.79
324	Transitional woodland/shrub	21.22	19.40	3.70	3.51	3.00
333	Sparsely vegetated areas	4.05	4.45	23.87	23.97	25.52
334	Burnt areas	0.00	0.00	0.00	0.00	0.00
	Total	128.23	128.23	128.23	128.23	128.23

), were tremendously large in proportion to the area burned, with dramatic results—especially in the 2018 event (Mati area).

The forest coverage of the area, regarding the coniferous (pine) forest (code 312), increased from 1945 to 1990, mainly due to the abandonment of the cultivated areas (Figure 5). After that year and until 2018, it constantly decreased (25%) because of the forest fires and the population growth. Analogously, the sclerophyllous vegetation and transitional woodland/shrub (code 324) areas were drastically decreased from 33.32 to 6.79 km² (79%), mainly converted to sparsely vegetated areas and urban fabric. A small portion of this change involved the transformation of some forest areas to forest urban area at the eastern part of the basin (Table 2 and Figure 6).

From the fire frequency map (Figure 5f), it was noted that approximately 13.68 km² (11.1%) and 26.58 km² (21.6%) of the total basin area have been burnt two and three times in the time span from 1995 to 2018, respectively, while 1.24 km² have been burnt up to five times.

3.2. LUCC-Urbanization

Concerning the two principal municipalities occupying 54% of the catchment area, for Rafina-Pikermi, a population increase of 1145% and 645% was estimated from 1940 to 2011 and 1981 to 2011 (from 1769 and from 3134 to 20,266 residents), respectively; while for Pallini, the increase for the periods 1940–2011 and 1981–2011 was 6032% and 1616% (from 902 and from 3367 to 54,415 residents), respectively. The overall population increase in Rafina basin since 1940 was estimated at approximately 1600% (from 7000 to 112,000 residents), while the growth since 1981 was around 620% (Figure 2). This vast urban expansion and the massive increase in construction activity, accompanied by forest vegetation coverage reduction, has drastically altered the landscape and hydrological characteristics of the basin.

The LULC maps analysis revealed a lot of crucial LUCC throughout the time span under consideration. The vast expansion of artificial surfaces (from 2.17 km² in 1945 to 45.11 km² in 2018), the transformation of forest areas into forest urban areas (3.1 km² in 2018), the decrease of the cultivated fields (from 78.70 km² in 1945 to 36.20 km² in 2018, 46%), and the significant increase of sparsely vegetated areas (from 4.05 km² in 1945 to 25.52 km² in 2018), as a result of the repeated forest fires that gradually degraded forested areas (Figure 6 and Figure 7a), are the most notable changes in the basin.

The comparison of Figure 2 and Figure 7b leads to the conclusion that there is a direct link between the artificial surfaces and sparsely vegetated areas increment with the population growth. As can be perceived from the fitted sigmoid curve in Figure 7b, the rapid increase of these areas coincides with the outset of rapid population growth (in the early 1980s). It should be noted that no data points are available between 1945 and 1990; however, the sigmoid curve seems to describe better the urbanization process, considering, moreover, the population growth and the information about the development of the studied area. From the analysis of Figure 7b, it is conspicuous that the rate of increase is reducing in the late 2000s; however, this reduction could be partly attributed to the financial crisis that affected Greece at the same period and up to now.

3.3. Runoff Response

The investigation of the impacts of urbanization, of forest fires, and their combined effect on runoff response was based on 13 LULC scenarios and 5 characteristic rainfall events using the SCS-CN direct runoff estimation method. To this end, the CN spatial distribution for each scenario was calculated based on the HSG map (Figure 3), on the corresponding LULC map (Figure 6) and, for the case of postfire conditions, on the corresponding burned area map (Figure 5). Characteristic examples of the resulted CN maps are presented in Figure 8. Specifically, the CN map for the beginning of the studied period that corresponds to preurbanization conditions is presented in Figure 8c and the CN map for the ending of the studied period that includes the effect of urbanization is presented in Figure 8a. CN values are generally higher in most parts of the watershed, mostly due to the extension of artificial surfaces (lower parts of the watershed) but also due to the degradation of the vegetation cover as a result of the repeated forest fires (north part of the watershed). The spatial average value of CN increased from 70 to 80.4 (

Table 3. Spatial average CN values and runoff response for all the investigated scenarios and storm events.

Fire Event	CNII	Runoff (mm)
Year	Spatial Average	Event 1 (160.7 mm)	Event 2 (113.8 mm)	Event 3 (103.9 mm)	Event 4 (91.3 mm)	Event 5 (72.6 mm)
1945	70.0	81.0	46.7	40.1	32.0	21.0
1990	74.8	93.7	57.0	49.8	40.8	28.3
1995 postfire	80.6	106.3	65.5	57.3	47.1	32.7
1998 postfire	80.1	105.0	64.3	56.2	46.1	31.9
2005 postfire	79.9	104.9	64.5	56.4	46.4	32.3
2000	79.5	103.6	63.4	55.4	45.4	31.4
2009 postfire	82.4	110.6	68.7	60.2	49.7	34.7
2012	79.9	104.7	64.4	56.3	46.3	32.2
2018	80.4	105.9	65.2	57.1	46.9	32.7
2018 postfire	80.7	106.8	66.0	57.8	47.6	33.2
1945 virtual fire	75.5	92.9	54.5	46.9	37.6	24.8
2018 virtual fire	82.6	111.3	69.4	60.9	50.3	35.4

). When it comes to the effect of forest fires on CN spatial distribution, CN increment is more prominent in the case of less urbanized conditions (Figure 8). As can be seen in Table 3, where all the results are summarized, the spatial average CN value for the same burned area (virtual fire) increases by 5.5 units for 1945 and only by 2.2 units for 2018. It should be noted that (as it is explained in detail in the methodology section) the runoff response was calculated in spatially distributed form and the spatial average CN values are only presented as an indication.

The calculated runoff response for all examined scenarios and all rainfall events are presented in Table 3. Four characteristic examples of the runoff response’s spatial distribution for the larger storm event (Event 1, 160.7 mm) are presented in Figure 9. The higher runoff depths are observed at the higher elevations due to the higher rainfall depths. However, these depths are limited in very narrow areas covered with bare rock or with sparsely vegetated areas combined with very sallow soils. The rest of the upper part of the watershed produces from negligible (Figure 9c) to low (Figure 9a) direct runoff due to the combination of natural vegetation (Figure 6) with permeable soils (Figure 3). The lower parts of the watershed produce from low (Figure 9c) to average (Figure 9a) direct runoff depths. As can be also seen in Figure 9, some parts of the watershed produce the same runoff in all cases. These parts are either agricultural areas or small urban areas that remained unchanged in the study period.

Compering the prefire and postfire runoff response in Figure 9 and Table 3 for 2018 and 1945, it can be observed that the impact of forest fires is much more prominent in 1945. As can be also seen in

Table 4. Runoff increase from the prefire conditions to the postfire conditions for all the forest fire events and for all the studied storm events.

Fire Event	Runoff Increase from the Prefire Conditions (mm (%))
Year	Event 1 (160.7 mm)	Event 2 (113.8 mm)	Event 3 (103.9 mm)	Event 4 (91.3 mm)	Event 5 (72.6 mm)
Virtual 1945	11.9 (12.8%)	7.8 (14.3%)	6.8 (14.6%)	5.6 (14.9%)	3.8 (15.5%)
1995	12.6 (11.9%)	8.4 (12.9%)	7.5 (13.1%)	6.2 (13.3%)	4.4 (13.5%)
1998	11.3 (10.8%)	7.3 (11.3%)	6.4 (11.4%)	5.2 (11.4%)	3.6 (11.2%)
1995 and 1998	13.6 (12.7%)	9.2 (13.8%)	8.1 (14%)	6.8 (14.3%)	4.9 (14.7%)
2005	1.3 (1.2%)	1.1 (1.7%)	1 (1.8%)	1 (2.1%)	0.8 (2.6%)
2009	7 (6.4%)	5.3 (7.7%)	4.8 (8%)	4.2 (8.5%)	3.3 (9.4%)
2018	2.1 (2%)	1.6 (2.4%)	1.5 (2.6%)	1.3 (2.8%)	1.1 (3.2%)
Virtual 2018	5.4 (4.9%)	4.2 (6%)	3.8 (6.3%)	3.4 (6.8%)	2.7 (7.6%)

, where the runoff increase from the prefire conditions to the postfire conditions is presented for all the forest fire events, the runoff increment was 11.9 mm for the virtual fire in 1945 but only 5.4 mm for the same event in 2018. As it is expected, the maximum runoff depths are observed in the postfire conditions in 2018 for both the actual and the virtual fire event (Table 3 and Figure 9b). Hence, the runoff response for the postfire conditions is already very high due to the urbanization of the watershed.

The same observations are also valid for all the storm events. However, as can be seen in Table 4, the percentage changes of runoff response increase as the rainfall depth decreases. For example, the runoff increments due to the 2009 wild fire are 6.4% for Event 1 and 9.4% for Event 5. Similarly, as can be seen in Table 3, the percentage runoff increments due to urbanization from 1945 to 2018 are 31% for Event 1 and 55% for Event 5. Considering also that smaller storm events are much more frequent, this fact may lead to more frequent flood events as the drainage network (paired by the public hydraulic works) in the area is considered inadequate to convey flash floods [98].

An additional important observation in Table 3 is that the runoff increments due to the urbanization of the watershed are much higher than those due to forest fires. Furthermore, the impact of forest fires is smaller in the urbanized watershed as the runoff response is already high before the fire.

The evolution of runoff response as a result of the urbanization evolved from 1945 to 2018 for two characteristic events (events 1 and 5, the larger and the smaller storm depths) are presented in Figure 10a,b. The evolution of runoff response considering the combined effect of urbanization and forest fires for the same events is presented in Figure 10c,d. As can be seen, there is a clear increasing trend in both cases. Interestingly, the slope of the linear fit is almost identical in the case only the effect of urbanization is considered and in the case of the combined effect of urbanization and forest fires. This is an additional indication that urbanization has a dominant impact due to its lasting effect. Furthermore, as can be seen in Figure 10a,b, even if the obtained determination coefficient value is very high, the linear fit does not seem to describe well the evolution of runoff response. Considering also the population growth (Figure 2), and the evolution of the artificial surfaces and sparsely vegetated areas in the studied watershed (Figure 7), the fitting of a sigmoid curve was also attempted. In fact, the determination coefficients indicate an almost perfect fit and the sigmoid curve seems to better describe the urbanization process considering also the population growth as well as information about the developments in the studied area. For example, the rapid runoff response increment during the 1980s and the 1990s is in agreement with the population growth and with the large-scale public constructions that took place in the study area during the same period. Moreover, the stabilization of runoff response since the late 2000s coincides with the financial crisis that affected Greece at the same period. However, it should be noted that more data points are required to obtain safer conclusions. Additionally, any feature projections would be very uncertain because many socioeconomic factors may influence future trends of urbanization and runoff response.

In order to obtain a clearer picture on the impact of urbanization, the relationship between the runoff response and the artificial surfaces area was analyzed (Figure 11a). A strong relationship between the two variables was obtained, however, by the visual inspection of the graph, the increase of artificial surfaces cannot describe entirely the increase of runoff response. Accordingly, the relationship between the runoff response and the sum of artificial surfaces and sparsely vegetated area was also investigated. (Figure 11b). The extension of the artificial surfaces due to the urbanization of the watershed and the increase of sparsely vegetated areas due to the urbanization and the repeated forest fires explains almost entirely the increase of runoff response.

4. Discussion and Conclusions

High and very high-resolution EO data of recent and historical satellite images and aerial photography have excellent potential for delineation of LULC and burn scars to a local scale. Highly accurate LUCC and vegetation features can be obtained from satellite and air photo images of complex landscapes, using various digital-image processing techniques, like photointerpretation and band ratioing. For small-scale catchments, with peri-urban heterogeneous cover affected from continuous wildfires, the high and very high-resolution data offer the advantage of being able to select precisely the information useful at the scale of the hydrological modelling units. In this study, detailed LULC maps, over the previous eight decades, of 1945, 1990, 2000, 2012, and 2018 were created and a thorough LUCC analysis was accomplished. Furthermore, the burned areas of 1995, 1998, 2005, 2009, and 2018 were extracted, and their relation to LUCC was estimated.

The combined use of remote sensing and GIS proved to be an effective tool for LUCC and their interrelationship with urban growth and wildfire events. The systematic and analytical examination of LUCC can be improved to determine the nature, rate, and location of the changes and their possible effects so as to reach an excellent potential for extending environmental management and land-use planning.

The studied watershed is a peri-urban area located approximately 25 km northeast of the Athens city center. This area underwent a significant population increase (Figure 2), a rapid increase of urban LU and a corresponding decrease of forested and agricultural areas (Figure 7a), especially from the 1980s to the early 2000s (Figure 7). Urbanization is considered by far as the fastest type of land use change in Europe, more than 10 times higher than any other LUCC, while peri-urban areas are growing much faster than other urban areas [19]. Urbanization is generally considered as one of the most pervasive anthropogenic LU changes [17]. In this study, it was found that urbanization resulted to a considerable increase of artificial surfaces as well as an increase of sparsely vegetated areas to the detriment of forested areas and of agricultural areas (Figure 7a).

Urban development is usually related with increased impervious surfaces and by extension, increased urban surface runoff [20,27]. The results of the current study indicate that the urbanization process in the studied watershed caused a considerable increase of direct runoff response, especially during the period of fast urban development (Figure 10a,b). It was also observed that percentage runoff increments due to urbanization are much higher for weaker storm events. Taking also into account that weaker storm events are more frequent, a devastating impact on the frequency of flood events is expected as the drainage network (paired by the public hydraulic works) in the area is considered inadequate to convey flash floods [98].

Direct runoff response increase in the studied area is strongly correlated with the extension of artificial surfaces (Figure 11a). However, the extension of sparsely vegetated areas due to the urbanization and the repeated forest fires also has an important contribution (Figure 11b). Hu et al. [24] reported that changes in surface runoff are most strongly correlated with changes in impervious surfaces. Miller et al. [21] observed that increases in impervious cover for rural catchments, as it is also the case of this study, may have a much greater impact on peak flows than for an existing urban area. Furthermore, the combination of increased impervious cover and faster runoff routing via storm drainage network of urbanized areas can lead to a much flashier response and larger peak flows [21]. However, as it is explained by McGrane [23], the heterogeneous urban landscape is complicating the evaluation of runoff response, because, while impervious surfaces exacerbate runoff processes, the variable infiltration dynamics of the remaining pervious areas are still uncertain. For this reason, future study on the refinement of the urban land use typologies used in similar studies may assist to a better representation of the hydrological processes’ alteration due to urbanization.

Apart from urbanization, forest fires are the most frequent and widespread disturbance of vegetation, particularly in Mediterranean ecosystems where the areas burned each year are extended and the fire frequency is high [130,131]. A key observation of this study is that the impact of forest fires is much more prominent in rural watersheds than in urbanized watersheds (Figure 9 and Table 3). However, the increments of runoff response are important during the postfire conditions in all cases. Many other studies also report that forest fires may have a critical effect on vegetated watersheds hydrological response and associated hydrological risks [35,36,38,40,41].

Generally, runoff increments due to urbanization seem to be higher than runoff increments due to forest fires affecting the associated hydrological risks. It should also be considered that the effect of urbanization is lasting and therefore, the possibility of an intense storm to take place is higher than in the case of forest fires that have an abrupt but temporal (generally, 2–5 years [36,42,43,44]) impact on runoff response. It should be noted though that the combined effect of urbanization and forest fires results in even higher runoff responses, which could be critical for the associated hydrological risks. An additional important observation is that the increase of sparsely vegetated areas has a profound impact as well. Accordingly, except of their direct effect, repeated forest fires may also have a more permanent effect through the gradual degradation of vegetation in forested areas and intensification of soil erosion [91]. Fire frequency seems to be positively related to urban growth, particularly in forest urban interface areas [132].

The SCS-CN method proved to be a valuable tool in this study, allowing the determination of the direct runoff response for soil, land cover, and land management complex in a simple but efficient way, through the estimation of the corresponding CN value, based on the sound documentation and the extensive relative literature [45,46,48,49,50,51,52,53,104]. Up to now, many studies have utilized SCS-CN to quantify the effect of LUCC, e.g., urbanization and forest fires [17,24,36,40,41,56,57,58,59,60,61,62,63,64], since there is a clear link between the value of its main parameter (CN) and LULC facilitating the evaluation of various scenarios.

However, as it is explained in detail by Soulis and Valiantzas [109], the method should be applied in spatially distributed form (i.e., area-weighted runoff estimates from spatially distributed CNs), otherwise, significant errors in runoff estimates can occur. On the positive side, SCS-CN method can be very easily implemented in GIS software tools in such form. An additional advantage is that in this manner, surplus information about the spatial distribution of runoff response is also provided, making it possible to evaluate the impact of specific LUCC at various parts of the watershed.

The obtained results depend on the selected CN values for the various soil/LC and management complexes as well as the CN values’ adjustments for the postfire conditions. In this study, the existence of the Lykorema experimental watershed in the upper part of the studied area provided the possibility to evaluate the CN values used for the various conditions based on real data [36,92,110]. As an example, previous studies on the estimation of the effect of forest fires on the hydrological response, increased CN by 5, 10, 15, and 20 units depending on the burn severity [60,121] based on general recommendations that are not based on real hydrological data. However, according to the results presented by Soulis [36] for the Lykorema experimental watershed, the CN values during the postfire period in permeable watersheds may increase by 29–36 units depending on the HSG.

This study focused on investigating the impacts of urbanization, forest fires, and their combined effect on runoff response using the SCS-CN direct runoff estimation method. However, the alteration of drainage pathways may also change significantly the storm runoff dynamics. It is reported by Miller et al. [21] that “routing of storm runoff via a storm drainage network can lead to a much flashier response and higher peak flows than would be attributed by increases in impervious cover or change in land use alone.” Qi and Liu [3] also reported that studies on impact of LUCC on flood peaks should also consider roughness variations. Accordingly, additional research is needed in order to investigate the effect of gradual and abrupt spatiotemporal LUCC on flood hydrographs and peak flows.

The analysis of the evolution of the urbanization process and the runoff response in the studied watershed may provide a better insight for the design and implementation of flood risk management plans. However, future projections are very uncertain because many socioeconomic factors may influence feature trends of urbanization and runoff response. Hence, it is clear that more policy attention at the regional level and the urban–rural interface is required [19]. Finally, the use of sustainable urban stormwater management solutions such as permeable pavements, green roofs, and other green infrastructure that can attenuate stormwater [133] should also be considered.