Assessment of Neotectonic Landscape Deformation in Evia Island, Greece, Using GIS-Based Multi-Criteria Analysis

Abstract

This study deals with the assessment and mapping of neotectonic landscape deformation in the northern part of the Evia Island (Central Greece). Multi-Criteria Decision Analysis (MCDA) utilizing Analytic Hierarchy Process (AHP) and Weighted Linear Combination (WLC) procedures were conducted for the calculation of the Neotectonic Landscape Deformation Index (NLDI). The study is based on the combination of morphotectonic, geomorphological and geological parameters. The GIS-based spatial MCDA led to the classification of the study area into five classes of neotectonic deformation (from very low to very high) and to a neotectonic deformation map. The results were compared with the outputs of a relative tectonic activity classification approach based on quantitative geomorphic analysis at a regional scale, including site-specific field observations. Areas of high and very high deformation are related to the major active faults of Dirfis, Kandili and Gregolimano–Telethrio. Other minor active normal faults of medium to high seismic risk level, affecting the northern and northeastern parts of the island, are also associated with areas of intense landscape neotectonic deformation.

1. Introduction

Geomorphic indices have been widely and successfully used as survey tools for the identification of areas experiencing rapid tectonic deformation [1,2,3,4,5,6]. Although many criteria for the classification of active tectonics in different areas have been proposed [7,8,9,10], there are only a few studies specifically dealing with the classification of landscape deformation. Among the latter is the article of Argyriou et al. [1], which proposed the Neotectonic Landscape Deformation Index (NLDI) and applied it successfully to assess the neotectonic deformation of the landscape in the active region of western Crete.

Landscapes and landforms in actively deforming areas are controlled by the interactions between tectonic movements and consequent surface processes that can lead to local erosion or deposition [11,12]. Neotectonic deformation along with climate and bedrock lithology is the major factor controlling the geomorphological development and evolution of the landscape. The landscape of Greece reflects in many ways the influence of horizontal and vertical crustal movements that have occurred in the geologically recent past and which are ongoing today.

Morphotectonic analyses along with studies demonstrating the active normal fault evolution scenarios (segmentation and linkage) can contribute to our knowledge about neotectonic deformation in regions with poor historical earthquake data [13,14,15]. A way to understand the sensitivity of the landscape is to consider the impact on it by a relatively small tectonically induced change [11]. The organization of the drainage networks and the shape of their corresponding catchments reflect this impact since in tectonically active areas fluvial processes are influenced by the activity of both regional and local faults [16,17]. For the assessment of the relative neotectonic activity of a region, several studies calculate geomorphic indices independently [18,19], while others combine multiple geomorphometric parameters [7,10,20,21]. In addition, there are some examples of GIS-based Multi-Criteria Decision Analysis (MCDA) studies using various morphometric variables as criteria, to detect landscape response to neotectonic activity and to characterize variations in active tectonics across the landscape of a given area [1,22].

The majority of the morphometric and neotectonic analyses combine several GIS techniques, such as spatial analysis tools and geoprocessing tools that include the use of Digital Elevation Model (DEM) mosaics, aerial photographs, satellite images and the Shuttle Radar Topography Mission (SRTM) elevation model [9,23,24]. Additionally, neotectonic activity may be studied and interpreted with the use of other surface methods such as optical remote sensing [25], which utilizes alternative software for image processing (among others, ERDAS, IMAGINE and ENVI) [26] in order to implement further digital techniques (image enhancement, contrast stretching etc.) [27]. Such studies also use new generation techniques—Light Detection and Ranging (LIDAR) [28], Digital Terrain Model (DTM) [29] and Airborne Laser Scanning DTM (ALS DTM) [30]—in order to get geographic data of higher resolution which will be used in spatial analysis and geostatistical modeling. Often the GIS and remote sensing techniques are combined with borehole profiles and field observations in order to obtain as much information as possible on the neotectonic activity of a specific area [29]. In other cases, the aerial photographs, DEMs and satellite images are visualized in Google Earth and interpreted with the use of field surveys [31]. Alternative resource data (e.g., QuickBird) and technologies (e.g., Permanent Scatterers Synthetic Aperture Radar Interferometry (Ps-InSAR)) along with 3D visualization methods can also be helpful in neotectonic interpretation [32,33]. Apart from the interpretation of the 2D simple DEM visualization and the remote sensing images, there are several different techniques for the 3D visualization of the landscape, including the relief-shading of multiple illumination directions that is widely used for further investigation in neotectonic studies [34] and maximizes the visualization of a surface for a better graphical display and analysis [35]—the combined viewing, surface derivatives, 3D viewing, spatial enhancements [36] and photorealistic and nonphotorealistic techniques [37].

GIS techniques including DEM analyses provide an opportunity to quantify the influence of the processes controlling landscape evolution [38] and have been extensively used in various deformation assessment methodological approaches [39,40,41,42,43]. GIS-based MCDA methods have become particularly attractive for spatial planning and management issues and can be considered as a significant tool for decision-makers, giving solutions to problems that require multiple criteria analyses. Such approaches are also quite common in other geoscience fields, such as landslide research (e.g., [44]) and flood hazard studies (e.g., [45,46]).

The primary objective of this paper is to assess and map the landscape neotectonic deformation of the northern part of the Evia Island (Central Greece), by applying a GIS-based MCDA approach which includes superimposing various thematic maps and implementing a detailed Analytic Hierarchy Process (AHP) followed by application of the Weighted Linear Combination (WLC) method. The study area is considered as a tectonically active area and there is an obvious influence of tectonic processes on the landscape evolution. However, only a few morphotectonic or neotectonic studies have been conducted in the study area and the fault zones are mainly offshore, which denotes that it is difficult to identify fault zone history and evolution utilizing other methods (e.g., paleoseismology). This means that the method used is a valuable tool, not only for the classification of the deformation of the study area, but also for the fault zone segmentation and consequently for the seismic hazard determination.

The study was based on the analysis and combination of nine morphotectonic, geomorphometric and geological factors, according to Argyriou et al. [1]. The parameters of the first two categories are useful in the identification of tectonically active areas and in the determination of the variation of the relative tectonic activity, while the geological factors are important in the sense that they affect the other parameters [1]. The wide and extensive use of these parameters in morphotectonic analyses (e.g., [27,47,48,49,50]) proves that their selection for the assessment of the landscape deformation is correct and rational. The most commonly used morphotectonic parameters are the Amplitude of relief (Ar), which is defined as the maximum difference in elevation within a 1 km² unit area [51], and the normalized Stream Length gradient index (SLk), which is linked to stream power—the ability of a stream to erode its bed and transport sediment from the highest part of its drainage basin to the lower point of the basin [14]. The geomorphological criteria include Drainage density (Dd), Drainage frequency (Fu), Slope gradient (S) and Topographic wetness index (Twi) and are indirectly related to the identification of morphotectonic features of the landscape. The Dd describes the spacing and the distribution of the drainage network in a catchment [52] and also provides information regarding the surface runoff potential, the ground surface steepness, the degree of landscape dissection, rock permeability and resistance to erosion. High Fu values are strongly correlated with high slope and stream channel steepness—indicative of high surface runoff [53]. The S is the rate of change of the Earth surfaces’ slope and is equal to the ratio of the vertical change of the altitude with respect to the horizontal distance, defined by a tangent level to a topographic surface [54,55]. The Twi, also known as compound topographic index (cti), was introduced by Beven and Kirkby [56] and is widely used to quantify topographic control on hydrological processes [57]. The geological criteria are the Lineament density (Ld), which is equal to the total length of tectonic lineaments that fall within each cell of the generated grid [58], the Lineament frequency (Lf), which corresponds to the total number of tectonic lineaments that fall within each cell of the generated grid [58] and the Lithology (Lth), which is used because it affects the spatial distribution of the other parameters.

The Neotectonic Landscape Deformation Index was computed as the sum of the weighted parameters, and a map of its spatial distribution across the study area was produced.

2. Study Area

The study area is the northern part of the Evia Island (Figure 1). Evia, having an area of 3654.53 km², is the second-largest Greek island. It has an elongated shape with a NW–SE direction and is situated between the Aegean Sea to the east and the neotectonic trench of the north Gulf of Evia to the west. To the north, the Oreos–Trikeri straits connects the north Gulf of Evia with the central-western Aegean Sea and separates the island from the mainland Greece. The topography of the study area can be described as a continuous succession of mountains and plains (Figure 1). The main features of the landscape are the mountain ridges of Dirfis (1736.8 m), Kandili (1236 m), Xiron (990 m), Telethrio (969 m) and Lichas (738 m), and the extensive low-lying Neogene basins of Istiea–Oreos and Paliouras.

North Evia consists of Alpine and Post-Alpine geological formations. The Alpine formations belong to the Sub-Pelagonian geotectonic zone and include carbonate rocks of Triassic–Jurassic age, the Upper Jurassic–Lower Cretaceous ophiolite complex and the Upper Jurassic–Lower Cretaceous schist–chert formation. During the period of Upper Jurassic–Lower Cretaceous, the zone emerged, Fe–Ni deposits were formed and subsequently were covered by the Middle–Upper Cretaceous tectonic nappe that consists of Upper Cretaceous (transgressive) limestones and ultrabasic rocks [59]. The Post-Alpine formations consist of marls and sandstones of Upper Miocene and Pleistocene age. Recent Holocene alluvial deposits occur along the major channels of the drainage networks and at the mouth of most of the rivers (Figure 2).

Neotectonic and stratigraphic studies of the Neogene and Quaternary formations of the study area [60,61] indicate the existence of traces of two Miocene–Pliocene extensional tectonic phases (NE–SW and NNW–SSE) that are responsible for the activity of the normal faults and a third Pliocene–Pleistocene phase associated with strike-slip faults (N100°–120°) with a prevailing counterclockwise movement. According to Chousianitis et al. [62], the N–S extension rates across the central part of the Evia Island, accommodated by the coastal fault system, equals 53 nstrain per year (ns/yr) and the extension drops about 50% farther southeast onshore of the Evia Island. The north Gulf of Evia together with the Gulf of Corinth are the most active structural features of Greece where faulting is localized today [63]. The fault system of the study area corresponds to two offshore NW–SE striking antithetic normal fault zones, which bound the mountainous relief of the northern part of the Evia Island and separate the land from the north Gulf of Evia and the Aegean Sea. These are the SW dipping Kandili–Telethrio fault zone in the west, and the NE dipping Dirfis zone in the east. Roberts et al. [64] suggest that the Dirfis fault (that runs along the northern side of Dirfis Mountain) is more active and is characterized by larger displacement. This fault zone continues offshore to the north, following the 200 m isobath, showing evidence of lower activity. According to Palyvos et al. [65], the Prokopi–Pelion fault zone at the SE flank of the Nileas depression acts as an obstacle to the propagation of the Dirfis fault zone. The offshore active normal faults along the west and east coasts of the north part of the island seem to control the evolution of the landscape. The characteristics of the drainage networks of the north Evia Island and their corresponding catchments [10,66] along with landforms such as alluvial fans, triangular facets, fluvial terraces [67] and uplifted marine notches [68,69] suggest active neotectonic uplift. The main faults of the study area, proposed by various authors [63,64,65,70] are depicted in Figure 2.

The north Gulf of Evia seems to be almost free of strong seismic events for the period after 1900 and seismic energy is released in the form of earthquake swarms [71,72], except for three relatively strong events: two earthquakes with local magnitude (ML) 5.2 happened on 17 November 2014, 25.6 and 26.2 km northwest of Chalkis, and a 5.3 magnitude earthquake occurred on 9 June 2015, 25.2 km northwest of Chalkis. However, several historic events have been recorded [73].

3. Materials and Methods

3.1. Geodatabase

For the purpose of this study, an ArcGIS spatial geodatabase was designed and implemented for the northern part of the Evia Island. The primary GIS layers include coastline, contour lines at 20 m intervals, elevation points, drainage networks, watersheds, lithology and faults. This dataset is derived from analogue data sources including topographic maps at 1:50,000 scale, published by the Hellenic Military Geographical Service (HMGS) and geological maps at 1:50,000 scale, published by the Institute of Geology and Mineral Exploration of Greece (IGME). These maps were scanned and georeferenced to a common coordinate system, the Greek Geodetic Reference System 1987 (GGRS’87) (during the conversion from Hatt to GGRS’87 an error of ±0.20 m occurred). The use of ArcGIS software and several extensions (e.g., Spatial Analyst) contributed to the extraction of several secondary thematic layers (e.g., DEM with 50 × 50 m cell size, hillshade map etc.) from the primary ones. The main secondary datasets (the thematic layers of the selected factors/criteria used for the assessment of the neotectonic landscape deformation and the final NLDI itself) were extracted from the DEM.

Published papers and fault databases [60,61,64,65,67,75,76,77,78] along with DEM analysis and Google Earth image interpretation were used to recognize and map the faults of the study area. Field control was necessary to correct and validate the preliminary main fault zone map that was initially drawn.

3.2. Multi-Criteria Decision Analysis (MCDA)

Using a GIS process entirely, spatial MCDA combines and converts the input geographical data (criteria map layers) into a spatial decision tool (final map) in order to support decision makers [79]. The stages of implementation of the MCDA carried out are similar to the general steps of other spatial methodologies using multiple criteria (e.g., [80,81]). Once the problem is determined (which in our case is the assessment of the neotectonic landscape deformation), the variables (criteria of the analysis) are selected and are reclassified, homogenized and quantified, and hence the decision table can be filled. Then the weighting coefficients are determined (weighting configuration), the weighted variables are combined, the consistency test is performed and finally the concluding proposals are presented and the validation is performed by comparing the results with the real conditions [82]. The weights, rank values and classes of the variables/criteria were assigned using AHP. Then using the WLC, their classes are multiplied with the corresponding weights and their values are summed up to produce the NLDI map. The methodological framework followed for the assessment of the neotectonic deformation of the north Evia landscape is shown in Figure 3.

Criteria Used

To assess the neotectonic landscape deformation of the northern part of the Evia Island, we used the NLDI proposed by Argyriou et al. [1], the calculation of which requires the analysis of nine morphotectonic, geomorphological and geological criteria (Figure 3).

Regarding the parameter Ar, many ways have been suggested for its calculation [83,84,85]. In this study it is calculated by subtracting DEMmin (the minimum of the input cells that are encompassed by the extent of the cell) from DEMmax (the maximum of the input cells that are encompassed by the extent of the cell) for each cell [1]. In the spatial analysis the centroid point of each cell is used. In order to obtain the suitable grid and the final spatial distribution map of Ar, the spatial interpolation method of kriging was used [47]. Ar provides information on regional uplift and on the recognition of faults and their kinematics, whereas its spatial distribution can lead to the identification of regions with differential uplift or to recent vertical displacements of uplifted or subsided blocks [47,86,87]. This parameter is also used in the assessment of the terrain morphological characteristics and of the degree of landscape dissection [88].

The Stream Length gradient index (SL) is given by the formula [89]:

$$\mathrm{SL}=(\mathsf{\Delta }\mathsf{H}/\mathsf{\Delta }\mathrm{L})\ast \mathrm{L}$$

(1)

where ΔH = change in elevation of the stream reach, ΔL = length of the reach, which is the horizontal distance corresponding in each case to ΔH, and L = the total channel length from the midpoint of the reach of interest upstream to the highest point of the channel [14].

Since the SL index is a sensitive variable to changes in channel slope, its spatial distribution represents the gradient changes along the stream channels of a drainage network. For this reason, SL index allows the evaluation of the relationships among possible tectonic activity, rock resistance and topography [14]. Abrupt changes in the gradient of a river can be associated with active tectonics. Troiani et al. [85] used this index in small catchments in central Italy and they found it to be a valid tool to detect the response of the topography to regional uplift processes, even if this response cannot be detected by other morphotectonic parameters. Additionally, SL index has been successfully tested for its usefulness in the detection of tectonic displacements in several tectonically active regions, e.g., [2,14].

Since SL index is strongly influenced by the length of the channel, when SL values of rivers of different channel length are compared, a normalization factor—that is, the graded river gradient (k)—must be used in order to avoid partiality [90]. The graded river gradient (k) outlines the steeper stream segment, and normalizes the effect of stream power and stream order. In this way the SLk index allows the comparison of stream gradient among rivers of different characteristics [2]. The SLk index is calculated by the following formula [91]:

$$\mathrm{SLk}=\mathrm{SL}/\mathrm{k}$$

(2)

where SL = Stream Length gradient index (Equation (1)) and k = graded river gradient, given by the formula:

$${\mathrm{k}=(\mathrm{C}-\mathrm{h}}_{\mathrm{f}})/\mathrm{lnLt}$$

(3)

where C = elevation of the river head, ${\mathrm{h}}_{\mathrm{f}}$ = elevation of the river mouth and $\mathrm{Lt}$ = total length of the river.

We followed the methodology proposed by Pérez-Peña et al. [90], which uses a constant spacing (150 m) for the main stream segments instead of the contour line equidistance that was proposed by Keller and Pinter [14].

From the SLk values an anomaly map was created using the ordinary kriging method (the spherical semivariogram model and a search radius of 15 points were used).

The parameter Dd is given by the formula [92]:

$$\mathrm{Dd}=\mathsf{\Sigma }\mathrm{L}/\mathrm{A}$$

(4)

where ΣL = total length of all the ordered streams and A = area of the drainage basin.

High values of the parameter Dd imply a dissected terrain, while low values indicate terrain with long hill slopes. Many studies suggest that this parameter reveals the effect of active tectonics [93,94] since low values indicate the prevention of a properly developed pattern that is often related to tectonic processes [1,95]. The spatial distribution of Dd was made using the line density tool of the ArcGIS toolset that calculates the density of linear features in units of length per unit of area.

The parameter Fu is given by the formula [92]:

$$\mathrm{Fu}=\mathrm{N}/\mathrm{A}$$

(5)

where N = total number of channels of all stream orders and A = area of the drainage basin.

The permeability of the geological formations as well plays a crucial role in the surface runoff. Hence higher values of the parameter Fu are related to impermeable formations [96]. However, in many cases the opposite is observed due to the effect of the tectonic activity on the development of the drainage network [97]. Uplifted regions can prevent proper development of drainage pattern and can be related to a lower number of stream channels [1]. For the construction of the Fu spatial distribution map of the study area, the ArcGIS Spatial Analyst extension (Kernel density toolset) was used.

The parameter S is estimated using the formula:

$$\mathrm{S}=\sqrt{{(\mathrm{G}}^2{+\mathrm{H}}^2)}$$

(6)

where G = east-to-west gradient and H = north-to-south gradient.

S values are given in degrees ranging from 0° to 90° or as percentage (100% corresponds to 45°). A slope map of an area highlights the effects of the tectonic activity on the morphology of the landscape [98]. S was calculated from the DEM with the slope tool of the ArcGIS surface toolset.

There are a number of methods for the calculation of the parameter Twi from a DEM. We calculated Twi by the following equation [99]:

$$\mathrm{Twi}=\mathrm{Ln}\left(\right(“\mathrm{facc}\_\mathrm{dem}”+0.001)/((“\mathrm{slope}\_\mathrm{dem}”/100)+0.001\left)\right)$$

(7)

where facc_dem = flow accumulation and slope_dem = slope gradient.

Its calculation leads to the determination of the spatial distribution of soil moisture due to the topographic changes or to the identification of places where the water tends to accumulate due to the topography, since the latter controls the spatial variation of hydrological conditions [57]. Twi is considered here to provide information on the characterization and identification of areas that might be influenced by structural control and not to explain hydrological processes and soil moisture distribution. Higher positive values indicate increased humidity, high moisture accumulation and characterize gently sloping surfaces, such as the foot of the slopes or areas composed of alluvial deposits where water tends to accumulate [100,101]. On the contrary, low negative Twi values indicate drier conditions and are detected near water divides, at the areas of increased runoff or where the slope gradient is high (e.g., incised V-shaped valleys) [1,102]. Hence, Twi is helpful for detecting possible fault-controlled topographic features [101].

To calculate the parameter Ld and to create its spatial distribution map, we used the line density tool of the ArcGIS toolset that calculates the density of linear features in units of length per unit of area. The lineaments used were the faults as shown in the geological map of the study area (Figure 2). These were the most appropriate data among the lineaments that were automatically extracted by the combination of two shaded relief images (at different sets of azimuth angle from 0° to 360° and a sun elevation of 45°) created from the DEM [103,104]. Areas of higher Ld values reflect high degree of shearing [105], highly permeable rocks [106] and deformation of greater intensity [107].

As in the case of Ld, for the calculation of parameter Lf we used the lineaments shown in the geological map of the study area (Figure 2). Its spatial distribution map was produced with the use of ArcGIS Spatial Analyst (Kernel density toolset). High values of Lf reveal regions undergoing high brittle rock fracturing, fault propagation and intense tectonic activity [108].

Lithology (Lth) is a geological criterion essential in the MCDA procedure since it discriminates between geological formations in terms of their response to the passage of seismic waves. The most “weak”, unconsolidated geological formations are characterized by reduced seismic wave velocities and thus by increased wave amplitude and seismic risk in comparison to the more resistant “strong” geological formations that do not produce this amplification indicating lower seismic risk. The geological formations of the study area were classified into four classes according to their lithology considering the Geological Strength Index (GSI) [109]. The first class (group 1) includes alluvial deposits, recent talus cones and scree, and semiconsolidated talus cones. The second class (group 2) consists of weak geological formations including marls, sandstones, schist–cherts and travertines. The third class (group 3) is characterized by medium resistance and consists of schists and phyllites, while the fourth class (group 4) includes hard rocks such as limestones, dolomites, greywackes, consolidated conglomerates, ophiolites, peridotites and quartzites (Figure 4i). The source for the geological formations of the study area was the geotechnical map of Greece at 1:500,000 scale (IGME [110]).

All of these parameters/criteria are related in a way to tectonic deformation and each is represented by a spatial distribution map of the classified values, reclassified to a regular scale (up to 10 classes) and then by a thematic map of the spatial distribution of the normalized values (with a single reference unit equal to 100). The analysis and interpretation of the thematic maps referring to each criterion can reveal zones of high and very high neotectonic activity since areas characterized by a specific range of classified values of each individual parameter can be linked to indications of neotectonic landscape deformation [1]. For most of the criteria used, including Ar, SLk, S, Ld and Lf, higher classified values correspond to more intense landscape deformation. However, for the parameters Dd, Fu and Twi, higher values characterize areas of lower landscape deformation. Regarding Lth, areas composed of weak rocks (group 1) are considered as prone to more intense relative tectonic deformation [1].

3.3. Analytic Hierarchy Process (AHP)

The AHP method [111] is flexible and well-structured, and simple to use. Its role is to deconstruct a problem, in a hierarchical context, at a level where the data are compared in pairs in order to assess the weight of each, at the next level. The method of comparisons per pair is the most commonly used process for the calculation of criteria weight coefficients in MCDA applications [112]. AHP is based on the allocation of weights to the criteria used, according to their importance. To determine the relative importance of each criterion (which according to Saaty [113] means how much more one element dominates another with respect to a given attribute), pairwise comparisons take place using the scale of absolute number of Saaty and Vargas [114]. This scale offers value judgments both quantitatively and qualitatively. Completion of the degree of significance between two variables requires field experience, and knowledge of the subject and/or the opinion of the experts [115]. In this way, a preference table is created which is the first step (a) of the AHP (

Table 1. Pairwise comparison matrix of the criteria concerning the goal—step (a) (αii: preference values).

a11	a12	…	a1n
a21	a22	…	a2n
	…
an1	an2	…	ann
Sum_1	Sum_2	…	Sum_n

).

After assigning the preference values (αii) (step a), the columns (Sum_i) (

Table 2. Table of values according to steps: (b) (preference values are divided by the column sum), (c) (sum of matrix rows) and calculation of the weights at the final step (d) (row sum divided by n).

					Weighting coefficient (wi)
a11/Sum_1	a12/Sum_2	…	a1n/Sum_n	Row_sum1	Row_sum1/n
a21/Sum_1	a22/Sum_2	…	a2n//Sum_n	Row_sum2	Row_sum2/n
	…				…
an1/Sum_1	an2/Sum_2	…	ann/Sum_n	Row_sumn	Row_sumn/n

) are summed and then (steps b and c) the preference values of each column (αii) in the table are divided by the sum of the corresponding column (Sum_i) and the values are summed again in rows this time (Row_sumi) (

Table 3. Calculation of the weighted sum vector (WVm) (where wi: weighting coefficient (Table 2) and αij: preference values (Table 1)).

WV1 = a11*w1 + a12*w2 + … + a1n*wn

WV2 = a21*w1 + a22*w2 + … + a2n*wn

…

WVm = an1*w1 + an2*w2 + … + ann*wn

). Finally (step d), the set of each line (Row_sumi) is divided by the number of variables (n) and the result equals the weighting coefficient of each criterion.

During the comparisons per pair, a subjective assessment of the importance of one criterion over the other is carried out; hence, a reduced number of criteria is more appropriate for the application of the method [116]. However, in the final stage the judgments are consistent as long as the Consistency Ratio (CR) (given by Equation (8)) is less than 0.1, tolerated if 0.1 < CR < 0.2 and rejected if CR > 0.2.

CR is calculated by the following formula:

$$\mathrm{CR}=\mathrm{CI}/\mathrm{RI}$$

(8)

where RI = a random index of specific values according to Saaty [112] and CI = the Consistency Index computed by the following equation:

$$\mathrm{CI}=(\mathsf{\lambda }\max -\mathrm{n})/(\mathrm{n}-1)$$

(9)

where n = number of items compared in the matrix and λmax = the maximum value of the eigenvalue that is obtained by the equation:

$$\mathsf{\lambda }\max =(1/\mathrm{n})\ast {\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}(\mathrm{WVi}/\mathrm{wi})$$

(10)

where ${{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}(\mathrm{WVi}/\mathrm{wi})=\mathrm{vector}\mathrm{coherence}\left(\mathrm{C}\right)$ and wi = the weighting coefficient estimated according to Table 2 and WVi = weighted sum vector (calculated according to Table 3).

3.4. Weighted Linear Combination (WLC)

In the final phase of the methodology, an overall priority score is created using the WLC technique and the overlay capabilities of the GIS software, in order to determine the composite map [116]. The data used were in raster format with cell size 50 × 50 m and a common projection (GGRS’87). Each of the raster layers represents a criterion. The values of each criterion were classified and standardized based to a common numeric range (100) and a final NLDI map came up after the aggregation of the weighted values of their cells [79], according to the following equation proposed by Argyriou et al. [1]:

$$\mathrm{NLDI}={\displaystyle \sum }\mathrm{Fw}\ast \mathrm{Fr}$$

(11)

where Fw = weight of each factor and Fr = standardized rank value of each class of the criteria.

The resulting map has the same range of values as the standardized criteria maps that were used [79], and the higher values of the index correspond to highly deformed areas or to areas more susceptible to neotectonic landscape deformation. The NLDI values of the study area were then classified into five (5) classes using the natural breaks classification method and the final neotectonic landscape deformation map was created. Class 1 of the lower values corresponds to very low deformation, class 2 corresponds to low, class 3 to moderate, class 4 to high and finally class 5 of the higher values corresponds to areas of very intense neotectonic landscape deformation. The natural breaks classification method was chosen since it sets the boundaries to data values at relatively large jumps [117].

4. Results

4.1. Conditioning Factors

Ar values in the study area range from 13.3 to 798 m. Values higher than 300 m are concentrated at the footwalls of the main fault zones of Dirfis, Kandili and Gregolimano–Telethrio, while lower values correspond to the low-lying alluvial plains, as shown on the Ar spatial distribution map (Figure 4a).

SLk index was calculated at 4311 points and the index values range between 0 and 127.72 m. An SLk anomaly map was produced from these points using the kriging method (Figure 4b). The values ranging from 0 to 24.99 were grouped in 10 classes using the natural break classification method. The points with the highest values that represent a trend away from other classes (anomalies) are recorded at the uplifting blocks of the Telethrio and Dirfis fault zones, while similar anomalies can also be found at the Kandili and Oreos straits fault zones without a specific directional trend. The influence of the tectonic features of the study area is reflected on the SLk index map since the values significantly increase at the major fault-controlled areas.

Dd and Fu values range from 0 to 5.23 and 0.21 to 4.52, respectively. The spatial distributions of these parameters are presented in Figure 4c,d. As expected, the drainage texture of the area is related to the spatial distribution of the various lithological types and to the topography. The highest values are concentrated at the southeastern part of the area where the terrain is dissected and the rough and steep relief of the Dirfis footwall increases the surface runoff.

Regarding S, the highest value is 67.21° and the spatial distribution of this criterion is depicted on the map of Figure 4e. The major part of the study area is characterized by low to medium slopes (0°–10° and 10°–20° respectively) and only 0.4% of the total area has slopes greater than 45°. Lower slopes are located at the north, east and southwest part of the study area where the Neogene drainage basins occur, while the highest slopes are found at the southeast, southwest and northwest part of the study area where the mountainous landscapes of Dirfis, Kandili and Telethrio Mts. exist. The significantly steeper morphology of these areas suggests intense deformation due to neotectonic activity.

Twi determines the spatial distribution of soil and surface saturation with regard to the influence of topography. The values of the index range between −7.62 and 12.96 in the study area. Since this index is related to slopes, as water tends to accumulate at the foot of slopes, high values are concentrated in areas prone to water saturation, such as the alluvial plains or along the gently sloping main channels of the Neogene basins in the southern, central and northern parts of the study area (Figure 4f). On the contrary, low values correspond to deeply incised V-shaped valleys, high relief surfaces with low moisture accumulation, and longitudinal ridges, which are indicative of an intensely deformed, possibly fault-controlled landscape, and characterize the mountainous parts of the study area.

Ld values vary from 0 to 12.48 and the spatial distribution map of this index is presented in Figure 4g. The higher values are in line with the presence of the main faults of the study area, whereas high Ld values also characterize the northwestern part of the island where several minor faults prevail.

Lf values range from 0 to 1.73 with the higher values concentrated at the southern termination of the Dirfis fault zone, at the Prokopi–Pelion fault zone and at the northern part of the study area where minor faults are common (Figure 4h).

Lth is based on the geological formations and surficial deposits cropping out in the study area. Figure 2 shows the spatial distribution of Lth after the classification of the geological formation and the surficial deposits into four groups according to their response to the passage of seismic waves. The southern part of the study area is dominated by hard lithologies (group 4), whereas the prevailing geological formations of the northern part are those of group 2 (weak rocks). Although Lth is not a direct indication of neotectonic deformation, it is integrated in the MCDA procedure since it plays an important role affecting the other factors.

4.2. Multi-Criteria Decision Analysis

The spatial distribution maps of the nine criteria shown in Figure 4a–i were the thematic maps used as inputs in the MCDA. All values of all the involved parameters are classified and normalized to 100. Higher classified ranking values correspond to areas of intense neotectonic deformation or to areas more susceptible to deformation, while lower values correspond to less neotectonically deformed areas.

The results of the AHP procedure, the extraction of criteria weights and the calculation of the CR (steps a, b–c, d of the methodology described) are shown in

Table 4. Pairwise comparisons, calculation of weighting coefficients and CR. (Ar: Amplitude relief, SLk: normalized Stream Length gradient index, Dd: Drainage density, Fu: Drainage frequency, S: Slope gradient, Twi: Topographic wetness index, Ld: Lineament density, Lf: Lineament frequency, Lth: Lithology, W: weights, WV: vector of weighted sum, C: coherence vector, λ: maximum value of eigenvalue, CI: consistency index, CR: consistency ratio).

Step a
		Ar	SLk	Ld	Lf	Dd	Fu		S	Twi	Lth
	Ar	1	1	2	2	3	3		3	6	5
	SLk	1	1	2	2	3	3		3	6	5
	Ld	0.5	0.50	1	1	3	3		3	5	3
	Lf	0.5	0.50	1	1	3	3		3	5	3
	Dd	0.33	0.3	0.33	0.33	1	1		1	3	4
	Fu	0.33	0.3	0.33	0.33	1	1		1	3	4
	S	0.33	0.33	0.33	0.33	1	1		1	2	4
	Twi	0.17	0.17	0.20	0.20	0.3	0.33		0.50	1	3
	Lth	0.20	0.20	0.33	0.33	0.25	0.25		0.25	0.33	1
	Sum	4.37	4.37	7.53	7.53	15.58	15.58		15.75	31.33	32.00
Step b–c
		Ar	SLk	Ld	Lf	Dd	Fu		S	Twi	Lth	Sum
	Ar	0.23	0.23	0.27	0.27	0.19	0.19		0.19	0.19	0.16	1.91
	SLk	0.23	0.23	0.27	0.27	0.19	0.19		0.19	0.19	0.16	1.91
	Ld	0.11	0.11	0.13	0.13	0.19	0.19		0.19	0.16	0.09	1.32
	Lf	0.11	0.11	0.13	0.13	0.19	0.19		0.19	0.16	0.09	1.32
	Dd	0.08	0.08	0.04	0.04	0.06	0.06		0.06	0.10	0.13	0.65
	Fu	0.08	0.08	0.04	0.04	0.06	0.06		0.06	0.10	0.13	0.65
	S	0.08	0.08	0.04	0.04	0.06	0.06		0.06	0.06	0.13	0.62
	Twi	0.04	0.04	0.03	0.03	0.02	0.02		0.03	0.03	0.09	0.33
	Lth	0.05	0.05	0.04	0.04	0.02	0.02		0.02	0.01	0.03	0.27
	Sum	1.00	1.00	1.00	1.00	1.00	1.00		1.00	1.00	1.00
Step d	CR Calculation
		W			W	WV			C	λ	CI	CR
	Ar	0.21		Ar	0.21	2.026			9.5350	9.4642	0.0580	0.0400
	SLk	0.21		SLk	0.21	2.026			9.5350
	Ld	0.15		Ld	0.15	1.423			9.6763
	Lf	0.15		Lf	0.15	1.423			9.6763
	Dd	0.07		Dd	0.07	0.684			9.4150
	Fu	0.07		Fu	0.07	0.684			9.4150
	S	0.07		S	0.07	0.647			9.3681
	Twi	0.04		Twi	0.04	0.339			9.2619
	Lth	0.03		Lth	0.03	0.279			9.2954
	Sum	1.00						Sum	85.178

. The procedure begins with the pairwise comparisons of the nine factors and the calculation of the weighting coefficients. All criteria might be controlled by lithology and thus take the values 5 or 4 in the comparison matrix (except for Ld, Lf and Twi, which take the value 3). In the same way, the factors Ar, SLk, Ld and Lf are of strong importance in comparison to the Twi index, while the factors Dd, Fu and S are affected by the Twi index and these criteria are considered as slightly to moderately more important in the decision matrix. Ar, SLk, Ld and Lf are slightly to moderately more important than geomorphological parameters (Dd, Fu, S). Ar and SLk are slightly more important in comparison to Ld and Lf. Equal importance, providing comparable information, is assigned to the groups of parameters: Dd–Fu–S, Ld–Lf, and Ar–SLk. The weighting coefficients, determined with the application of the comparisons per pair of the AHP in decreasing order follow: for SLk and for Ar the weight is 0.21, for Ld and Lf the weight is 0.15, for Dd, Fu and S it is 0.07, for the Twi the weight is 0.04 and the smallest weight equal to 0.03, is calculated for Lth.

The calculated CR (Table 4) in this study is 0.0400 (lower than 0.1) within the range for the acceptance of the consistency of judgments in the pairwise comparison matrix.

As described in the methodology, the neotectonic landscape deformation was calculated following the WLC procedure, using the classified and normalized values of the nine criteria involved and their weighting coefficients (

Table 5. Table of the classified (using the natural breaks method) and standardized (100) values of each criterion used in the analysis (Ar: Amplitude relief, SLk: normalized Stream Length gradient index, Dd: Drainage density, Fu: Drainage frequency, S: Slope gradient, Twi: Topographic wetness index, Ld: Lineament density, Lf: Lineament frequency, Lth: Lithology, F: Factors, W: Weights, C: Classes, R.V.: Rank Values, S.R.: Standardized Rating).

F.	W.	C.			R.V.	S.R.	F.	W.	C.			R.V.	S.R.
(Ar)	0.21	500.1	–	798	10	100	(S)	0.07	42.45	–	67.21	10	100
		450.1	–	500	9	90			35.33	–	42.44	9	90
		400.1	–	450	8	80			30.59	–	35.32	8	80
		350.1	–	400	7	70			26.10	–	30.58	7	70
		300.1	–	350	6	60			21.89	–	26.09	6	60
		250.1	–	300	5	50			17.67	–	21.88	5	50
		200.1	–	250	4	40			13.45	–	17.66	4	40
		150.1	–	200	3	30			8.97	–	13.44	3	30
		100.1	–	150	2	20			4.23	–	8.96	2	20
		13.3	–	100	1	10			0	–	4.22	1	10
(SLk)	0.21	22.5	–	24.99	10	100	(Twi)	0.04	(−7.62)	–	(−4.98)	10	100
		20	–	22.49	9	90			(−4.97)	–	(−1.64)	9	90
		17.5	–	19.99	8	80			(−1.63)	–	1.4	8	80
		15	–	17.49	7	70			1.41	–	2.62	7	70
		12.5	–	14.99	6	60			2.63	–	4.04	6	60
		10.01	–	12.49	5	50			4.05	–	5.76	5	50
		7.51	–	10	4	40			5.77	–	7.69	4	40
		5.01	–	7.5	3	30			7.70	–	9.92	3	30
		2.51	–	5	2	20			9.93	–	12.96	2	20
		0	–	2.5	1	10			12.97	–	18.23	1	10
(Dd)	0.07	0	–	0.53	10	100	(Ld)	0.15	11.24	–	12.48	10	100
		0.54	–	1.05	9	90			9.99	–	11.23	9	90
		1.06	–	1.57	8	80			8.74	–	9.98	8	80
		1.58	–	2.09	7	70			7.5	–	8.73	7	70
		2.1	–	2.62	6	60			6.24	–	7.49	6	60
		2.63	–	3.14	5	50			5	–	6.23	5	50
		3.15	–	3.66	4	40			3.71	–	4.99	4	40
		3.67	–	4.18	3	30			2.5	–	3.7	3	30
		4.19	–	4.71	2	20			1.21	–	2.49	2	20
		4.72	–	5.23	1	10			0	–	1.2	1	10
(Fu)	0.07	0.21	–	0.93	10	100	(Lf)	0.15	1.56	–	1.73	10	100
		0.94	–	1.34	9	90			1.39	–	1.55	9	90
		1.35	–	1.68	8	80			1.22	–	1.38	8	80
		1.69	–	1.97	7	70			1.05	–	1.21	7	70
		1.98	–	2.24	6	60			0.87	–	1.04	6	60
		2.25	–	2.54	5	50			0.40	–	0.86	5	50
		2.55	–	2.90	4	40			0.53	–	0.39	4	40
		2.91	–	3.34	3	30			0.36	–	0.52	3	30
		3.35	–	3.86	2	20			0.18	–	0.35	2	20
		3.87	–	4.52	1	10			0	–	0.17	1	10
							(Lth)	0.03	group 1			4	100
									group 2			3	75
									group 3			2	50
									group 4			1	25

).

After the reclassification of the NLDI values (Figure 5) into five categories with the Jenks natural breaks classification method, the final map was produced (Figure 6) showing the northern part of Evia Island classified into five zones of different degree of deformation (from very low to very high):

• Very low (17.15–30.32)
• Low (30.33–37.15)
• Moderate (37.16–43.74)
• High (43.75–53.25)
• Very high (53.26–79.35)

On the final maps (Figure 5 and Figure 6), one can notice that—apart from the areas of high to very high landscape neotectonic deformation that are directly affected by the Dirfis, Kandili, Gregolimano–Telethrio and Oreos straits active major fault zones—high values are also observed at the northeastern part of the study area. This area is characterized by several local minor normal faults of moderate to high seismic risk [75], which are probably responsible for the high NLDI values and for the rough landscape.

5. Discussion

The main goal of this study was to assess the neotectonic landscape deformation of north Evia Island through the estimation and mapping of the NLDI, an index proposed by Argyriou et al. [1], through a GIS-based MCDA approach using the AHP and WLC. The analysis led to the identification of both areas of high and very high neotectonic deformation characterized by a rugged relief with steep slopes that have been affected by active faults, and areas of low tectonic activity characterized by a smooth relief. The areas identified as of high and very high deformation are related to the activity of the main normal fault zones of Dirfis, Kandili and Gregolimano–Telethrio, or are affected by minor local active faults located at the northern and northeastern parts of the island.

The verification of the spatial distribution of the neotectonic deformation zones included field observations and the comparison of the results with the findings of other relative tectonic activity classification approaches performed in the study area (Figure 7, Figure 8 and Figure 9).

The final NLDI map was overlaid by the Index of Relative Tectonic Activity (IRTA) map [10]. IRTA takes into account five morphometric parameters, including the Hypsometric integral (Hi), the Asymmetry factor (Af), the drainage basin Slope (Sb), the Relief ratio (Rh) and the Melton’s ruggedness number (M), and are calculated for 189 drainage basins in an attempt to assess the relative tectonic activity within the study area (Figure 7). Additionally, in an attempt to better evaluate the relationships between the results of the two indices, five topographic sections crossing different parts of the island were constructed approximately perpendicular to the main fault zones (Figure 7 and Figure 8). Along each topographic profile, the diagrams of both NLDI and IRTA results are also depicted. In general, the results of the GIS-based MCDA approach are in good agreement with the spatial distribution map of the IRTA. In most cases, the areas of high landscape deformation correspond to catchments of high relative tectonic activity, whereas low relative tectonic activity catchments are related to areas of low neotectonic deformation. Landscapes of the highest values of both NLDI and IRTA indices correspond to regions affected by the main active normal fault zones. The map (Figure 7a) that was produced by the overlay of the NLDI and IRTA maps along with the perpendicular to the fault zones cross sections can lead to the following inferences.

• The results of the neotectonic landscape deformation assessment revealed that the landscape of the most extensive catchments that drain the uplifting block of the Dirfis fault zone (catchments 85 and 101) are characterized by very low values of NLDI, indicating that this offshore fault zone is not continuous but segmented into three discrete portions (the segmentation points/areas are marked with black arrows in Figure 7b). The east segment extends from the eastern termination of the fault to catchment 85 (Figure 9e, location 4), the middle segment lies between the catchments 85 and 101, while the western segment extends from catchment 101 to the western end of the fault zone. This segmentation limits the magnitude of a potential earthquake caused by this fault zone. The Dirfis fault zone terminates to the north at its intersection with the Prokopi–Pelion fault, which has a roughly NE–SW trend and acts as a barrier separating the Dirfis fault zone from its probable extension to the north. North of this area the uplift is limited and the landscape is smoother, indicating low neotectonic deformation. The landscape along the cross-section of the Dirfis fault zone (Figure 8e) and the higher values of the indices at the east and west coasts are indicative of horst morphology in this part of the island.
• The Kandili fault zone also seems to be segmented into two portions (green arrow in Figure 7b). The north segment has a length of 22.73 km, whereas the length of the south segment is 8.64 km. A fault that crosses the Kandili fault zone (Figure 9b, location 1) (near drainage basin 60) seems to be responsible for this segmentation. The scarp of this crossing fault has been verified during the fieldwork (Figure 9c, location 2). On the tectonic map of the north Gulf of Evia published by Sakellariou et al. [70], this fault extends northwestward, crossing the Gulf reaching the Ag. Konstantinos fault zone at the opposite coast. According to Palyvos et al. [65], this fault extends eastward up to Kireas stream, north of the Prokopi–Pelion fault zone. The Kandili fault zone cross-section (Figure 8d) shows higher values of both NLDI and IRTA indices at the east and west coasts, which supports the view of a horst morphology at this part of the island. The same stands for the Lichada Peninsula (Figure 8a).
• The Gregolimano–Telethrio fault zone consists of two segments and both seem to be very active based on the degree of the deformation of its footwall landscape (Figure 9d, location 3). The Gregolimano and Telethrio topographic profiles (Figure 8c,d) indicate that both relative tectonic activity and landscape neotectonic deformation increase significantly from east to west. This means that the Gregolimano–Telethrio fault zone shows higher levels of relative neotectonic activity.
• The landscape of the northeast part of the study area is classified as highly deformed. The area south of the drainage basin 160 is affected by a large number of smaller faults of moderate seismic risk level.

Intense neotectonic activity at the areas of high and very high NLDI values is verified by the presence of uplifted/submerged tidal notches. A single uplifted tidal notch and in some places a second higher older one are preserved at the base of the coastal cliffs of the Dirfis, Telethrio and Kandili fault zones, suggesting recent tectonic uplift. Another submerged tidal notch is found at the southwestern part of the study area and is related to neotectonic subsidence also responsible for the submerged ancient constructions of this area. Both uplift and subsidence have been related to recent tectonic movements [68], confirmed as well by geophysical observations [118]. Other landforms/markers of active tectonics such as V-shaped valleys, scarp/dip topography and stream piracies are also identified at the footwalls of the main fault zones.

The analysis of the paper has some limitations. The criteria weights were established based on pairwise comparison that is among the subjective weighting methods. It would be useful to calculate the weights using mathematical functions, which correspond to the objective weighting methods that do not include the interference of the analyst. Other limitations are related to the validation method. Except from the comparison of the results with the values of IRTA, the calculation of other indices of active tectonics, e.g., Iat (Index of relative active tectonics) [7] or Riat (Relative index of active tectonics) [119] and the use of their values for comparisons with the results, would be valuable during the validation phase. Furthermore, the landforms/markers of active tectonics, used to confirm the spatial distribution of deformation, were mapped at a regional scale. More detailed geomorphological mapping of the areas of high to very high deformation would be useful in further interpretation of the results. Despite these limitations, the results are reliable and the maps of the neotectonic landscape deformation can be used for further research and in seismic hazard analysis.

6. Conclusions

This study focused on the determination of the neotectonic landscape deformation of the northern part of the Evia Island based on the values and the spatial distribution of the Amplitude of relief (Ar), the normalized Stream Length gradient index (SLk), the Drainage density (Dd), the Drainage frequency (Fu), the Slope gradient (S), the Topographic wetness index (Twi), the Lineament density (Ld), the Lineament frequency (Lf) and the Lithology (Lth). The map of the spatial distribution of the Neotectonic Landscape Deformation Index (NLDI) values, which has been obtained by the combination of these nine morphotectonic, geomorphological and geological indices through a GIS-based multi-criteria decision analysis approach using the Analytic Hierarchy process and the Weighted Linear Combination, confirms that the landscape of the study area clearly reflects the effect of recent tectonic activity. The estimated landscape deformation is high to very high in the vicinity of the major active normal faults. More specifically, high values of deformation are identified at the footwall of the Dirfis fault zone (especially at its NW and SE termination), and at the footwalls of the fault zones of Kandili (mainly at its southern end), and Gregolimano–Telethrio. In addition, the landscape of the northeastern part of the study area is highly tectonically deformed as it is controlled by local minor active normal faults.

Apart from mapping of faults, this approach was also useful in mapping the terminations of the fault zones and in detecting the degree of their segmentation. In addition, it allowed us to detect fault segments and boundaries. The analysis revealed that the Kandili fault zone consists of two portions, whereas the Dirfis fault zone seems to consist of three segments. This finding is of crucial importance as the potential maximum earthquake magnitudes expected by these fault zones are significantly lower since each fault segment ruptures independently of the other. The Gregolimano–Telethrio fault zone seems to be more active since its footwall landscape is characterized by very high values of deformation. The extension of the Dirfis and of the Oreos straits fault zones proposed by other studies were not confirmed by the results of this approach. The findings of this study are in good agreement with the results of the application of a composite index of relative tectonic activity named IRTA. The results were also verified by means of field observations.