**Temporal and Spatial Analysis of Flood Occurrences in the Drainage Basin of Pinios River (Thessaly, Central Greece)**

#### **George D. Bathrellos 1,\*, Hariklia D. Skilodimou 1, Konstantinos Soukis <sup>2</sup> and Efterpi Koskeridou <sup>3</sup>**


Received: 11 August 2018; Accepted: 7 September 2018; Published: 11 September 2018

**Abstract:** Historic data and old topographic maps include information on historical floods and paleo-floods. This paper aims at identifying the flood hazard by using historic data in the drainage basin of Pinios (Peneus) River, in Thessaly, central Greece. For this purpose, a catalogue of historical flood events that occurred between 1979 and 2010 and old topographic maps of 1881 were used. Moreover, geomorphic parameters such as elevation, slope, aspect and slope curvature were taken into account. The data were combined with the Geographical Information System to analyze the temporal and spatial distribution of flood events. The results show that a total number of 146 flood events were recorded in the study area. The number of flood events reaches its maximum value in the year 1994, while October contains the most flood events. The flood occurrences increased during the period 1990–2010. The flooded area reaches its maximum value in the year 1987, and November is the month with the most records. The type of damages with the most records is for rural land use. Regarding the class of damages, no human casualties were recorded during the studied period. The annual and monthly distribution of the very high category reaches the maximum values, respectively, in the year 2005 and in June. The analysis of the spatial distribution of the floods proves that most of the occurrences are recorded in the southern part of the study area. There is a certain amount of clustering of flood events in the areas of former marshes and lakes along with the lowest and flattest parts of the study area. These areas are located in the central, southern, south-eastern and coastal part of the study area and create favorable conditions for flooding. The proposed method estimates the localization of sites prone to flood, and it may be used for flood hazard assessment mapping and for flood risk management.

**Keywords:** historic flood data; old topographic maps; GIS; temporal and spatial distribution of flood events; marshy areas and lakes; flood hazard assessment

#### **1. Introduction**

Natural hazards are physical events that can cause significant damages to the natural and human environment. Endogenic or exogenic processes such as active tectonics and climate changes are capable of changing landforms and triggering natural hazards, which in some cases control human activities [1–11].

Floods are physical phenomena active in geological time and the result of excess runoff. When rivers overtop their banks, the excess water goes to the floodplain. Floodplains represent favorable sites for man to settle because they are fertile, level, easy to excavate and near water. These features have contributed to increased development and urbanization of floodplains, thereby increasing the chances of the flood occurrences that cause disasters. Despite the construction of flood control works such as dams, levees and channeling rivers, the flood damage has increased [12–16].

Floods occur frequently in some parts of the world. While for some areas, yearly flooding is necessary to sustain crops, for other areas, flooding spells disaster. Especially, in urban areas, floods are considered among the most dangerous natural hazards due to the increasing number of events. Their consequences are not only environmental, but social and economic, as well, since they may cause damages to urban areas and agricultural lands and may even result in the loss of lives [17]. Worldwide, flood events have caused the largest amount of deaths and property damage during the last decades [18]. In Europe, the annual average flood damage in the last two decades was about €4 billion per year [19]. In the Mediterranean countries, floods tend to be greater in magnitude compared to the inner continental countries, and they frequently cause catastrophic damages [20].

Flood occurrence can be estimated in space and time through a sound basis of knowledge acquired by the scientific use of a large number of historical documents [21]. It is possible to gather useful information about historical floods that occurred during a period ranging from 1–2 centuries by using historical data [22]. Payrastre et al. [21] used discharges of historical floods that occurred during a period ranging from 1–2 centuries in four watersheds of the French Mediterranean area. They estimated the flood peak discharges and the associated low and high bound of its possible values for the main floods of each studied area. Historic data are useful for floodplain management, flood insurance rating, emergency planning and flood risk management. A careful search for data and their suitable arrangement may result in a powerful tool for flood disaster prevention and land use planning [23–25].

One essential step in any flood hazard analysis and estimation is the flood hazard recognition. This statement includes the record of possible geomorphologic evidence of flood activity, historic accounts and records [26]. Hydrologic models using time series flood data can be applied to assess flood peaks, depths, volumes and to map flood hazard areas [27–31]. However, these methods require data that are often unavailable [32]. Alternatively, researchers have used historical documents and historic flood data to estimate flood hazard [22,33]. Tropeano et al. [22] analyzed historical documents for a statistical elaboration of the frequency of landslide and flood events in Northern Italy in the last five centuries. Diakakis et al. [33] used an extensive catalogue of flooding phenomena during the last 130 years to examine flood events in Greece. According to their statistical and spatial analysis, urban areas tend to present higher flood recurrence rates than mountainous and rural ones and an increasing trend in reported flood event numbers during the last few decades. Moreover, a diachronic appraisal of the landscape evolution throughout the study of old topographic maps may include information on ungauged historical floods and paleofloods [34,35]. Skilodimou et al. [35] used old topographic maps to investigate the causes of flooding generation in the southwestern coast of Attica in Greece. They proved that the former wetlands and the lagoon of the study area have been dried up and covered by buildings, causing flood occurrences.

This paper aims to provide recognition of flood hazard by using historic data in the drainage basin of Pinios (Peneus) River, in Thessaly, central Greece. For this purpose, a catalogue of historical flood events that have caused damages was studied. These data, along with old topographic maps and geomorphic parameters such as elevation, slope, aspect and slope curvature were evaluated to record the temporal and spatial distribution of flood events of the study area. The analysis, processing and the evaluation of the old maps, geomorphic parameters and flood events were performed in a GIS environment.

#### **2. Study Area**

The drainage basin of the study area is located in Thessaly, central Greece (Figure 1a), and it covers an area of about 11,200 km2. The elevations of the study area vary from 0–2678 m a.s.l. (Figure 1b). The hydrologic basin is drained by Pinios River, which has an approximate length of 205 km. It is the third longest river in Greece and crosses a large part of the eastern part of central Greece discharging into the Aegean Sea, where it forms a delta.

**Figure 1.** (**a**) Location map of the study area; (**b**) the elevations of the study area, the drainage network, the road network and the main settlements (modified from [36]).

Its course starts at the northwestern part of the Thessaly plain, from the confluence of the Ion and Malakasiotis rivers. It is surrounded by mountainous areas, which enclose its drainage basin and form its watershed (Figure 1b). To the north are the Titaros Mt. (1837 m) and the Kamvounia Mt. (1615 m); to the northeast are the Olympos Mt. (2917 m) and the Ossa Mt. (1978 m); to the east is the Pilio Mt. (1548 m); to the south is the Orthrys Mt. (1726 m); and finally, to the west are the Pindos Mt. (2204 m) and the Koziakas Mt. (1901 m).

The major tributaries of Pinios River are the Malakasiotis, Portaikos, Pamisos and Enippeas rivers to the west and south and the Ion, Lithaios, Neochoritis and Titarisios rivers to the north, which all drain large, heterogeneous areas, through extensive hydrographic networks (Figure 1b).

The plain area crossed by Pinios River is divided by the presence of a low-lying hill area into two parts (Figure 1b): (a) a western part (Trikala-Karditsa), whose altitude varies between 80 and 200 m, and (b) an eastern part (Larisa), whose altitude varies between 45 and 100 m. The mountainous and hilly area crossed by Pinios River has a quite rugged surface with altitudes exceeding 200 m. The most important gorges are the ones of Kalamaki (internal low-lying hill area), Rodia and Tempi (between Olympos and Ossa). Through these areas, the water flow velocity is high, while erosion is

intense. All the above, in conjunction with the fertility of the plain areas, necessitated the construction of drainage and irrigation projects.

Fifty five percent (55%) of the study is mainly rural, semi-urban area. Irrigated agricultural land comprises 56% of the total cultivated area in the basin. The intense and extensive cultivation has led to a remarkable water demand increase, which is usually covered by the over-exploitation of groundwater resources [37].

The climate in the western and central side of the drainage basin of Pinios River is continental with cold winters and hot summers with a large temperature difference. The coastal area of the basin has a typical Mediterranean climate. Summers in the study area are usually very hot and dry, and in July and August, temperatures can reach up to 40 ◦C. The mean annual precipitation over the basin is about 779 mm. Rainfall is distributed unevenly in space, and it varies from about 360 mm at the eastern part of the basin to more than 1850 mm at the western mountain peaks. Generally, the rainiest months are November and December, while rainfall is rare from June–August. The mountain areas receive significant amounts of snow during the winter months. The mean annual total flow of the river is 3500 × <sup>10</sup><sup>6</sup> m3. The hydrologic regime of the river is affected by both winter rainfalls and spring snowmelt. The river flow regime at the main stem is perennial. The mean average flow near its delta ranges from more than 150 m3/s in February and March to about 10 m3/s in August and September. The water of Pinios River is used primarily for irrigation [37,38].

The Pinios River crosses alpine formations of various units and post-alpine formations. More specifically, the alpine formations at its western part belong to the Koziakas unit, at the central part to the Sub-Pelagonian and Pelagonian units and at its eastern part to the Pelagonian and Olympos-Ossa units. At the western sub-basin, a great part of its length passes through molassic formations of the Meso-Hellenic trough (Eocene-Pliocene) and Quaternary formations, while at the eastern part, it crosses mostly Neogene and Quaternary formations. The permeable rocks and sediments of the study area are divided into two major units, the karstic one consisting of carbonate formations (limestones, marbles and dolomites) and the granule materials, which include Neocene and Quaternary sediments. Semi-permeable formations are comprised of loose to semi-coherent Quaternary and Neocene deposits along with thin-bedded limestones and ophiolites. Impermeable formations consist of schists and fine-grained Quaternary and Neocene deposits [39–41].

The plain of Thessaly has suffered from flooding since ancient times, when several structures had been constructed in order to control the Pinios River 2500 years ago. Today, despite the construction of flood protection works, floods remain an important problem of the area [42].

#### **3. Materials and Methods**

#### *3.1. Data*

In the present study, the historical flood database of Greece was used [43]. It refers to floods that occurred during the period 1979–2010. This database was required by the Directive 2007/60/EC concerning floods in Europe. In Article 4 of the Directive 2007/60/EC, the construction of a historic flood data database is assigned to every state-member as a part of the preliminary assessment of flood risk, as well as finding the sensitive to flooding areas in each country, for the database will be used for the construction of flood risk maps [44]. The database contains information such as the location of flood (county and village or locality), geographical details, dates of occurrence, characteristics of the flood (e.g., flash flood, snow melt flood, etc.), the sources of flooding (e.g., pluvial, fluvial, groundwater, etc.), the causes of flooding, the mechanism of the flood, the extent of flooded area (in km2), the maximum distance of the flood (in km), the type of damages, the total cost of the damage caused by the flood event, the number of human casualties and the class of the total damage caused by the flood (very high, high, medium or low).

Four old topographic maps of 1881 with the scale 1:200,000 were used in the current work. The used sheets are: Ioannina-Metsovon-Grevena-Kozani-Servia [45], Larissa-Elasson-Katerini [46],

Arta-Trikala-Karditsa [47] and Volos-Farsalos-Lamia [48]. Figure 2 shows a mosaic of the four old topographic maps used.

**Figure 2.** Old topographic maps of 1881, sheets: (**a**) Ioannina-Metsovon-Grevena-Kozani-Servia; (**b**) Larissa-Elasson-Katerini; (**c**) Arta-Trikala-Karditsa; (**d**) Volos-Farsalos-Lamia.

Finally, thematic maps were used. These maps are the digital elevation model (DEM) of the study area, which is presented in Figure 1 and the geological map of the study area. They were acquired from previous work [36]. The DEM map was used to produce the elevation, slope, aspect and slope curvature maps. The geological map of the study area was used to identify the semi-permeable and impermeable formations of the study area.

#### *3.2. Statistical Analyses*

Several works have used historic flood data to analyze flood frequency, flood peak discharges and estimate flood hazard maps [21,27,29,30]. This kind of flood frequency analysis requires hydro-climatic data, such as flow and rainfall. In many cases, this is difficult to do, because detailed data are unavailable. On the other hand, researchers have analyzed the temporal and spatial distribution of floods based on the count of events [22,33]. The scope of the present work is to recognize flood hazard and examine the temporal and spatial distribution of flood events in the study area, rather than to create a flood hazard map. Additionally, high resolution data such as discharge or rainfall during flood events were not readily available. For these reasons, the statistical analysis depends on the location of floods, dates of occurrences, the flooded area, the type of damages and the total damage caused by the floods.

For the description and assessment of historical floods, information was sourced from various national and regional authorities, scientific reports and newspaper articles. In general, the accuracy of the data used was considered acceptable as multiple checks between different sources were carried out to ensure consistency of the data. The location of events shows some uncertainty. According to [49], when there was no reference to a particular community, but the geographical determination was different (i.e., reference to river or torrent), location was determined on the basis of the other descriptive information. Thus, in some cases, as an event location is given the center of the municipal department. The historic database includes information about the causes of flooding, the flood mechanism and characteristics along with the maximum distance of the flood. The data of these records were not systematic and accurate, and they were not taken into account in the present work.

As a first step, a temporal statistical analysis was performed to identify the number of events and flooded area distribution during this period. The statistical significance of these variables was tested by *p*-values. Additionally, to identify the seasonal distribution of flood events, these parameters were analyzed for each month.

Flood events that occurred in the study area had adverse consequences and caused damages to economic activity. The type of damage was categorized into three categories: economic (positions of infrastructures, industrial and commercial centers), rural land use (positions of farmland with significant economic value production) and property [49]. The distributions of the number of each type of damage from 1979–2010 along with their monthly distribution were examined.

The total damages caused by the floods were categorized into five classes taking into account the number of human casualties, amount of monetary compensation (state compensation for damages to agriculture and for damages to settlements) or the size of flooded area [50]. The classes of the total damage and the limits of the above-mentioned parameters are presented in Table 1.



The annual and monthly distributions of the two above-mentioned parameters were studied during the period 1979–2010.

#### *3.3. Spatial Analyses*

Furthermore, a spatial analysis of the flood occurrences was carried out. The study area was divided into seven individual drainage basins. The spatial distribution of flood events in each drainage basin was examined. A spatial database was created, and ArcGIS 10.3 software was used to process the collected data.

Old maps may describe or simply locate areas that have flooded in the past. In many cases, former wetlands and lagoons that have been dried up cause flood events. Thus, the study of old topographic maps is very helpful to examine the causes of flood genesis [22,34,35]. The comparative observation of the old topographic maps of 1881 and the recent topographic maps of the study area led to a mapping of the previous wetlands. Marshy areas and lakes were digitized from the older edition survey maps, and they compared with the spatial distribution of current flood events. Finally, bibliographic data [36] were used to examine the association of the paleo-environment of Pinios River with the spatial frequency of floods.

The geomorphological setting of an area affects flood occurrences. Lowland morphology, gentle slopes and the flat areas create favorable conditions for flooding [51]. For this reason, the elevation, slope, aspect and slope curvature of the study area were examined. The selection of these geomorphological parameters and the determination of the class numbers, as well as their boundary values was based on the literature [51–54]. The elevation map was classified into two

categories: (i) <100 m and (ii) >100 m. Similarly, the slope was divided into two categories: (i) <5◦ and (ii) >5◦. The aspect map shows the slope direction and identifies the flat areas that have values = −1◦. Thus, the aspect map was categorized into two classes: (i) = −1◦ and (ii) >−1◦. In the case of curvature, negative curvatures represent concave, zero curvature represents flat and positive curvatures represent convex, respectively. Figure 3 shows the elevation, slope, aspect and the curvature maps and their categories. The areas with elevation < 100 m, slope < 5◦, aspect = −1◦ and curvature = 0 were combined to identify the lowest and flattest areas of the study area. This was accomplished by the intersection tool of ArcGIS 10.3.

**Figure 3.** Thematic maps of the four geomorphological parameters: (**a**) elevation; (**b**) slope; (**c**) aspect; (**d**) curvature.

#### **4. Results**

#### *4.1. Temporal Distribution of Flood Events*

A total number of 146 flood events were recorded in the drainage basin of Pinios River during the period 1979–2010. Figure 4 shows the temporal distribution of flood occurrences during this period.

The number of flood events reached its maximum value (38) in 1994. The value of flood events was high (36) in the year 1987 and was relatively high (19) in 2002. In 1994, flood events were two-times higher than in 2002. In the 1980s, a total number of 40 flood events were recorded, while in the 1990s and 2000s, 55 and 47 were observed, respectively. Thus, there was a trend of increasing flood occurrences over the last two decades. The statistical analysis of the number of events showed that *p* = 0.009. Thus, the data used have values *p* < 0.05 and are statistically significant.

Figure 5 presents the seasonal distribution of flood events. The statistical analysis proved that October contained the most flood events (42). High values of flood occurrences were observed in March (35), in June (11) and in November (10). In October, the number of flood events was four-times higher than in June and November.

**Figure 4.** The annual distribution of flood events from 1979–2010.

**Figure 5.** The monthly distribution of flood events during 1979–2010.

On the other hand, substantial differences appeared in the temporal distribution of the flooded area (Figure 6). It reached its maximum value (87.4 km2) in 1987. Additionally, relatively high values were calculated in the years 1994 (69.5 km2) and 1997 (64.8 km2). The *p*-value was found to be 0.01, and thus, the data used were statistically significant. Regarding the seasonal distribution of flooded

area (Figure 7), it reached its maximum value (136.5 km2) in November. It presented relatively high values in March (83.2 km2) and October (73.1 km2).

**Figure 6.** The monthly distribution of the flooded area during 1979–2010.

**Figure 7.** The monthly distribution of the flooded area during 1979–2010.

#### *4.2. Temporal Distribution of Damages*

Concerning the type of damages, their temporal and seasonal distributions are presented, respectively, in Figures 8 and 9. It is important to note that no human casualties were observed during the studied period. The values of economic damages were low during the studied period and reached the maximum value (3) in the years 1987, 1990 and 2010 (Figure 8). Similarly, their monthly distributions had relatively low values. The maximum value (4) of economic damages was recorded in March (Figure 9).

**Figure 8.** The annual distribution of type of damages in the study area during 1979–2010.

**Figure 9.** The monthly distribution of the type of damage.

On the contrary, the values of rural land use were the highest among all the types of damages (Figure 8). It reached its maximum value (37) in 1994, while it had relatively high values in the years 1987 (33) and 2002 (17). October was the month with the most records (34) in this type of damage (Figure 9). Additionally, high values were observed in March (12) and June (11). The maximum value (10) of the property damages was recorded in 2005 (Figure 8). March included the most number of property damages (17), while relatively high values (Figure 9) were recorded in October (six).

Figures 10 and 11 show the temporal and monthly distribution of categories of damages in the study area. The maximum value (10) of the very high category was observed in 2005, whilst its relatively high values were recorded in the years 1979 (three) and 1987 (three). Moreover, this category reached its maximum value (10) in June; a relatively high value (four) was calculated in November. The high class of damages reached its maximum value (six) in 1987 and a relatively high value (three) in the year 1994. March contained most records (five) of this class, and a high value was computed in October. The medium class reached its maximum value (24) in 1994, and a high value (12) was observed in 1987. The month with the most records (28) of this class was October; while a high value (13) was calculated in March. The annual and monthly maximum values of the low damages were observed in the year 1987 (14) and in March (13). Finally, the maximum value of the unknown class was recorded in the year 2002 (17).

**Figure 10.** The annual distribution of categories of damages in the study area during 1979–2010.

**Figure 11.** The monthly distribution of categories of damages during 1979–2010.

#### *4.3. Spatial Distribution of Flood Events*

The drainage basin of Pinios River were divided into seven sub-basins, which were: (a) the upper reaches of Pinios River; (b) Trikala-Neochoritis River; (c) Karditsa-Enippeas River; (d) Larisa-Karla; (e) Tempi Valley; (f) Titarisios River; and (g) the lower reaches of Pinios River. The spatial distribution of flood events in each drainage basin was examined. Table 2 shows the spatial distribution of flood events expressed as the number of events per drainage basin during the studied period. The most of floods were located in the drainage basin of Karditsa-Enippeas River.

**Table 2.** Spatial distribution of flood occurrences in each drainage sub-basin of the study area during the period 1979–2010.


The mapping of the marshy area and lakes of the study area was based on the topographic map of 1881. They covered an area of about 161 km2. Moreover, according to [36], two paleo-lakes existed in the study area during the Quaternary. The former marshy areas and lakes along with the boundaries of two paleo-lakes were compared with the historic flood events (Figure 12).

**Figure 12.** Marshy areas, lakes, paleo-lakes and lowest and flattest areas vs. flood occurrences in each drainage sub-basin of the study area during the period 1979–2010.

In the drainage basins of Trikala-Neochoritis River and Larisa-Karla, many flood events (12 events) were located in marshy areas, Lake Karla and paleo-lakes. In the drainage basin of Karditsa-Enippeas River, many floods (eight events) were located in or near marshy areas and Lake Xinias. In the drainage basins of Titarisios River and the lower reaches of Pinios River, few floods (two events) were located in or near marshy areas and Lake Askouris. Finally, in the drainage basins of the upper reaches of Pinios River and Tempi, no floods were associated with marshy areas.

The intersection of areas with an elevation <100 m, slope <5◦, aspect = −1◦and curvature = 0 is shown in Figure 12. The lowest and flattest areas were mainly located in the drainage basins of Trikala-Neochoritis River, Larisa-Karla and the lower reaches of Pinios River. In the drainage basins of Trikala-Neochoritis River and Larisa-Karla, they were mainly located in the areas where the former marshy areas and lakes were located. Regarding the spatial distribution of flood events, many flood events (61 events) were located in these areas.

#### **5. Discussion and Conclusions**

In the present study, the historical flood database of Greece and the old topographic maps of 1881 were used. These data were combined in a GIS environment to analyze the temporal and spatial distribution of flood events. The study area was the drainage basin of Pinios River, in Thessaly, central Greece.

The database contains floods that occurred between 1979 and 2010. A total number of 146 flood events were recorded in the drainage basin of Pinios River during this period.

The number of flood events reaches its maximum value (38) in the year 1994 (Figure 4). This may be related to flood events that occurred on 21–22 October 1994. No information about the causes of flooding, the flood mechanism and characteristics along with the maximum distance of the flood were included in the historic database. According to [55], this flood was caused by a severe storm. This intense storm caused flooding along the Pinios River and its tributaries, as well as inundation of agricultural and residential areas.

The statistical analysis showed that the data used are statistically significant. The number of flood occurrences shows a rising trend during the period 1990–2010. The study of the monthly distribution of the events demonstrates that October contains the most flood events, while March is the second month with many flood records (Figure 5).

Regarding the flooded area, it reaches its maximum value (87.4 km2) in the year 1987 (Figure 6). This is attributed to flood event that occurred on 24–27 March 1987. The flood of March 1987 was caused by intense rainfall that produced direct runoff, as well as snowmelt. The water level of the Pinios River measured during the flood event exceeded 6.3 m in Kalamaki Valley and 8.0 m at Tempi Valley when the normal water level was less than 1 m in both locations (Figure 1b). Owing to the narrow riverbed at both locations (Kalamaki and Tempi), significant parts of the Thessaly plain upstream of each location were inundated [56].

Although October is the most flood-prone month, substantial differences appear in the monthly distribution of the flooded area. It reaches its maximum value (136.5 km2) in November (Figure 7). This difference may be related to several reasons such as the difference in storm characteristics.

In the study area, the economic damages have low values. The type of damage with the most records is the rural land use. The maximum value is recorded in the year 1994, while the month with the maximum value is October (Figures 8 and 9). According to [55], during the flood of 21–22 October 1994, more than 70 houses in about 20 communities were totally destroyed by the flood, more than 200 suffered severe damage and another 90 minor damage, whereas 80 km2 of agricultural land (cotton fields) were flooded. Concerning the property damages, the maximum value is observed in the year 2005, while March contains the most events (Figures 8 and 9). A flood event occurred in the study area on 16 June 2005, which caused serious damage to infrastructures and properties.

Regarding the class of damage, no human casualties were recorded in the study area between 1979 and 2010. According to [57], the most destructive events regarding human victims occurred in the city of Trikala in 1907. This flood caused at least 200 fatalities. Thus, the classification of the total damages is based on the amount of monetary compensation and flooded area of the historical flood events. The annual maximum value of the very high class is observed in the year 2005, while its monthly maximum value is recorded in June (Figures 10 and 11). This fact is related to the flood event of 16 June 2005. The amount of monetary compensation was 1,027,146 euros and was given for damages to infrastructures and properties [43]. The maximum value of the high class is recorded in the year 1987, and the maximum value of its monthly distribution is observed in March. This may be related to the flood event of 24–27 March 1987. Finally, the medium class has its maximum values in the year 1994 and in October, and this is related to the flood event of 21–22 October 1994.

The study area was divided into seven sub-basins to examine the spatial distribution of flood events (Table 2). The most occurrences are recorded in the southern part of the study area and specifically in the drainage basin of Karditsa-Enippeas River. Moreover, increased clustering of flood events is observed in the drainage basin of Larisa-Karla.

Old topographic maps of 1881 and bibliographic data were used to examine the relation of the paleo-environment of Pinios River with the spatial distribution of floods. The old topographic maps show earlier marshy areas and lakes, which nowadays, have dried up. According to [36], paleo-lake systems were located in the drainage basin of Pinios River during the Quaternary. One paleo-lake was located in the western part of the study area. Another, lower altitude paleo-lake, which was not connected to the previous one, was in the eastern part. During the Quaternary, the weathering of the carbonate rocks led to the connection of independent paleo-environments, and the Pinios River emerged. There is a certain amount of clustering of flood events in the areas of former marshes, lakes and paleo-lakes, which are located in the drainage basins of Trikala-Neochoritis River, Larisa-Karla and Karditsa-Enippeas (Figure 12).

In the study area, drainage and irrigation works along with levees were built 90 years ago, to protect it from flooding. However, incidents of annual flood events are reported in specific areas, and still, floods remain a big problem in the region [58]. Figure 12 presents the sites with earlier marshes

and lakes, which are located in the lowest and flattest parts of the study area. The locations with earlier marshes and lakes are characterized by gentle slopes with low elevation and form the bottom of the drainage basin of the drainage network. These areas are developed mainly over semi-permeable or impermeable formations. The underlying geology and soil type affect the quality and the rate of infiltration. More specifically, the sediments within these wetlands may consist of materials such as clays and fine sands. These materials create a saturated layer on the surface during intense rainfall. Therefore, the specific runoff may be unable to further infiltrate the surface water into the subsurface. Consequently, these specific locations are sites where surface water runoff naturally collects and can lead to floods. Although, the wetlands and lakes have dried up, the morphology of the relief in these areas has not changed, causing flooding. According to studies carried out in other similar regions [35,52], the drainage basins that were 15% covered by wetlands had a flooding supply of 60–65% more than those not bearing wetlands. Consequently, drying up of the marshy areas and lakes of the region has strongly favored the flood events in the drainage basin of Pinios River.

The low and flat areas are located in the central, south-eastern and coastal part of the study area (Figure 12), and they contain many flood events. Usually, when a flood occurs, the lowest and flattest areas will be flooded first. Thus, these locations are areas prone to flood. This lowland morphology of the plain and the narrow passages along the river course such as Kalamaki, Rodia and Tempi gorges are the main reasons for the flooding [58]. According to [55], other reasons favoring the flood genesis due to human activities are some bridges with inadequate heights that are across the river and the construction by the farmers of "handy" barriers in the river channel for storage of irrigation water.

Endogenic processes such as active tectonics of the study area probably will not affect the flooding hazard in the near future. On the other hand, exogenic processes such as climate changes are capable of influencing future flood occurrences. Possible climate changes will cause sea level rise. Different scenarios were produced by the Intergovernmental Panel on Climate Change (IPCC), which predicts sea level rise from 0.3 to about 1.0 m until 2100 [59]. The possible sea level rise will increase the flood risk in the coastal area of the study area.

The historic flood data and old topographic maps provide valuable information for land use planning at a regional scale, leading to the determination of the safe and non-safe areas for urban activities [60]. The proposed methodology estimates the localization of sites prone to flood, and it may be utilized for flood hazard assessment mapping and for flood risk management. In areas prone to flooding, the appropriate land use planning along with the selection of the proper constructions are essential to prevent and mitigate the consequences of flood hazard occurrence. Consequently, proper land use for specific areas (i.e., parks, residential areas, etc.) may be determined. Additionally, construction of flood control works such as dams and levees, ponds, lakes and lagoons may be selected. Thus, planners, engineers and policy makers may use the applied approach for new and existing land use planning projects and for floodplain management.

**Author Contributions:** G.D.B conceived of the research. H.D.S. and G.D.B. designed the research and the data analysis. K.S. prepared and analyzed the data. G.D.B., H.D.S., K.S. and E.K completed the field work. H.D.S. and E.K. created the figures. G.D.B and H.D.S. wrote the paper.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Flood Hazard Mapping of a Rapidly Urbanizing City in the Foothills (Birendranagar, Surkhet) of Nepal**

**Sushila Rijal 1,\*, Bhagawat Rimal 2,3 and Sean Sloan <sup>4</sup>**


Received: 26 March 2018; Accepted: 3 May 2018; Published: 5 May 2018

**Abstract:** Flooding in the rapidly urbanizing city of Birendranagar, Nepal has been intensifying, culminating in massive loss of life and property during July and August 2014. No previous studies have monitored underlying land-cover dynamics and flood hazards for the area. This study described spatiotemporal urbanization dynamics and associated land-use/land-cover (LULC) changes of the city using Landsat imagery classifications for five periods between 1989 and 2016 (1989–1996, 1996–2001, 2001–2011, 2011–2016). Areas with high flood-hazard risk were also identified on the basis of field surveys, literature, and the Landsat analysis. The major LULC changes observed were the rapid expansion of urban cover and the gradual decline of cultivated lands. The urban area expanded nearly by 700%, from 85 ha in 1989 to 656 ha in 2016, with an average annual growth rate of 23.99%. Cultivated land declined simultaneously by 12%, from 7005 ha to 6205 ha. The loss of forest cover also contributed significantly to increased flood hazard. Steep topography, excessive land utilization, fragile physiographic structure, and intense monsoonal precipitation aggravate hazards locally. As in Nepal generally, the sustainable development of the Birendranagar area has been jeopardized by a disregard for integrated flood-hazard mapping, accounting for historical land-cover changes. This study provides essential input information for improved urban-area planning in this regard.

**Keywords:** urbanization; flood; remote sensing/GIS; Birendranagar; Nepal

#### **1. Introduction**

Hazards, defined as natural or human-induced activities that elevate the probability of material, social, or natural loss [1], are typified by the nexus of uncontrolled urbanization in contexts susceptible to natural flood, landslides, and earthquakes in Nepal. Here, urbanization is understood as processes leading toward increased population density, socioeconomic activities, and expanded built up areas and associated infrastructures [2]. Natural hazards and urbanization can interact to amplify land-use change, such as that negatively affecting agricultural areas [3,4]. In hazardous areas, including newly formed urban areas, land management and planning that would enhance resilience depend therefore on an understanding of those land-use/land-cover change (LULC) patterns that accentuate latent hazards [5].

Floods are the most common and devastating natural hazards [6] on a global scale and have been increasingly frequent and devastating since the mid-20th century [4]. Of all flood events recorded between 1950 and 2011, most have occurred during recent decades. Some 2% occurred during the 1950s, rising rapidly each decade thereafter to 3.9%, 6.6%, 13.2%, 21.9%, and 52.2% for the 1960s, 1970s, 1980s, 1990s, and 2000s, respectively [7]. Flood events also accounted more than two-third of all hydrological disasters during 2012–2014 [8] and 90.9% during 2015–2016 [9]. Flood hazards are global occurrences, yet associated economic damage and fatalities concentrate disproportionately among the continents [7]. Due to favorable geomorphological, metrological, and anthropogenic factors [10], more than 60% of the total economic and human losses during the period 1950 to 2011 were concentrated in Asia [7]. South Asia accounts for 33% of all Asian floods, 50% of associated fatalities, and 38% of the effected regions. As a proportion of South Asian totals, Nepal accounts for 7.2% of fatalities, 7.4% of the total victims, and 3.1% of economic losses, thereby ranking the country after India, Bangladesh, Pakistan, and Afghanistan [11]—8th globally in terms of flood fatalities [12] and 30th globally in terms of flood-hazard risk [13]. Of the five million Nepalese effected by natural hazards between 1971 and 2007, 68% were flood-related [14]. Between 1982 and 2014, nearly 9000 Nepalese lost their lives in flood and landslide events [15,16]. The Koshi flood event during 2008 affected 65,000 people and 700 ha of fertile land in the eastern lowland Tarai region [17]. In August 2017, floods affected 35 districts and destroyed 190,000 homes, fully or partially, with total economic losses estimated at \$584.7 million [18]. While flooding in Nepal is triggered by monsoonal rainfall, steep and erosive topography, and wide catchments [4,19], its frequency and intensity increased largely due to increasing anthropogenic factors, namely, improper land-use, poorly-planned urbanization, deforestation, and settlements along river banks [20].

The study area, Birendranagar city of Nepal's Surkhet district, is rapidly urbanizing due to high rates of migration. Locals have migrated to the city in the quest of better quality of life and opportunities in the formal and informal markets [21], resulting in a 720% urban population growth rate for the period 1982–2015 (Table 1). It is a major socioeconomic hub and administrative center and an important gateway to Karnali zone, as well as a migrant-receiving area, mainly from the Dailekh district and the Karnali zone. Urban development has been haphazard due to deficient urban plans/policies and weak implementation. New and expanding urban settlements are characterized as spatially dispersed or discontiguous, and frequently arise over prime farm lands surrounding historic Birendranagar city. These urban areas, collectively defining Birendranagar city, are frequently devastated by annual flood events affecting transportation networks, and other infrastructure (e.g., bridges, markets), livelihoods, fatalities (numbering 38 in 2014 [16]), and soil erosion and transport to neighboring districts (27,092.94 Mt) of soil annually [22]. Despite their frequency, or perhaps because of their frequency, such events struggle to sustain the attention of policymakers [4]. Resilience is a vital tool with which to reduce the vulnerability [23]; however, the integration of flood-hazard risk management in regional plans and policies is hindered due to the lack of routine-based researches. Urban planning and flood hazard management to address these issues is inadequate, in part due to the lack of reliable, updated spatial databases.



To this end, this study uses remote sensing and GIS techniques to observe LULC changes and identify flood-risk hazards underlying the urbanization of Birendranagar city. This research will remain as an important benchmark for Nepalese planners/policymakers and land-change researchers, since its insights and outputs may serve as essential inputs for sustainable land-use plans and strategies for flood-hazard mitigation.

#### **2. Method**

#### *2.1. Study Area*

Nepal, like the focal region of this study, is increasingly urban. Nationally, the urban population grew from 2.9% of the general population in 1952/54 to 17.1% in 2011 [24] to more than 50% by 2017 [25]. Urban centers similarly increased in number from 10 to 58 over this period, and to 292 following after the local level reconstruction in 2017. Urbanization has been driven largely by internal migration, in turn sustained by regionally unequal development and economic opportunities [21,26]. Nepal experienced decade long political armed conflict during 1996–2006 [27,28], which displaced or otherwise spurred the migration of rural dwellers to growing cities in search of security. As the conflict subsided by 2006, several development activities were advanced throughout the country [29], again concentrating largely in select urban areas. Coincidentally, a peri-urban land market boom resulted in rapid, disorganized settlement expansion at the expense of arable lands [21,29].

The study area, Birendranagar city is the capital city of Karnali Province (Figure 1). After the reclassification of local administrative units as mandated by the New Constitution of Nepal, 2015, the former Village Development Committees (VDCs) were consolidated as Birendranagar city. This city spans 24,582 ha area and is divided amongst 16 wards housing 100,458 residents as of 2015 (Table 1) [24].

**Figure 1.** Location Map of the study area.

The Birendranagar city spans Nepal's sub-Himalayan and lesser-Himalayan zones, which are characterized by a warm-moist temperate, hot-dry sub-tropical, warm-dry sub-tropical, and cool-moist temperate climates, with annual temperatures ranging from 10 ◦C to 30 ◦C [22]. Elevation ranges from 380 to 2224 MASL. Birendranagar city receives intensive monsoonal precipitation, mainly between June and September, causing flooding and soil erosion, including along river banks. The study region is composed of fluvio-lacustrine sediments (sand, silt, clay, cobbles, and pebbles) deposited from the northern and southern parts of Siwalik Hills. Hard rocks in the region are comprised of sedimentary, meta-sedimentary, and metamorphic rock.

#### *2.2. Extraction of LULC Change*

Remote Sensing (RS) and Geographic information System (GIS) are successfully applied in various fields. They are effectively used for the LULC analyses [4,30–34], and flood area mapping. Multi-temporal time series of Landsat images [35] is capable of exploring accurate [36] LULC dynamics in specific time and space [21].

A city level land-use/land-cover (LULC) change analysis was realized by classifying time-series Landsat satellite imagery for 1989, 1996, 2001, 2006, 2011, and 2016. Atmospherically-corrected and maximal cloud-free Landsat images were sourced from United States Geological Survey (USGS) data portal (https://earthexplorer.usgs.gov). All scenes were verified for geometric accuracy, and data were projected on to UTM 44N (WGS 1984). For land-cover classification, six images with path/row 144/40 were selected for analysis. A "region of interest" (ROI) boundary representing municipality level study area was delineated for the classification. ENVI software was used to stack, subset, and classify the Landsat images via a maximum-likelihood (ML) algorithm [21,37,38]. A supervised approach with the ML-classifier algorithm was applied for the extraction of LULC. The land-cover classification system as recommended by Anderson et al. 1976 [39] was taken into consideration to classify the LULC and mainly seven LULC classes: urban/built-up, cultivated land, forest, shrub, sand, water, and others were identified (Table 2). Google Earth images (http://earth.google.com) and topographical maps from Survey Department, Government of Nepal (scale 1: 50,000) [40] were used to verify the LULC results. LULC transitions between each observation year were assessed using land change modeler (LCM) of TerrSet software developed by Clark Lab (www.clarklabs.org). Socioeconomic information on the drivers of urbanization were obtained from Central Bureau of Statistics [24] and use of Surkhet district profile [22].


**Table 2.** Land-cover classification scheme.

#### *2.3. Identification of Flood Hazards*

Areas of relatively acute flood hazard were further identified through field verification and informal surveys of local people during a field visit in 2014. High-hazard areas were geo-located using a GPS, and post-hazard observations were realized using high-resolution images of Google Earth. Contextual information was derived from the disaster reports of the National Planning Commission (NPC) [41] and the Ministry of Home Affairs (MoHA) [13,18], with additional information were garnered from the District Development Committee (DDC) [22]. Precipitation data for Birendranagar station were sourced from Department of Hydrology and Meteorology [42].

#### *2.4. Urban/Built-Up Area Expansion Rate*

Total urban area was estimated each observation year to measure the rate of urban expansion [21,37]. The rate of urban expansion describes the average annual growth of urban area during a given period, thus

$$\text{BER} = (B\_2 - B\_1) / (T\_2 - T\_1) \* 100\tag{1}$$

in which BER implies the urban expansion rate (*ha/year*), and *B*1, *B*<sup>2</sup> refers to the urban area (*ha*) between times *T*1, *T*<sup>2</sup> in years.

#### *2.5. Classiifcation Accuracy Assessment*

The verification of each Landsat classifications was realized with reference to high-resolution satellite imagery in Google Earth (http://earth.google.com), 1:50,000 scale topographical maps from the Nepalese Survey Department [40], and GPS points collected during the field visit, with the mixture of these reference data varying by classification year. Over the all-time periods of 1989–2016, some 350 sample points were interpreted across the city. Points were distributed via stratified random sampling method, with at least 50 points for each land-use/cover class. The overall accuracy, user's accuracies, and producer's accuracies were obtained for each observation year [21].

#### **3. Result and Discussion**

#### *3.1. Land-Use Land-Cover Change between 1989 and 2016*

The overall accuracies obtained for respective years were 83%, 82%, 85%, 85%, 84%, and 86% in the city. Over the 27 years of observation, the predominant LULC changes were (a) the rapid increase in urban cover after 2001, and more gradual increase in shrub lands; (b) the simultaneous losses of cultivated lands, as well as the steady but lesser decline in forest cover (Table 3, Figures 2 and 3).

During 1989–2016, urban area increased 571 ha, from 85 ha to 656 ha, with an average annual growth rate of 23.99% (Figure 2). Virtually all urban growth occurred after 2001; urban area was steady until 2001, but more than doubled between 2001 and 2006, from 113 ha to 289 ha. The decrease in cultivated lands paralleled this growth of urban areas in both timing and areas (Figure 2). Cultivated lands have declined 800 ha since 1989, from 7005 ha to 6205 ha. The increases in sandy areas and water after 2011 are mainly associated with the major flood events during 2013 and 2014.

Cultivated land was the major source of the newly expanded urban area. During 1989–1996, with the annual urban growth rate of 0.67%, urban area increased slightly by 4 ha, of which 75% (3 ha) was previously cultivated land. Subsequently, during the period 1996–2001, 24 ha new urban area was established (5.4% growth rate), 98% of which was cultivated previously (Figure 3). Between 2001 and 2006, all 176 ha of new urban cover occurred over cultivated lands (31% growth rate). The period 2006–2011 experienced unprecedented 295 ha of urban growth (20% growth rate, with 97% sourced from cultivated). This conversion rate continued for the last period, 2011–2016, despite its urban growth rate dropping to only 72 ha (2.5%).

Similar to the replacement of cultivated lands by expanding urban areas, forests have been steadily replaced by expanding shrub lands as economic development proceeded, aggravating flood hazard. The inflow of migrants to Birendranagar city and encroachments upon forest areas accelerated after the eradication program of malaria regionally in 1958. Forests were felled to supply resources for the development of Indian railway line, and subsequent development activities widely provoked deforestation and forest degeneration [43]. Also, devastating practices of deforestation under the Rana regime (1846–1951) were continued under the subsequent Panchayat system (1960–1990). More recently, the protection of forest cover to safeguard environmental integrity and ecological functions such as hydrological flow and flood protection has been prioritized in vulnerable regions like Birendranagar. The national government has launched various community-based forest management plans and

President Chure-Tarai conservation program to maintain current forest cover. These efforts likely contributed to the cessation of forest loss after 2006 in Birendranagar following earlier losses (Figure 2).


**Table 3.** LULC change of Birendranagar city during 1989–2016 (Area in ha).

**Figure 2.** LULC change trend during 1989–2016. (**a**) Urban/built-up, (**b**) cultivated Land, (**c**) forest, (**d**) shrub (**e**) sand and water body, and (**f**) others.

**Figure 3.** LULC change Maps of 1989 to 2016.

Nepal's New Constitution of 2015 mandated the reconstruction of all administrative areas and their reclassification within the federal administrative hierarchy. Birendranagar was declared the capital of Karnali Province. Hence, its recent history of urbanization is expected to continue as new governmental investments in economic and infrastructure development attract additional migrants. A ring road around the city, currently under construction, is expected to enhance the industrial, commercial, and development activities and thus to provoke urbanization and LULC change widely in the near future [44].

#### *3.2. Flood Area Analysis*

The Siwalik (Churiya) range and the southern slopes of the Mahabharat range are highly prone to geo-disasters [14,45] due to their fragile geology, steep slopes, and intense monsoonal precipitation during June through September (Figure 4), causing regular landslides and debris flow along creeks and steep slopes [20]. The highest average monthly rainfall was 488.77 mm and 428.64 mm for July and August, respectively, between 2000 and 2015 [42]. The southern belt of Churiya region and mid-hills are especially subject to intense rainfall from monsoonal low-pressure systems entering from the Bay of Bengal. The topographical slope of the study area ranges from below 4 degrees in the valley plain to 48 degrees in the Mahabharat range (Figure 5). The rocky soils cannot absorb intense rainfalls, resulting in major overland runs-off carrying soil and debris that have caused significant economic and human losses as urban expansion and deforestation have proceeded [41]. New settlements and urban expansion along river banks have disturbed river channels and drainage system. Deforestation and sand/gravel mining on the Siwalik range and upstream river beds have aggravated landslide hazards during the dry seasons. The major erosive and destructive forces of swollen rivers, which locals claim have increased in recent decades, dissipate only once steep riverine channels give way to gentler slopes, frustrating potential geo-engineering solutions to flood disasters.

*Land* **2018**, *7*, 60

**Figure 4.** Trend of monthly mean rainfall of Birendranagar station during 2000–2015.

**Figure 5.** Slope Map of the Study Area.

Surkhet district is affected annually by monsoonal flood events, the impacts of which have been increasingly devastating in recent decades. During August, 2014, more than 12,385 people from 2327 families were displaced [46], 1581 houses were washed away, 301 houses were damaged, 15 government schools were destroyed, and 31 more were damaged, resulting in the economic loss of NRs 6 billion. Additionally, 411 irrigation projects, 99 drinking water schemes, 11 child development centers, and 663 ha of forest were swept away [43]. Latikoili, Ghatgaun, Satakhani, Chinchhu, Lekhparajul, Hariharpur, Babiyachaur, Tatapani, Taranga, Dharapani, and Kunathari settlements were the most affected areas [46]. Birendranagar city was particularly highly affected by the flooding [44] on Khorke River (Figure 6) and Itram River (Figure 7). This 2014 flood not only impacted the settlements and cultivated area, but also claimed human lives.

**Figure 6.** River bank erosion and sediment transportation by Khokre River (photo taken by author, 2014).

**Figure 7.** Affected infrastructures from the floods in the Itram River (photo taken by author, 2014).

Our survey of flooding events highlighted that the floods along the Neware river, Gagre river, Geruwa river, Tuni river, Dwari river, and Dundra river also effected the nearby settlements and farm lands. This study has identified several settlements and cultivated lands along these river banks that are at high risk of flood hazards, which are presented in Figure 8.

**Figure 8.** Flood hazard map.

#### **4. Conclusions**

The study describes the LULC changes and urbanization process surrounding Birendranagar city during 1989–2016 using Landsat imagery and extensive consultations with local residents and officials. It highlights the rapid expansion of urban/built-up area from 2001. Urban areas have expanded most aggressively over fertile farm lands and become agglomerated within the flood-prone valley, particularly along Surkhet-Jumla-Karnali Highway and cross sectional roads. Forest conservation programs have helped slow deforestation during the last two decades. While wider reforestation measures are probably required to counter soil erosion and downstream transport during monsoons, this alone is probably insufficient to reduce flood hazard to reasonable levels. New settlements and cultivated lands established along river banks have especially high monsoonal flood hazards. Therefore, careful sustainable urban planning, migration control, and redirection measures, as well as flood hazard management and monitoring programs, are also required. Flood hazard management and monitoring strategies are not yet explicitly incorporated into local urban development plans, let alone regional land-use plans [47]. Indeed, much recent urban development has been informal or otherwise not subject to specific plans and zoning. Addressing the challenges of flood hazard mitigation in Nepal is thus fundamentally a challenge of instituting good governance practices based on solid empirical foundations, whereby development is subject to local zoning plans and laws.

**Author Contributions:** S.R. and B.R. designed the project and wrote the manuscript. B.R. contributed to remote sensing and GIS part including field data acquisition and verification. S.S. contributed to writing and English editing. All authors revised and edited the final manuscript.

**Acknowledgments:** The authors thank all the scientists and local participants who have participated in the establishment of database.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Physical and Anthropogenic Factors Related to Landslide Activity in the Northern Peloponnese, Greece**

#### **Hariklia D. Skilodimou 1, George D. Bathrellos 1,\*, Efterpi Koskeridou 2, Konstantinos Soukis <sup>3</sup> and Dimitrios Rozos <sup>4</sup>**


Received: 27 June 2018; Accepted: 15 July 2018; Published: 19 July 2018

**Abstract:** The geological, geomorphic conditions of a mountainous environment along with precipitation and human activities influence landslide occurrences. In many cases, their relation to landslide events is not well defined. The scope of the present study is to identify the influence of physical and anthropogenic factors in landslide activity. The study area is a mountainous part of the northern Peloponnesus in southern Greece. The existing landslides, lithology, slope angle, rainfall, two types of road network (highway-provincial roads and rural roads) along with land use of the study area are taken into consideration. Each physical and anthropogenic factor is further divided into sub-categories. Statistical analysis of landslide frequency and density, as well as frequency and density ratios, are applied and combined with a geographic information system (GIS) to evaluate the collected data and determine the relationship between physical and anthropogenic factors and landslide activity. The results prove that Plio-Pleistocene fine-grained sediments and flysch, relatively steep slopes (15◦–30◦) and a rise in the amount of rainfall increase landslide frequency and density. Additionally, Plio-Pleistocene fine-grained sediments and flysch, as well as schist chert formations, moderate (5◦–15◦) and relatively steep slopes (15◦–30◦), along with the amount of rainfall of >700 mm are strongly associated with landslide occurrences. The frequency and magnitude of landslides increase in close proximity to roads. Their maximum values are observed within the 50 m buffer zone. This corresponds to a 100 m wide zone along with any type of road corridors, increasing landslide occurrences. In addition, a buffer zone of 75 m or 150 m wide zone along highway and provincial roads, as well as a buffer zone of 100 m or 200 m wide zones along rural roads, are strongly correlated with landslide events. The extensive cultivated land of the study area is strongly related to landslide activity. By contrast, urban areas are poorly related to landslides, because most of them are located in the northern coastal part of the study area where landslides are limited. The results provide information on physical and anthropogenic factors characterizing landslide events in the study area. The applied methodology rapidly estimates areas prone to landslides and it may be utilized for landslide hazard assessment mapping as well as for new and existing land use planning projects.

**Keywords:** landslides; geographic information system (GIS); frequency ratio; density ratio; human activities; land use planning

#### **1. Introduction**

Landslides are physical phenomena, active in geological time. They have affected the natural environment and existing biota since even before the appearance of man on Earth. Nowadays, they are considered as one of most significant natural hazards worldwide. Frequently, their associated consequences have adverse long-term effects. When these consequences have a major impact on human life and activities, they become natural disasters [1].

Mass movements are identified as the movement of a mass of soil, rock, debris or earth down a slope [2]. They can include falls, topples, slides, spreads and flows. These phenomena are part of the process of hill slope erosion which is responsible for the introduction of sediment into streams, rivers, lakes, reservoirs and finally the oceans [3].

Landslides can be triggered by a variety of external factors, such as intense rainfall, earthquakes and rapid stream erosion [3–9]. Additionally, human activities such as deforestation of slopes, removal of slope support in road cuts, alteration of surface runoff paths, have become important triggers for landslide manifestation [10,11].

Every year, landslides kill people and cause huge property damage in mountainous areas of the world [11]. Kumar et al. [12] reported that in the Himalayas over 2000 landslides killed more than 5000 people during 2013. In the USA, landslide events cause an estimated US\$1–2 billion in economic losses and about 25–50 fatalities annually [13]. In Italy, at least 263 people were killed by mass movements during the period of 1990–1999 [14].

Much of Greece consists of hilly and mountainous terrain subject to landslide manifestation. Therefore, landslide events are common phenomena and cause significant damage to road networks, sections of urban areas and cultivated land [15–17].

Several physical process factors and human activities influence landslide activity. Landslides in the future will most likely occur under geomorphic, geologic, and hydrologic conditions that have produced past and present landslides. Types of bedrock, slope steepness, and precipitation zones represent, respectively, geologic, geomorphic and hydrologic factors. Weak, incompetent rock is more likely to fail than strong, competent rock; general steeper slopes have a greater chance of landsliding, while rainfall is considered as an important factor in slope stability, almost as important as gravity [2,3,11,15]. In addition, the frequency of landslide events commonly increases, when man-made structures induce changes in mountainous environments [3]. Human settlements, cut-and-fill construction for roads, the construction of buildings and railroads, changes in land cover, and terracing for agriculture contribute to the conditions that lead to slope failure [18,19].

However, in most cases the influence of physical and anthropogenic factors in landslide manifestations is not well known. Determination of the actual landslide zone induced by geologic, geomorphic parameters and rainfall, along with human activities such as the road network, urban and cultivated areas is an important tool for engineers, planners, and environmental managers. This procedure is very useful to identify areas prone to landslide and assess landslide hazard, and it is also a necessary step for land use and government urban planning policies worldwide [20–23]. Moreover, the determination of the landslide zone of disturbance is vital for road engineers who must deal with the costly and sometimes life-threatening problems caused by road-induced landslides. Realistic assessment of the impact of proposed transportation construction must be quantifiable as road costs may be significantly increased by landsliding during construction [12,18,24].

During the last decades geographical information systems (GIS) and Earth observation (EO) data have become integral tools for the evaluation of natural hazard events. These current geospatial technologies are very useful for assessing future hazard occurrences and identifying the vulnerability of communities to hazards. In this sense, GIS is an excellent tool in the spatial analysis of multi-dimensional phenomena such as landslides [25–30].

The scope of the present study is to identify the influence of physical and anthropogenic factors on landslide occurrences. To accomplish this, existing landslides, lithology, slope angle and rainfall in conjunction with two types of road network and land use were taken into consideration. Statistical analysis and GIS are applied to process and evaluate the landslides and factors. Thus, the spatial distribution of landslide frequency, magnitude and the association of lithology, slope, rainfall, roads and land use with landslide occurrences in mountainous terrain were determined. The case study area was a mountainous part of the northern Peloponnesus in southern Greece.

#### **2. Study Area**

The study area is located in the northern part of the Peloponnesus, in southern Greece (Figure 1a). The region covers an area of about 194 km2 and its altitude varies from 0 to 1400 m. The morphology of the area comprises southern mountainous land with very steep slopes reaching an altitude of 1400 m a.s.l.; intermediate semi-mountainous land with lower altitudes and steep slopes; and a northern coastal area of limited extent with low altitudes and gentle slopes. The drainage networks of the area flow with a main direction from SW to NE and discharge into the Gulf of Corinth (Figure 1b).

**Figure 1.** (**a**) Location map of the study area; (**b**) the elevations of the study area, the drainage network, the road network and the main settlements.

Climatologically, the study area is classified as Mediterranean, because of its coastal extension, without considerable temperature variations. The annual precipitation is about 800 mm in the mountainous part of the study area and about 550 mm at the lowlands.

The study area is located on the southern part of the Gulf of Corinth, and its landscape evolution is controlled by the neotectonic action of the graben, which forms the gulf. The presence of Pleistocene marine terraces in the northern part of the area proves that it has been affected until today by significant neotectonic uplift. This fact has been clearly expressed by the intense seismicity of the Gulf of Corinth [31].

The geological formations that crop out in the study area comprise both Neocene–Quaternary deposits and alpine formations. Alpine formations from two Hellenic geotectonic zones (Olonos–Pindos and Gavrovo–Tripolis) participate in the southern part of the study area. The study area is constructed by: (a) Holocene alluvial deposits, talus cones and scree; (b) Plio-Pleistocene sediments consisting of conglomerates such as clayey marls, marls, silty sands and sandstones; (c) schist-chert formations (Upper Cretaceous–Eocene); (d) Cretaceous limestones; and (e) schist-chert formations (Jurassic–Lower Cretaceous) [32].

Therefore, the study area mainly consists of Neocene deposits and, as a part of the Corinthian graben, it is characterized by intense neotectonic activity. Consequently, it is affected by many landslide events. It is an appropriate case study area because it comprises small urban areas, where there is still a potentiality of future urban growth and planning. Moreover the study area has often suffered from the consequences of several landslides, which usually cause serious damage in inhabited areas, road networks, and cultivated areas [33].

#### **3. Materials and Methods**

This study was carried out using the following:


A spatial database was created, and ArcGIS 10.0 software was used to process the collected data.

#### *3.1. Landslide Inventory Map*

The landslide inventory map compilation has involved the following steps: (a) landslides were recorded from previous works [33]; (b) landslide locations were recognized on satellite images; (c) the manifestations of landslide were verified and mapped by field work. The main scarp of every recorded landslide during the field work was depicted in topographic maps at a proper scale and then digitized as a polygon layer. According to Yilmaz [34], the scarp sampling strategy gives better results than the point one. The landslides were used for the compilation of the landslide inventory map (Figure 2). A number of 270 sites of landslide manifestation were examined throughout the study area, with a varied size from 4130 m<sup>2</sup> to 91,000 m2, having affected a total area of about 6.5 km2.

**Figure 2.** The landslide inventory map of the study area.

#### *3.2. Physical and Anthropogenic Factors*

The occurrence of landslides is largely a function of the interaction of several factors such as lithology, tectonic settings, geomorphological settings, earthquake, rainfall and human activities. These factors, which are directly or indirectly related with the occurrences of landslides, are commonly known as landslide-related factors [35–39]. In landslide susceptibility mapping studies, it is believed that the accuracy of the results increase when numerous landslide-related factors are included in the analytical process [40]. In many cases, this is usually difficult to do, because detailed data is hard to find. The aim of the present work is to determine the influence of physical and anthropogenic factors on landslide occurrences rather than the production of a landslide susceptibility map. Moreover, data for all the aforementioned factors was not available. For these reasons, analyses in this study depend on the physical factors which are lithology, slope angle and rainfall, while the anthropogenic factors are the road network and land use.

Each factor is separated into sub-classes. The determination of the classes' number as well as their boundary values were based on: (a) the literature review [3,11,15–21,33–39]; (b) the extended field observations in the framework of this study; and (c) personal experience from previous studies.

#### 3.2.1. Lithology

Lithology is the main predisposing factor controlling landslide development. Different lithological units have different engineering geological behaviors and they are very important in providing data for landslide-related factors studies [3,36,39,40]. For the study area the geological formations were digitized and unified according to their engineering geological behavior [20,32,33,35], in relation to landslide manifestation. Thus, lithology includes five classes as follows: (a) fine, fine-coarse to coarse and loose to semi-coherent Quaternary formations; (b) cyclothematic formations (Plio-Pleistocene fine-grained sediments and Flysch sediments); (c) Plio-Pleistocene coarse-grained sediments; (d) thin-bedded schist chert formations; and (e) moderate to thick bedded limestones (Figure 3).

**Figure 3.** The lithology of the study area consisting of: Quaternary formations (fine, fine-coarse to coarse and loose to semi-coherent sediments), Cyclothematic formations (Plio-Pleistocene fine-grained sediments and Flysch sediments), Plio-Pleistocene coarse-grained sediments, thin bedded schist chert formations, and moderate to thick bedded limestones.

#### 3.2.2. Slope Angle

The slope angle has an effect on slope stability, increasing the landslide hazard [35,39,40]. Contours with 20 m intervals and height points were digitized from topographic map (scale 1:50,000) and saved as line and point layer correspondingly. A digital elevation model (DEM) was derived

from the digitized elevation data using the 3D Analyst extension of ArcGIS, and the slope layer was extracted from DEM. The slopes were classified into five classes: (a) <5◦, (b) 5◦–15◦, (c) 15◦–30◦, (d) 30◦–45◦, and (e) >45◦ (Figure 4).

**Figure 4.** Map showing the spatial distribution of slopes of the study area.

#### 3.2.3. Rainfall

Precipitation is a triggering factor for landslide occurrences [35,39]. The mean annual rainfall of the area varies between 550.7 mm and 789.8 mm. The intensity of rainfall was not analyzed due to lack of data. For the necessities of this study, the precipitation map was produced, using the data of the main meteorological stations in the area and applying the Inverse distance weighted (IDW) interpolation method. The distribution of representative rainfall depends upon the spatial distribution of the stations and elevation [41]. The stations used are uniformly distributed in the study area both hypsometrically and territorially giving accurate precipitation distribution.The precipitation map was separated into three classes, i.e.,: (a) <600 mm, (b) 600 mm–700 mm, and (c) >700 mm (Figure 5).

**Figure 5.** Map showing the spatial distribution of rainfall of the study area.

#### 3.2.4. Road Network

The artificial and natural parts of the slopes around a road network are more sensitive in landslide manifestations [35]. The road network of the study area was digitized as a polyline layer using the topographic map. The complete road network included about: 24 km of highway, 322 km of provincial roads, 389 km rural roads, and totaled about 735 km in length (Figure 6).

The proximity to the road network is often related to an increase of landslide occurrences, so disturbance zones were created around the road network of the area. The creation of disturbance zones is referred to using the GIS term 'buffer' zone along a road network using GIS software. The classes of the buffer zones used in this analysis were: 10, 25, 50, 75, 100, 150, 200, 300 m and >300 m (6). The buffer lengths were measured on both sides of road network and represent a zone that is twice as wide as the value listed above. Therefore, a 100 m buffer describes a zone or swath 200 m wide along the road. Two sets of the aforementioned buffer zones were created along the different road types: one for highways and provincial roads and the other for rural roads.

**Figure 6.** Buffer zones of increasing length along the highway, provincial, and rural roads of the study area.

#### 3.2.5. Land Use

The land use of the study area was taken from the CORINE 2012 Land Cover (CLC) map, Copernicus Program [42]. The program contains land cover data for Europe including land cover class description at scale 1:100,000 published by the European Commission. The CORINE land use map was classified as follows: (a) urban area, (b) cultivated area, (c) forest, (d) shrubby area, and (e) bare area (Figure 7). The land use of the area was saved as polygon layer.

**Figure 7.** The land use of the study area.

To predict landslides, it is necessary to assume that landslide occurrence is determined by landslide-related factors, and that future landslides will occur under the same conditions as past landslides [3,36]. For this reason, as a first step, the relative frequency of landslides was calculated. The relative frequency is given from Equation (1):

$$FL = Ln/\Sigma Lm \tag{1}$$

where *FL* = the relative frequency of landslides, *Ln* = the number of landslide that located in each category of the landslide-related factors, and *ΣLn* = the total landslide events of the study area.

The landslides were transformed from polygon layer to point layer to calculate the landslide frequency. Thus, the number of landslide events, which were located in the classes of each factor was then determined by using the Spatial Analyst Tools of ArcGIS 10.0.

The next step was the estimation of the landslide frequency by applying a frequency ratio statistical analysis. This approach is based on the relationships between spatial distribution of landslides and each class of the involved factors. According to Lee and Pradhan [36], the frequency ratio is the ratio of the probability of a landslide event to a non-event for a given zone. The frequency ratio is estimated by the following Formula (2):

$$FR = FL / AC \tag{2}$$

where *FR* = the frequency ratio, *FL* = the relative frequency of landslides expressed in percentages, and *AC* = the area ratio for each category to the total area in a percentage form.

Furthermore, to determine the magnitude of a landslide area, the relative density of landslides was computed. The relative density of landslides is expressed by the Equation (3):

$$DL = La / \Sigma La \tag{3}$$

where *DL* = the relative density of landslides, *La* = the landslide area within each class, and *ΣLa* = the total landslide area.

Consequently, the area of landslides within the classes of each factor was computed by using spatial analyst capabilities. Similarly to landslide frequency, a density ratio statistical analysis was applied to specify the relationship between landslide area and the landslide-related factors. The density ratio is calculated according to the following mathematical operator (4):

$$DR = DL / A\mathbb{C} \tag{4}$$

where *DR* = the density ratio, *DL* = the relative density of landslides expressed in percentages, and *AC* = the area ratio for each category to the total area in a percentage form.

#### **4. Results**

The statistical analysis described above allows discrimination of landslide frequency and magnitude in each category of the adopted factors. The physical and anthropogenic factors, their categories and the results of the statistical analysis are presented in Table 1.

Regarding the physical factors, Figure 8 shows the relative frequency (*FL*) and density (*DL*) of landslide distribution in each lithological formation of the study area. The maximum of *FL* (64%) and *DL* (61%) values are observed in the area underlain by cyclothematic formations consisting of Plio-Pleistocene fine-grained and flysch sediments. Additionally, high *FL* (24%) and *DL* (28%) values are presented in the Quaternary formations which consist of fine-grained to coarse-grained loose sediments. These values are almost two times lower in Quaternary formations than ones in cyclothematic formations. In Plio-Pleistocene coarse-grained sediments, the *FL* (10%) and *DL* (10%) values are six times lower than in cyclothematic formations. Schist chert formations have very low *FL* (1%) and *DL* (1%) values. The minimum of *FR* and *DL* values are observed in moderate to thick bedded limestones.

**Table 1.** The frequency (*FR*) and density (*DR*) ratio values of landslides into each category of physical and anthropogenic factors, A = area of each category in m2, *AC* = the area ratio for each category to the total area in a percentage form, *Ln* = number of landslide in each category, *FL* = the relative frequency of landslides expressed in percentage, *La* = the landslide area within each class, and *DL* = the relative density of landslides expressed in percentage.


Apart from the *FL* and *DL*, the landslide frequency ratio (*FR*) and density ratio (*DR*) were calculated. According to Lee and Pradhan [36] the FR value of 1 is an average value. Thus, ratio values greater than 1 indicate a strong relationship between landslides and the given factor, and ratio values smaller than 1 indicate a poor relationship between landslides and the given factor. As demonstrated in Table 1, the schist chert formations have the maximum *FR* value (7.1), indicating a strong relationship between this lithological formation and landslide occurrences. Similarly, in cyclothematic formations the *FR* value is 1.5, showing a strong association with landslide events. For the remaining lithological formations of the study area, the *FR* values are found to be less than one, indicating a poor relation with landslides.

**Figure 8.** Relative frequency (*FR*) and density (*DL*) values within each category of lithology in the study area. The *FR* and *DL* values are expressed in percentages.

Alike the frequency ratio, the greater the density ratio values above one, the stronger the relationship between landslide occurrence and the given factor; the lower the ratio values below one, the lesser the relationship between landslide occurrence and the given factor. The maximum value of *DR* (9.8) is observed in the schist chert formations, while in the cyclothematic formations and Quaternary formations the *DR* values are 1.4 and 1.0 correspondingly (Table 1). This fact proves a strong relationship between the aforementioned formations and landslide manifestation.

Concerning the slope angle of the study area, Figure 9 shows the *FL* and *DL* distribution in each category of slope. The majority of landslide events (50%) and magnitudes (52%) are located in the class of slopes 15◦–30◦. Very high *FL* (46%) and *DL* (44%) values are observed in the category of medium slopes 5◦–15◦. The minimum *FR* and *DL* values are located in the class of very steep slopes >45◦ (Table 1).

**Figure 9.** *FL* and *DL* values distribution in each class of slope angle. The *FR* and *DL* values are expressed in percentages.

The maximum *FR* and *DR* values are found to be 1.2 in the class of steep slopes 15◦–30◦. *FR* and *DR* values greater than one (1.1) have the class 5◦–15◦. Thus, these slopes are strongly related to landslide occurrences. The classes of slopes <5◦, 30◦–45◦and >45◦ have *FR* and *DR* values <1, indicating a poor relationship between them and landslides (Table 1).

Figure 10 illustrates the *FL* and *DL* values distribution in each category of rainfall. The maximum values of *FL* (57%) and DL (52%) are attributed in the category of 600 mm–700 mm. Moreover, very high values of *FL* (37%) and *DL* (43%) are observed in the category >700 mm.

**Figure 10.** *FL* and *DL* values' distribution in each class of rainfall. The *FR* and *DL* values are expressed in percentages.

The maximum *FR* values 1.1 and 1.0 represent, respectively, the categories of rainfall >700 mm and 600–700 mm. These two categories are strongly related to landslide manifestation. On the contrary, the maximum *DR* value (1.3) is observed in the category >700 mm, indicating a strong relationship between this category and landslide events (Table 1).

The applied statistical analysis identifies the distribution of landslide frequency and magnitude in nine zones of increasing length measured around to the road network. The analysis was performed for two different types of road network. Figure 11 shows the distribution of *FL* and *DL* values in within each buffer zone of the highway and provincial roads. The relative frequency and density of landslides are expressed in percentages.

**Figure 11.** *FL* and *DL* values within each buffer zone extending 0 to 300 m from highway and provincial roads and >300 m in the study area. The *FL* and *DL* values are expressed in percentages.

The *FL* reaches its maximum value within the first 50 m of distance from any given highway and provincial road and decreases in distances beyond 50 m from roads. The *FL* value is greatest (9%) within 50 m buffer zone and is relatively high as far as 100 m from roads. The value of *FL* within buffer zone of 50 m is about two times higher than in distances: (i) approximately 10 m (4%) and (ii) between 75 and 100 m (5%) from roads. Moreover, the *FL* values become relatively high in distance beyond 150 m (13%) from roads. This increase is not related to roads but is likely due to other factors influencing landslide manifestation (Figure 11 and Table 1).

As in the case of *FL*, the *DL* reaches its maximum value (9%) within the buffer zone of 50 m and shows a gradational decline with buffer length (Figure 11). The *DL* value is relatively high in distances between 50 and 100 m from roads. The *DL* value in area of 10 m buffer zone is about two times lower (4%) than in 50 m buffer zone and for area of 100 m buffer zone is one and a half times lower (6%) than in 50 m buffer zone.

In the buffer zones of 10 and 25 m, the *FR* values are 1.2 and 1.3 respectively, indicating a strong relationship between these zones and the occurrence of landslides. For the distances of 50 m and 75 m from roads, the *FR* values were found to be 1.1 and 1.0, respectively, showing a strong association with landslide events. The *FR* values are <1 in distance beyond 100 m, proving a poor relation with landslide occurrences (Table 1).

The maximum *DR* value (1.3) is found to be at distance of 50 m from roads. High *DR* values (1.1) are observed at distances of 10, 25 and 75 m from roads. Thus, these zones are strongly correlated with landslide occurrences. The ratio values in areas beyond 100 m of roads are <1, showing a low probability of landslide occurrences (Table 1).

With regard to the rural roads, the maximum *FL* value (18%) is attributed to a 50 m buffer zone and is relatively high as far as 100 m from roads (Figure 12). The *FL* value within a buffer zone of 50 m is about three times higher than in distance approximately 10 m (6%) and two times higher than in distance between 75 and 100 m (9%) from roads. Similarly to *FL* values, the *DL* values have its maximum within 50 m buffer zone (13%). The *FL* and *DL* values are increased in distance beyond 150 m from rural roads (Figure 12).

**Figure 12.** *FL* and *DL* values within each buffer zone extending 0 to 300 m from rural roads and >300 m in the study area. The *FL* and *DL* values are expressed in percentages.

The *FR* values reach their maximum (1.9) in 50 m buffer zone and its values are >1 in the area of the 10 m, 25 m, 75 m and 100 m buffer zones. Likewise, the *DR* values are >1 at distance 10 m, 25 m, 50 m, 75 m, and 100 m from rural roads. These observations indicate a strong relationship between landslides and the aforementioned distances (Table 1).

Concerning the land use of the study area, the calculation *FL* and *DL* provides an estimate of the influence of land use in landslide manifestation. The statistical analysis of *FL* shows that the vast majority of landslide events (70%) are located in cultivated areas, and a high value of *FL* within shrubby areas (21%). In urban and forest areas, it is seventeen and a half times lower (4% in each one) than in cultivated areas. The minimum value of *FL* is observed (1%) in a bare area (Figure 13 and Table 1).

**Figure 13.** *FR* and *DR* values within each category of land use in the study area. The *FR* and *DR* values are expressed in percentages.

The *DL* values follow the distribution of *FL* values in each category of land use. More specifically, the maximum value of *DL* (64%) is attributed to cultivated areas. A relatively high density value (28%) is observed in shrubby areas, while the minimum density (1%) value is calculated in bare areas. The *DL* values in urban and forest areas (4% in each one) are low (Figure 13 and Table 1).

The *FR* was found to be 1.4 in cultivated areas, which proves a high probability of landslide manifestation. In the rest of the categories of land use, the ratio values are <1, indicating a strong relation with landslide manifestation (Table 1). The highest *DR* value is 1.4 and refers to the cultivated areas, showing a strong association between this land use and landslide occurrences. The remaining classes of land use have ratio values <1, indicating a low probability of landslide manifestation (Table 1).

#### **5. Discussion**

Hydrological parameters such as precipitation and human activities are capable of changing landscape characteristics. On the other hand, the geological, geomorphic conditions of a mountainous environment remain unchanged or vary little from a human perspective. These factors influence where landslides occur [3,19–24]. In many cases, the association of physical factors and human activities in landslide occurrences is not well defined. According to Lee and Pradhan [36], statistical approaches based on the observed relationship between each factor and the spatial distribution of landslides is very useful to reveal the correlation between landslide locations and factors.

In the present study physical and anthropogenic factors are analyzed and evaluated to determine their association with landslide occurrences. Lithology, slope angle, rainfall, two types of road network and land use are correlated with the existing landslides by using *FL*, *DL*, *FR* and *DR* statistical analysis and GIS. The case study area is a mountainous part of the northern Peloponnese in Greece. The results provide information on the physical the anthropogenic factors characterizing landslide events in the study area.

The statistical analysis proves that the lithological formations of the study area influence the occurrence of landslides. The events and magnitude of landslide areas reach their maximum in the area underlain by cyclothematic formations consisting of Plio-Pleistocene fine-grained sediments and flysch (Figure 8). Additionally the *FR* and *DR* values are >1 in areas where cyclothematic sediments and shist-chert formations outcrop (Table 1). The latter is possibly due to the limited extent of shist chert formations in the study area. Plio-Pleistocene sediments are fine-grained with a variety of lithological horizons. They consist of clays, marls, alternating sands of a varying degree of diagenesis and/or their mixed phases. Flysch are strongly folded sediments because of the tectonic action (nappes and upthrusts), which in many places results to the formation of thick weathering mantle. Schist chert formations consisting of alternations of cherts, siltstones, thin plated limestones and sandstones while volcanic tuffs rarely participate at places. A thick weathering mantle is formed, mainly in the cases of the surface occurrence of siltstones. The aforementioned soils, hard soil to soft rocky sediments, are weak and incompetent formations and are prone to landslides [4,32,33]. Consequently, the cyclothematic sediments and schist chert formations of the study area are strongly related to landslide occurrences.

Regarding the slope angle, the maximum *FL* and *DL* values are observed in the category of relatively steep slopes 15◦–30◦. High values of *FL* and *DL* are found to be in the category of moderate slopes 5◦–15◦ (Figure 9). Additionally the *FR* and *DR* values are >1 in the classes of 5◦–15◦ and 15◦–30◦, indicating a strong relationship between them and landslide events. *FR* and *DR* values <1 are observed in the remaining classes of gentle slopes (<5◦–30◦), steep slopes (30◦–45◦) and very steep slopes (>45◦), showing a poor relationship between them and landslides (Table 1). According to Rozos et al. [33] this peculiar condition can be explained easily, as in nature slopes consisting of soil or hard soil to soft rocky formations (like those of the study area), and having high angles, fail almost immediately after their formation, resulting in lower slope angles. Finally the slopes with an inclination of around the angle of friction are those which fail after the action of triggering factors. On the other hand, rocky slopes are stable even in high angles suffering only from rock falls, wedge failures etc.

The *FL* and *DL* values increase with increasing in amount of precipitation (Figure 10). The maximum *FL* and *DL* values are attributed in the category of 600 mm–700 mm. FR and DR values >1 are observed in the category of rainfall >700 mm, showing a strong correlation with landslide events (Table 1). Therefore, the rainfall is an important factor in triggering manifestations of landslides [33,35].

Concerning the anthropogenic factors, nine buffer zones of increasing length measured around to two types of road network and five categories of land use are evaluated. The landslide events and magnitude of landslide areas are increased within a distance of approximately 0 to 50 m from any given road type (highway, provincial and rural road) and decreased in distances beyond 50 m from roads (Figures 11 and 12, Table 1). The relatively high values of *FL* and *DL*, which are measured in the 150 m buffer zone, are probably due to other landslide-related factors such as type of lithology or slope angle (Table 1). As already mentioned above, landslides of the study area tended to occur in locations with relatively steep slopes or consisting of cyclothematic sediments.

Consequently, the frequency and magnitude of landslides increase in close proximity to roads. According to Bathrellos et al. [35] this may be due to the fact that the road network sometimes destabilizes adjacent marginally balanced slopes, mainly by removing natural support for the upper part of the slope through undercutting the base of the slope during its construction and by adding extra weight on them. For the study area, landslide disturbance is associated with roads at distances of as much as 50 m from where they are constructed. Further away the frequency and magnitude of landslide occurrences gradually decreases. However, the 50 m buffer length, or 100 m swath, represents a sizeable area of the land surface. These results indicate that a 100 m-wide zone along road corridors increase the landslide activity.

When planning a new road, the alignments of the route should be carefully determined considering the probability of landslides during and after construction. The roads are likely to face unexpected problems of slope stability [12]. The quantification of landslide frequency, magnitude and distribution in mountainous terrain is an important consideration in calculating the cost of new roads construction and the maintenance of existing road networks [18]. Therefore, road engineers

and planners in similar mountainous environments may utilize the buffer zone of 50 m during the construction of a new road or to protect the existing road network from future landslide occurrences.

The *FR* and *DR* values (Table 1) show that buffer zones with distances smaller than/or equal to 75 m from highway and provincial roads are strongly associated with landslide manifestation. For buffer zones with distances >75 m, the ratio is <1, indicating a poor relation with landslide occurrences. Conversely, the *FR* and *DR* values (Table 1) are >1 at distances lower than/or equal to 100 m from rural roads, showing a strong relationship between landslides and these distances. In the study area, a 150 m wide zone along highway and provincial road corridors or a 200 m wide zone along rural road corridors indicates a strong association with landslides. This difference may be related to the fact that the rural road network is long and dense in the study area. Thus, the closer the distance is to the road network, the greater is the relationship with landslide manifestation.

In the case of land use, the majority of frequencies and magnitudes of the landslide area are located in cultivated areas (Figure 13). The *FR* and *DR* values are >1 only in cultivated land areas, showing a strong relation to slope failures (Table 1). The variations of the vegetation in an area constitute an important parameter affecting the slope failures, as slope stability is very sensitive to changes in vegetation [33,39]. The soil cohesion is modified depending on the type of vegetation and, thus, cultivated or sparsely vegetated areas are more prone to landslide processes. Additionally, crops can increase the moisture in the soil and alter ground water conditions. In the study area, crops cover an area of about 90 km2 and this represents 42% of the total area (Table 1). The extensive cultivated land combined with the altered ground water conditions is capable of causing landslide problems.

By contrast, the *FL* and *DL* values are limited in urban areas (Figure 13). Furthermore, the *FR* and *DR* values indicate a poor relationship between urban areas and landslide occurrences (Table 1). Most of the urban areas are located in the northern coastal part of the study area where there are few landslides.

#### **6. Conclusions**

In the present study, statistical analysis and GIS was applied to determine the relation of physical and anthropogenic parameters with landslide activity in a mountainous terrain.

Concerning the physical factors, the lithology of the study area controls landslide activity. The cyclothematic formations consisting of Plio-Pleistocene fine-grained sediments and flysch increase *FL* and *DL*. These sediments along with schist chert formations are strongly associated with landslide occurrences. Since *FL* and *DL* increase in relatively steep slopes, geomorphologic factors such as slope angle influence landslide occurrences. Moderate (5◦–15◦) and relatively steep slopes (15◦–30◦) are strongly associated with landslide events. Rainfall is an important factor in triggering the manifestation of landslides, as their frequency and density rise with an increase in the amount of rainfall. An amount of rainfall of >700 m is strongly related to landslide events.

In terms of the anthropogenic factors, the statistical analysis proves that the zone of landslide disturbance associated with roads is extensive. The frequency and magnitude of landslides increase in close proximity to roads. The maximum *FL* and *DL* values are observed within the 50 m buffer zone. This zone (100 m wide) along with any type of road corridors increase landslide occurrences. On the other hand, the buffer zone of 75 m or a 150 m wide zone along highway and provincial road corridors is strongly related to landslide manifestation. Since the rural road network is long and dense in the study area, the buffer zone of 100 m or a 200 m wide zone along rural roads is strongly connected with landslide events. The extensive cultivated land of the study area leads to an increase of *FL* and *DL*. This land use is strongly associated with landslides. Urban areas are poorly related to landslide activity, because most of the urban areas are located in the northern coastal part of the study area where landslides are limited.

The proposed methodology reveals a relatively simple and quick way of determining the association of physical and anthropogenic parameters with landslide activity in a mountainous terrain. Topographical, geological, hydrological, transportation and landslide location data can easily be found

and their analysis and evaluation are simple and rapid. In regional studies, the applied procedure can be used for the localization of sites prone to landslides and for landslide hazard assessment mapping. Therefore, engineers, planners, decision-makers and environmental managers may utilize the proposed methodology in new and existing spatial planning projects. Additionally, it may be used by the local authorities to guide them in the adoption of policies and strategies aiming at landslide hazard mitigation.

**Author Contributions:** H.D.S. and D.R. conceived the research; H.D.S. and G.D.B. designed the research and the data analysis; K.S. prepared and analyzed the data; H.D.S., G.D.B., K.S., E.K. and D.R. completed the field work; H.D.S. and E.K. created the figures; H.D.S., G.D.B. and K.S. wrote the paper.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Land* Editorial Office E-mail: land@mdpi.com www.mdpi.com/journal/land

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18