Investigating a Large-Scale Creeping Landmass Using Airborne and Subsurface Remote Sensing Techniques

1. Introduction

Landslides represent complex geological phenomena composed of various formations, characterized by abrupt or gradational variations in their physical properties. Recognized as some of the most catastrophic natural hazards [1], with notable frequency in Greece. Their complex nature necessitates comprehensive investigation, including several aspects such as the delineation of the geomorphological regime, assessment of the in-situ stability conditions, characterization of the subsurface strata and rock properties, as well as determination of the kinematic parameters [2].

In the case of a fast-moving landslide event, direct monitoring and investigation of the physical and mechanical parameters that control slope failure is difficult. However, landslides characterized by low rates of movement, ranging from millimeters to multiple meters per year, offer an ideal opportunity for detailed investigation of the landslide processes under near-stable conditions [3]. Despite not being immediately disastrous, slow-moving landslides can have lasting harmful effects on local communities by building damaging, degrading infrastructure [4,5] and agricultural lands [6]. Moreover, their systematic movement influences the morphological evolution of the landscape through erosional processes [7].

This type of movement typically occurs in weak materials such as unconsolidated soil deposits with high clay content or highly weathered, fractured and altered rocky formations [7,8]. Additionally, they are situated in areas characterized by high levels of seasonal precipitation [9] with a complex subsurface structure and hydrogeological characteristics [10]. It is noteworthy that similar slow deformation patterns have also been documented in highly weathered metamorphic rocks, with no clearly defined rupture surface at depth [11,12].

Slope stability is conventionally assessed using the factor of safety (FoS), defined as the ratio of the resisting forces (shear strength parameters) to the driving forces (shear stress) along the base of the potential failure surface [13,14]. External forces that reduce the FoS can either initiate failure or accelerate the movement of an already initiated landslide.

Fundamental to slope stability analysis is the application of critical state soil me-chanics and particularly the Mohr–Coulomb failure criterion [15]. This principle asserts that failure occurs when the combined application of shear and normal stresses on any plane within a material reaches a critical value, that is determined by the material’s shear strength parameters such as cohesion and the angle of internal friction.

Hydrological processes play a critical role in affecting the shear strength parameters of landslide materials [16]. The infiltration of water can alter pore-water pressure and consequently the effective normal stress within the slope [17]. In nearly fully saturated conditions, the increased pore-water pressure causes a loss of cohesion and reduction of the effective normal stress, potentially resulting in shear strength degradation [18]. Additionally, weathering and alteration processes in the landslide material initiate the formation of weakened zones, prone to failure [19,20].

The origin of the water that infiltrates into the landslide materials may arise from direct precipitation or anthropogenic irrigation and can infiltrate throughout the matrix of the landslide [21,22], or through preferential pathways (e.g. fissures and cracks) created by the landslide movement [23,24,25]. Additionally, land-use changes such as deforestation can diminish the stabilizing effects provided by vegetation growth. Vegetation not only promotes transpiration, which reduces subsurface water content, but also enhances mechanical stability via an extensive root network [18].

Beyond the aforementioned hydrological and mechanical factors, decreases in material strength may also be driven by dynamic stress variations induced during seismic events [26].

A variety of methodologies have been developed to investigate landslide prone areas [27], regarding their geomorphological characteristics, the physical and mechanical properties of the subsurface formations, the stability conditions, the rate of displacement etc.

Unmanned Aerial Vehicles (UAVs), more commonly referred to as drones, have become an indispensable tool in the collection of geospatial data, particularly in the domains of photogrammetry and Light Detection and Ranging (LiDAR). These technologies have advanced the study of landslides significantly, enabling precise mapping, monitoring, and analysis of terrain dynamics [28,29]. UAV photogrammetry captures high-resolution imagery, while LiDAR provides precise elevation data, even in densely vegetated areas. This combination enables the identification of potential landslide triggers and the assessment of slope stability. Additionally, UAVs play a significant role in monitoring temporal changes in landslide dynamics. By conducting repeated surveys, researchers can detect soil displacement, erosion, and other changes in terrain over time [30,31]. LiDAR, with its high degree of accuracy, is particularly effective in detecting subtle variations in elevation, which are crucial for understanding the progression of landslides [32,33]. Furthermore, the analysis of LiDAR data can provide information on the volume of debris and slope stability, which helps in the planning of mitigation measures [34].

Geophysical investigations provide crucial insights into subsurface structures, by measuring the spatial distribution of physical parameters across various geological formations [35,36,37]. More specifically, an integrated approach of electrical resistivity and seismic techniques is widely applied in complex environments such as landslide prone areas [2,38,39,40]. The integrated application of these techniques facilitates the delineation of the landslide mass, by defining its geometry, spatial distribution and the depth of the rupture zone.

The Electrical Resistivity Tomography (ERT) technique has proven to be a particularly advantageous method in landslide studies due to its non-invasive nature and its ability to generate high-resolution results of subsurface heterogeneities. Its contribution is valuable in identifying zones of varying water/moisture content, which are critical for understanding water's role in landslide dynamics [41]. Changes in resistivity can also reveal the interfaces between different lithological units, such as sliding material and stable bedrock [42]. These interfaces are critical for delineating the boundaries of landslides, their internal structure and assessing their stability [43,44,45]. Additionally, ERT has been found to be invaluable in mapping zones of high-water content or elevated pore pressure. These factors have been shown to play critical roles in the destabilization of slopes and the triggering of landslides [46,47].

In addition to monitoring, ERT has been extensively utilized for evaluating landslide susceptibility and assessing mitigation measures [45,48]. The integration of ERT data with geological and geotechnical information facilitates the development of robust models that predict areas at high risk of failure [41]. Furthermore, ERT has been utilized to assess the efficacy of stabilization techniques, including drainage systems and retaining structures [49]. By observing changes in resistivity before and after the implementation of mitigation measures, researchers can assess their impact on subsurface conditions.

Time-lapse electrical resistivity tomography (ERT) is a particularly powerful tool for monitoring active landslides [49,50,51]. By conducting repeated resistivity measurements over time, researchers can track changes in subsurface properties associated with variations in moisture content, pore pressure, and material deformation. This approach enables the detection of early warning signs of impending slope failure, such as increased water infiltration or the development of shear zones [23,50,52].

Integrating analysis of ERT data with other geophysical techniques, such as seismic methods or ground-penetrating radar, enhances the reliability of subsurface models [53,54,55] and provides a comprehensive understanding of landslide mechanisms. Seismic refraction tomography (SRT) and multichannel analysis of surface waves (MASW) are widely utilized techniques for subsurface characterization of landslides, through the calculation of the P-wave and S-wave seismic velocities [56,57,58] helpful to the delineation of the bedrock’s depth and the geometry of the shear zone surface. Additionally, seismic velocities are correlated with the elastic properties and mechanical behavior of the subsurface materials such as their strength and elasticity, which are fundamental to assessing slope stability. They can be used to define their spatial variability in order to identify weak zones, often supplementing geotechnical data and improving soil slope stability models [17].

2. Study Area

The study area is located in the northern central Evia Island (Greece), NE of Dirfy mountain, at the proximity of Stropones community (Fig. 1). The broader region exhibits a rugged morphological relief, with slopes exceeding 30° locally. The geomorphological conditions and the subsurface hydrogeological and tectonic regime probably led to the formation of an unstable elongated mass in the area, with an SSW-NNE long axis that dips to the NNE, and with sloping sides bounded by the two main streams in the area (Loutsa and Miloi). According to the classification of landslides proposed by Varnes [59], this mass has the characteristics of a pre-existing translational landslide with soil creep in several parts..

During the past decades this slide phenomenon has caused damages to infrastructure (e.g. local road cracks) and building cracks and fractures [60,61,62,63], especially in the proximity of Stropones village, as shown in Figure 1.

2.1. Geology

According to the updated geological map (Figure 2) by Fotiadis and Karras [64] the study area mainly comprises Alpine formations that belong to the Pelagonian geotectonic unit of the Internal Hellenides. The greatest part of the study area is structured by a highly heterogeneous tectonic mélange formation (Pc.mg,fl). This formation is primarily composed of phyllites, schists, meta-sandstones, meta-graywackes, quartz and granite conglomerates of the Pelagonian Paleozoic basement. Moreover, within the tectonic mélange formation, several tectonic bodies of Paleozoic granites (Pz.Cs.γ), Liassic marbles (Ji.k,c) up to 1.5km sized, commonly known as Stropones marbles and smaller bodies, fragments of Carboniferous-Permian black micritic limestones, are also observed. The tectonic bodies of granites (Pz.Cs.γ) are strongly mylonitized and fragmented, while the Stropones marbles (Ji.k,c) either as massive or thick-bedded formations, include small Megalodontidae fossils. This tectonic mélange was formed during post-Eocenic compressive phases of the internal Hellenides, linked with the Alpine orogenetic cycle [65]. It is the structurally lowest formation in the area of interest, underlying all the other geological formations of the same unit.

The tectonic mélange formation is overlain by Lower to Middle Triassic platy limestones (Ti-m.k), which appear to be folded and underlying a volcano-sedimentary series (Tm.tf) [66]. This series consists of lavas and pyroclastics alternating with red silex siltstones, followed by turbiditic formations rich in olistostromes of pumiceous to amygdaloid and porphyritic lavas of andesitic and dacitic composition. The volcano-sedimentary series are overlain by a Late Triassic to Early Jurassic neritic carbonate platform. The base of this carbonate platform is comprised by white-gray dolomites and grey reefal dolomitized limestones (Ts-Ji.d), which transit to unbedded and medium-thick-bedded, gray neritic limestones (Ts-Jm.k). During the Middle Jurassic period, this continental margin was overthrusted by an ophiolitic nappe, primarily consisting of schistose serpentinites (σ). This overthrust is linked to the Eohellenic orogenetic phase. Subsequently, the ophiolitic nappe and the carbonate platform were overlain by Upper Cretaceous neritic limestones, followed by the deposition of the Middle-Upper Maastrichtian calcareous flysch formation [65]. These last formations are not observed in the proximity of the study area.

Regarding the post-Alpine sediments, they are mainly observed in the northern central part of the study area, in the proximity of the Stropones community. They primarily consist of talus cones (H.sc), comprised by unconsolidated coarse-grained angular and fine-grained loose Holocene sediments, which transit to conglomerates near the streams, originating from the surrounding formations.

2.2. Climate Characteristics

The climatic characteristics of a region with a landslide occurring represent a critical factor that must be considered. The temperature of a region, the atmospheric precipitation and the evaporation regime have a direct impact on the rate of infiltration in an area and thus the rate of aquifer recharge. Atmospheric precipitation determines also to a large extent the type and density of vegetation cover, which has a major influence on the rate of infiltration [67]. On the other hand, the relative humidity of the atmosphere plays a significant role in dictating the nature and extent of weathering of the bedrock [68].

In the context of the study of the hydrogeological conditions of the wider area of Stropones community, and due to the unique climatic conditions, a meteorological station was installed at the area in order to record the atmospheric precipitation, relative humidity, air temperature, wind speed and direction, barometric pressure, UVA and PAR. In the diagram of Figure 3, the cumulative monthly rainfall is presented (in mm), as recorded by this station during the period from March 2023 to March 2025.

2.3. Hydrogeology

The Stropones watershed occupies an area of almost 10 km². The development of the hydrographic network is controlled by the tectonic regime of the area, which also determines the conditions for the development of the existing underground karstic aquifer in the carbonate formations to a very significant extent. The carbonate formations, in terms of hydrogeology, generate high-capacity aquifers due to their considerable spatial development and high-water permeability. Furthermore, in the study area, due to the intense tectonism and the resulting fragmentation of the formations, secondary porosity develops also in the Middle Triassic volcanic-sedimentary series and the main mass of the tectonic mélange. The observed, relatively sparse hydrographic network of the area indicates the high rate of water infiltration to the subsurface geological formations. This process creates preferable water flow paths enhancing the weathering and alteration processes of the geological formations. The area is also characterized by the presence of multiple water springs, one of which is located to the South (Figure 4), close to Stropones community and has a significant flow

Accordινγ το the standards established by the International Association of Hydrogeologists, the geological formations of the area are categorized into three primary categories and six secondary categories, according to their degree and type of hydro-permeability, their spatial distribution, the efficiency of the water table they develop, their lithology and their stratigraphic structure [69]. These are (Figure 4):

• Cohesive (mainly carbonate) formations with secondary permeability.

• Extensive and high productivity aquifers. This subcategory includes all the carbonate formations in the study area (Ji.k,c, Ti-m.k, Ts-Ji.d, Ts-Jm.k).

2.

Granular or fissured formations with limited or no groundwater concentration.

a)

Locally significant groundwater, primarily in zones of fracturing and weathering of cohesive formations. This subcategory includes the volcanic-sedimentary series of the Middle Triassic (Tm.tf), the Middle Jurassic ophiolitic nappe (σ) and the main mass of the tectonic mélange (Pc.mg,fl). In this area, due to the intense fragmentation, the development of secondary porosity and atmospheric precipitation, the formations of this subcategory are also aquiferous.

b)

Strata with essentially no groundwater resources. This subcategory includes Palaeozoic granites (Pz.Cs.γ).

These amounts of water infiltrate the local formations, gradually saturating them, influencing the water-induced instability by increasing the pore-water pressure, alternating the effective stress regime and thereby further degrading the shear strength parameters of the formations.

3. Methodology

3.1. Remote Sensing Techniques

To monitor changes in the broader study area over time and systematically document alterations occurring at both residential and environmental levels, a series of historical aerial photographs was acquired from the Hellenic Military Geographical Service.

The selection of aerial photographs was based on the availability of the flight dates, image resolution (determined by flight altitude, scale and cloud cover) and the overall suitability for the analysis. Consequently, aerial photographs for the years 1945, 1971, 1988 and 1996 were chosen. For each year, at least three aerial photographs were selected with sufficient overlap. Using the Erdas Imagine software and by applying the photogrammetry technique, orthophotomosaics and digital terrain models (DTMs) were generated for each selected year.

Subsequently, an unsupervised classification was applied to each orthophotomosaic using ArcGIS Pro. The classification process involves assigning each pixel to a specific land use or land cover category based on its spatial location and radiometric value (solar reflectance in numerical form). In the case of the unsupervised classification, a specialized algorithm identifies distinct spectral patterns within the image, assuming that these patterns correspond to different land cover classes. The classification results were subsequently manually validated, refined and interpreted, incorporating field observations to accurately assign spectral clusters to meaningful land cover categories.

Further analysis of the classification outputs, combined with specialized image processing techniques, enabled a quantitative assessment of land cover changes over time. This was achieved by computing differences in the number of classified pixels per category and evaluating changes in class areas across different years.

The detailed investigation of the area of interest was achieved through the integration of UAV-based photogrammetry and LiDAR techniques, in order to generate high resolution orhophotomosaics and digital elevation models (DEM), revealing topographic features underlying the vegetation.

UAV photogrammetry, a process that involves the capture of overlapping images of the Earth's surface using cameras mounted on UAVs, is a critical component of this technological advancement. These images are processed through sophisticated techniques such as Structure from Motion (SfM), a method that utilizes image-based information to create georeferenced orthophotos, Digital Elevation Models (DEMs), and 3D point clouds [28,31]. This method provides high-resolution spatial data, making it invaluable for landslide investigations. Conversely, LiDAR utilizes laser pulses to measure distances and generate detailed 3D models of the terrain, capable of penetrating vegetation to reveal underlying topographic features [32,33]. The integration of these technologies has proven to be particularly effective in landslide studies, offering complementary datasets that enhance the accuracy of susceptibility models and risk assessments [28,34].

The DJI Matrice 350 RTK, equipped with integrated GNSS (Global Navigation Satellite System) that supports GPS (Global Positioning System), GLONASS (Global Navigation Satellite System), BeiDou, and Galileo satellite arrays, was utilized in this study. The aircraft's position was determined through the use of high-precision corrections in real time, facilitated by simultaneous connection via GPRS mobile telephony to a network of permanent reference stations. This approach enabled accurate determination of the aircraft's position. The Matrice 350 RTK is equipped with the DJI Zenmuse L1 camera, which consists of the LIVOXTM LiDAR module with three returns of the signal (laser beam), as well as one RGB camera with 20 MP resolution and a focal length lens of 8.8/24 mm.

The UAV was programmed to maintain a constant height of 100m from the ground (Follow Terrain mode). Additionally, the camera angle was set at 60° and 45° from the horizon, for the photogrammetry and LiDAR flights correspondingly, acquired in April 2024. The flight speed during data collection was set at 5.0 m/s. The percentage of front and side overlap was set at 70% and 50%, respectively.

According to the aforementioned characteristics, the ground sampling distance (GSD) was determined to be 2.73 cm/pixel for the Ortho and 3.86 cm/pixel for the Oblique. The total area covered was 0.82 km² and the density of the final point cloud was determined to be 4405 points/m².

The photogrammetric processing of the RGB aerial photographs and the LiDAR data was conducted with the DJI Terra software. The photogrammetric processing included alignment of the photos based on their spatial information and generation of the point cloud and the final orthophotomosaic (resolution 2 mm). Moreover, the accuracy of the final photogrammetric products was further enhanced by incorporating sixteen ground control points (GCPs), which were measured on the field through an RTK GNSS system. Regarding the LiDAR data processing, a point classification was conducted based on the returns of the LiDAR pulse (Figure 5), providing the ability to remove vegetation in the final DEM of the scanned area. The generated DEM was further processed using the slope analysis tool of the ArcGIS Pro software, enabling to highlight areas of steep slopes, prone to landslide activation.

3.2. Surface Geophysical Investigation

3.2.1. Electrical Resistivity Tomography

The underlying principle of ERT is that distinct electrical resistivities are exhibited by different geological materials based on their physical and chemical properties, such as porosity, saturation, and conductivity. Electrodes placed along the ground surface inject electrical currents into the subsurface, and the resulting potential differences are recorded at other electrode locations. These measurements are then converted into apparent resistivity values, which are processed using advanced inversion algorithms to produce two-dimensional (2D) or three-dimensional (3D) resistivity models [48].

The ERT measurements were acquired across twelve (12) sections (Figure 6), located at several places in the study area, based on the accessibility of the area, since there was dense vegetation and building covering the area. The electrode spacing varied between 5 to 10 m, due to the available spatial parameters. The ERT measurements were carried out with an IRIS Instruments Syscal Pro switch 48 resistivity unit, using 48 electrodes. The data was collected employing three distinct electrode configurations-the Wenner, Wenner-Schlumberger and Dipole-Dipole-in order to assess their differential sensitivity to variations in the subsurface electrical resistivity distribution. The sequences for the various electrode configurations for each ERT profile were then generated using the Electre Pro software, with this sequence determined by the geometric characteristics of each profile. Thereafter, the sequences were imported into the resistivity meter to initiate the automated data acquisition process.

The geoelectrical data underwent processing utilizing the Res2DInv software developed by Geotomo. Prior to the primary processing, a thorough quality control assessment of the data was conducted. This assessment entailed the identification and removal of points characterized as noisy or inconsistent, i.e., bad points. For solving the inverse problem, a modified version of the Gauss-Newton least-squares equation was applied [70]. This modification was achieved through the implementation of the smoothness-constrained least-squares method [71]. Meanwhile, the finite element method was implemented to calculate the forward model [72].

3.2.2. Seismic Method

In the context of the present study, P-wave and S-wave velocities were determined by applying the Seismic Refraction Tomography (SRT) and Multichannel Analysis of Surface Waves techniques. SRT technique relies on the measurement of direct and refracted travel times of seismic waves generated by a controlled seismic source (e.g. sledgehammer), at multiple locations (shotpoints) along the survey line. The generated waves propagate in all directions through the subsurface and are either transmitted through the surface layer producing direct arrivals, or are refracted at deeper layers, where they return to the surface due to the critical refraction phenomenon occurring at different lithological interfaces, producing refracted arrivals. Seismic wave propagation is recorded at the surface using an array of receivers (geophones), evenly spaced along the survey line. Direct and refracted first arrival times are analyzed to construct travel-time curves, based on which an initial velocity model is created. Subsequently, an inversion algorithm is applied, allowing for the iterative refinement of the model’s parameters [73]. This process continues until an optimal velocity model is achieved, minimizing the discrepancy between the measured and calculated travel times.

Multichannel Analysis of Surface Waves (MASW) technique utilizes the dispersive characteristics of surface waves, where their propagation velocity and penetration depth are determined by their frequency component [74]. Therefore, the frequency-dependent dispersion of surface waves reveals the material properties of soil layers at varying depths. Typically, MASW technique involves the application of a wavefield transform method for the transformation of seismic records from the time-distance domain to the frequency domain. In this manner, an image of the energy density distribution of surface waves at various frequencies is generated and the fundamental dispersion curve is identified based on the maximum energy amplitudes. Finally, through mathematical inversion of the dispersion curve, the 1D S-wave velocity distribution is derived, providing an average characterization of the subsurface covered by the geophone spread, corresponding to the midpoint of the receiver array [75].

Seismic data acquisition was conducted along a 57.5 m seismic line using 24 geophones with a natural frequency of 4.5 Hz, spaced at 2.5 m intervals. Seismic waves were generated by a 6.5 kg seismic sledgehammer at thirteen shot points distributed along the seismic line. Two of the shot points (outshots) were positioned at 30 m offset on either end of the seismic line, in order to increase the investigation depth. Additionally, two shot points were placed at a near offset distance of 5 m from the first and last geophones, to facilitate Rayleigh wave recording for the application of the MASW processing technique. The remaining nine shot points were placed in between the active spread at a 7.5 m shot interval, ensuring an adequate ray coverage through the subsurface for the reliable application of the SRT technique. Seismic waves were recorded by a 24-channel Geometrics SmartSeis seismograph at a 0.250 ms sampling interval and 512 ms record length.

The seismic data intended for the SRT technique was processed using the DW Tomo module of the Geogiga Seismic Pro software. All the seismic records were preprocessed by incorporating low-pass Butterworth filters and signal gain adjustments in order to enhance the first arrivals of the P-waves. The first break picking was conducted manually for all seismic records for the generation of the travel-time curves (Figure 7a). Afterwards, a gradient velocity model was calculated by applying the conventional intercept-time method [76], which was then employed as the initial model for the tomographic inversion approach. For the inversion of the seismic data, the smoothing-constrained regularized inversion approach [77] was implemented. The horizontal and vertical smoothing lengths were set to 10.00 and 3.12m respectively and a total number of 26 iterations were selected for the inversion. For the forward modeling solution, the shortest path method [78] was applied, with horizontal and vertical model cell dimensions equal to 1.250x0.625 m correspondingly.

The seismic record selected for MASW analysis corresponded to the shot point located at 5m near-offset from the last geophone. A specific processing window was selected in the seismic record in order to isolate surface waves from body waves (Figure 7b), optimizing the effectiveness of the wavefield transform method. The selected window was then transformed from the time-distance domain to the frequency-phase velocity domain by the application of the Fourier Transform. In the generated dispersion image (Figure 7c), the maximum energy amplitudes of the fundamental mode of the surface waves was picked in a frequency range of 7-22 Hz, composing the dispersion curve. Afterwards, an initial model was automatically calculated by incorporating the thickness option of the software. This option automatically determines the 1^st layers thickness according to the highest peaked frequency of the dispersion curve, while the thickness of the other layers is determined using the following equation: $H_{i+1}=H_i\times R$, where $H_i$ the thickness of the i^th layer and R a ratio equal to 1.2. Finally, the genetic algorithm inversion approach [79] was implemented to invert the dispersion curve, resulting in the calculation of the 1D Shear wave velocity distribution, corresponding to the midpoint of the seismic array.

3.3. Density Determination

The density of the geological formations of the study area is essential in order to calculate the elastic moduli. Geological samples were collected from outcrops during field campaigns and especially across Seismic Line 402 in order to conduct the essential laboratory measurements and establish the density of the geological formations [53,80,81]. Three weight measurements were taken for each sample, including: (i) the dry weight in air (W1); (ii) the saturated weight in air (W2); and (iii) the saturated weight while submerged in water (W3). From these measurements, the dry bulk density (ρ_d), saturated bulk density (ρ_s) and granular density (ρ_g) of the samples were derived. A total of 247 geological specimens were weighed using the Kern EMB 100-3 precision scale with an accuracy of 0.001 g (Figure 8). Measurement errors were assessed by analyzing multiple samples from each geological formation, considering the median characteristic value. The calculated densities along with the measurement errors for each geological sample are presented in Table 1.

4. Results and Interpretation

4.1. Remote Sensing

The results of the unsupervised classification process applied to the orthophotomosaics from 1945 and 1996 are presented in Figure 9. The classification process categorized the landscape into four distinct classes: Class 1 and Class 2 correspond to areas with different types of vegetation, including forested areas, shrubland and riparian vegetation. Class 3 represents areas with low vegetation cover and Class 4 encompasses non-vegetated surfaces, such as exposed rocky slopes, uncultivated fields, roads and urban areas.

The comparison of the classification between the two orthophotomosaics (Figure 9) reveals a significant increase in the frequency of occurrence of the areas classified as Class 3 and Class 4. The numerical processing revealed that inside the region of the watershed (blue dashed line) of the study area the extent of Class 1 has decreased, from 1945 to, by 20,890 m² probably due to forest fires that have burned the area and deforestation processes.

Additionally, the DEM generated from the LiDAR survey in 2023 is presented in Figure 10. After the removal of the vegetation, the morphological relief was highlighted in great detail, revealing the boundaries of the unstable mass. The study area was divided into four discrete morphological blocks, each exhibiting relatively consistent elevation levels. These blocks were separated by steep slopes (greater than 45°), which resulted in a general degradation of elevation towards the North. This classification was conducted subsequent to the slope analysis in the LiDAR-generated digital elevation model (DEM) and on-field observations (Figure 10).

4.2. Geophysical Investigation

4.2.1. ERT Results

The results of the ERT technique are presented in three different fence diagrams (Figure 11, Figure 12 and Figure 13), one for each block as recognized by the morphological analysis explained in section 4.1. All the ERT profiles presented in fence diagrams are conducted on the highly heterogeneous tectonic mélange formation (Pc.mg,fl) of the creeping mass. In the first and second block the subsurface is characterized by a wide range of resistivity values (20-5000 Ohm.m), in contrast with the third block, where a lower resistivity value range is observed (20-400 Ohm.m). For that reason, a different color scale was used for the third block.

Concerning the ERT fence diagram of the first block (Figure 11), high resistivity values (>600 Ohm.m) dominate throughout the subsurface. These resistive formations are attributed to the cohesive blocks of the Paleozoic granites (Pz.Cs.γ) and/or the Stropones marbles (Ji.k,c), located as parts of the tectonic mélange that develop the morphological outcrop of the third block. On the contrary, lower resistivity zones (<400 Ohm.m) are interpreted as the main mass of the tectonic mélange, primarily comprised by phyllites and schists. Across section 301 (Figure 11), the conductive zone seems to occupy the largest part of the investigated surface in contrast to the other sections of the block. The resistive zones seem also smaller in extent across this section. These conductive zones and especially their near-surface part (< 10m depth) may be characterized by a high percentage of moisture content, based on the field observations.

Moving to the North, a similar subsurface structure is observed in the ERT fence diagram of the second block (Figure 12), with relatively resistive formations (>600 Ohm.m), related to blocks of granites and marbles of varying sizes, alternating with the more conductive formations (<400 Ohm.m) corresponding to the main mass of the tectonic mélange. In this case, the resistive formations that occupy the morphological outcrop of the second block seem constrained in spatial distribution, in contrast to the first block. The more conductive formations dominate, especially at areas of lower topography, on the lateral sides of the morphological outcrop that extend towards the two streams of the area, Miloi and Loutsa, located on the eastern and western part of the area correspondingly. Therefore, it seems that no resistive blocks of marbles or granites are located in these areas, which could affect the rate of the creeping mass. The conductive zone that prevails in the area of this block, could be related to the enhanced instability observed in situ during the field measurements.

Farther to the North, closer to Stropones village, the subsurface regime of the third block (Figure 13 and14) appears to be different in comparison to the two southern blocks described above. The resistivity values range is constrained to 20-400 Ohm.m, indicating that the subsurface is mainly dominated by the more conductive formations, related to the main mass of the tectonic mélange. The resistive zones of this block do not practically correspond to the values of the two previous blocks and, therefore, probably not to blocks of marbles and/or granites, since the color scale used here is narrowed down. Furthermore, the third block is located directly downwards of a water spring that infiltrates water to the subsurface and consequently downgrades the resistivity of the subsurface formations. Blocks of marbles and granites are not expected here based on the resistivity values or have been masked due to the increased subsurface water.

A diachronic investigation between the wet and dry periods (April – October 2023) was conducted with the acquisition of time-lapse data across three sections (101#201, 102#202, 103#203). An indicative example of the time-lapse results is illustrated in Figure 15. This specific section was selected due to its vicinity to the prementioned spring, expected to affect the most the subsurface status. Positive changes indicate reduction, while negative changes an increase in the subsurface resistivity values, at the end of the dry period. The most significant percentage changes (plus 20-80%) in the resistivity values are observed only in the central and very near-surface part of the section and are not evaluated as important due to their small and surficial extent. The main part of the subsurface area presents changes in resistivity up to a maximum of 10%. This zone (Figure 15) has already been related to the main mass of the tectonic mélange. The very small variations (close to zero) of the resistivity indicate a relatively stable hydrogeological regime, with a constant presence of subsurface water throughout the year.

On the other hand, there is a near-surface (max 10 m depth) zone with an increase in resistivity values (minus 10-40%), which is relatively expected due to the dry period. Most of the resistivity variations are observed in the upper 10m of all time-lapse sections and can be related to pore-water content variations. Thus, the increase in pore-water pressure, during the wet period, can cause a loss of cohesion and reduction of effective normal stress, potentially resulting in shear strength degradation of this surficial zones, which destabilizes the creeping mass.

4.2.2. Seismic Method Results

In Figure 16 the SRT-402 profile and the 1D shear wave velocity distribution of the MASW technique (black line), corresponding to the midpoint of the profile, are presented. The reliable part of the SRT profile is defined by the active spread length covered by the geophones (black symbols on the surface). Based on the distribution of seismic velocities, three seismic layers can be identified, indicating the presence of subsurface formations with different characteristics. The first, near-surface seismic layer is characterized by seismic velocities of Vp = 300-700 m/s and Vs = 260 m/s, with a thickness ranging from 2.5 to 4.0 m. The second seismic layer exhibits seismic velocities of Vp = 900-1800 m/s and Vs = 400 m/s and 2.5-6.5 m thickness. Finally, the third seismic layer is defined by seismic velocities of Vp > 1800 m/s and Vs > 600 m/s, with average values of Vp = 2100 m/s and Vs = 730 m/s, from the top of the layer until the maximum investigation depth. This formation is encountered at depths greater than 9.5, 5.0 and 7.5 m correspondingly at the beginning, mid and end part of the active spread length.

Taking into consideration the geological condotions (Figure 2) and field observations conducted during the data acquisition, the seismic layers are interpreted as follows. The first, near-surface, seismic layer is attributed to the presence of loose, unconsolidated materials of the colluvial soil deposits, due to the low seismic velocity values that characterize it. Subsequently, the second and third seismic layers are interpreted as the main mass of the highly heterogeneous tectonic mélange (Pc.mg,fl) formation. More specifically, the second seismic layer may be attributed to the weathering layer of this formation, predominantly composed of highly weathered schists, while the third layer, characterized by higher seismic velocities, corresponds to the more compacted and slightly weathered section of the tectonic mélange formation (Pc.mg,fl).

Based on the calculated P-wave and S-wave velocities and the density determination from the laboratory measurements (section3.3), the elastic moduli of the investigated geological formations across seismic line 402 have been calculated. In Table 2, the equations used for the elastic moduli calculation are presented.

Poisson’s ratio (ν) represents the ratio of transverse to axial strain in response to applied stress. It is a dimensionless parameter, with values typically ranging between 0.1 for compact formations and 0.5 for loose formations, strongly affected by factors such as saturation and in-situ stress conditions. Young’s modulus (E) quantifies the material’s stiffness by expressing the ratio of longitudinal stress to longitudinal strain under uniaxial compression or tension. Similarly, the shear modulus (G) expresses the ratio of shear stress to shear strain. Both Young’s modulus and the shear modulus serve as fundamental parameters in soil mechanics, describing the material’s stiffness. They are widely utilized in soil consolidation studies and numerical modeling of landslide processes. Additionally, the bulk modulus (K) describes a material’s resistance to uniform compression by relating volumetric strain to applied hydrostatic pressure.

The elastic moduli were calculated for the second and third seismic layers, which correspond to the weathering layer and the more compacted section of the tectonic mélange formation (Pc.mg,fl), respectively (Table 3). Due to the presence of the subsurface water, mainly originating from the nearby springs and the high levels of precipitation, the saturated bulk density was considered to be more appropriate. At this point we remind that the tectonic mélange comprises multiple lithologies, including Stropones marbles (Ji.k,c), tectonic bodies of Paleozoic granites (Pz.Cs.γ), schists (Pc.mg,fl), phyllites (Pc.mg,fl) and black micritic limestones (Pc.mg,fl). Therefore, to ensure a representative density value for the whole mélange formation, the average saturated density of these lithologies (Table 3) was adopted, while for the weathering layer of the formation, the calculated saturated density of the weathered schists was used (Table 3). Regarding the seismic velocities, the average Vp and Vs values of these two formations were incorporated for the calculations. The elastic moduli of these two layers are presented in Table 3.

4. Conclusions

The interpretation of the acquired surface geophysical data provided new information regarding the movement of the mass, including the hydrogeological aspects. For example, several subsurface water flow paths have been revealed. According to the meteorological data, rainfall is observed almost all times of the year, apart from a couple of months during the summer, which means that there is a relatively continuous, subsurface water flow. The depth of the rupture zone of the creeping mass is delineated at the first five to ten meters from the surface, especially from the difference in diachronic resistivity change.

Apart from that, the deforestation of the surrounding mountainous area, due to forest fires or anthropogenic interference has been revealed that probably plays an important role in the landslide /creeping evolution.

The determined elastic moduli of the formations could be useful in new constructions in the area.. Moreover, the adumbration of existing cohesive zones, interpreted as blocks of marbles and granites, revealed the structure diversity of the mass.

As a first step in mitigating the effects of surface creeping, water regulation in terms of natural surface water flows, irrigation practices and wastewater management (from homes and other uses) could be considered. The next step should be to drill a few tens of metres deep to confirm the results of the remote sensing geophysical techniques.