Morphological Changes Detection of a Large Earthﬂow Using Archived Images, LiDAR-Derived DTM, and UAV-Based Remote Sensing

: In mountainous landscapes, where strongly deformed pelitic sediments outcrop, earthﬂows can dominate denudation processes and landscape evolution. This paper investigated geological and geomorphological features and space-time evolution over a 65-year time span (1954–2019) of a large earthﬂow, representative of wide sectors of the Apennine chain of southern Italy. The landslide, with a maximum length of 1.85 × 10 3 m, affects an area of 4.21 × 10 5 m 2 and exhibits two source zones: a narrow and elongated transport zone and a lobate accumulation zone. Spatial and temporal morphological changes of the earthﬂow were assessed, comparing multi-source and multi-temporal data (aerial photographs, Google Earth satellite images, Light Detection and Ranging (LiDAR) and Unmanned Aerial Vehicles (UAV) system data). Geomorphic changes, quantiﬁed using Digital Terrain Models (DTMs) of differences, highlighted an extensive lowering of the topographic surface in the source area and a signiﬁcant uplift at the landslide toe. Moreover, the multi-temporal analysis showed a high increase of landslide surface (more than 66%) during the last 65 years. The volumetric analyses showed that different sectors of the earthﬂow were active at different times, with different rates of topographic change. Overall, the used approach highlighted the great potentiality of the integration of multi-source and multi-temporal data for the diachronic reconstruction of morphological landslide evolution.


Introduction
Earthflows are among the most common and widespread mass movement phenomena in nature; they occur in many of the world's hilly and mountainous areas [1][2][3] where they represent the primary agent of erosion and contribute large amounts of sediment to streams and rivers [4]. Such type of landslides can cause extensive damage to property and infrastructure [5][6][7].
Earthflows generally exhibit a wide, amphitheatre-like source zone, a narrow and elongated transport zone, and a lobate depositional area characterized by multiple compressional features [2,5,8,9]. Most earthflows move primarily in the source zone by sliding on discrete basal and lateral failure surfaces [10]. They occur in fine-grained soils that consist dominantly of plastic silt or clay as well as rocky soils that are supported by a plastic silt-clay matrix and progress gradually by brief episodic movements or periods of sustained, relatively steady movement [10][11][12][13]. Long-term persistence of slip surfaces at residual strength contributes to earthflow reactivation [14]. Mechanical and hydrologic isolation of earthflows by clay-rich layers at their boundaries also contribute to the persistent instability of the phenomena [15]. Typically, pulses of movements characterize the temporal behavior of the earthflows, which exhibit relatively short periods of intense activity followed by longer quiescent periods of reduced or no measurable movement [5,16,17]. Earthflows typically move at slow to moderate rates ranging from 10 −10 m/s to 10 −3 m/s; in a few instances, earthflows have surged briefly to rates ranging from 10 −2 m/s to about 10 −1 m/s [2].
The insight into the space-time activity of earthflows is essential for accurately recognizing their morphological evolution [18] and for improving knowledge of the mechanisms of deformation of mass movements and can contribute to improved risk assessment and mitigation [19,20].
During the last decades, the integration and development of several photogrammetric and of geographic information system (GIS) techniques have increased the potential use of remotely sensed data for mapping and assessing the space-time evolution of landslides [21][22][23]. In this respect, historical aerial photographs proved to be a useful tool, not only for qualitative analysis of the geomorphic processes but for a quantitative assessment as well. For instance, the multi-temporal interpretation of photographic images provides the systematic examination of the progressive evolution of landslides and may lead to a better understanding of their causal factors [24][25][26].
An important step forward was the use of multi-temporal Digital Terrain Models (DTMs) to detect the topographic changes of landslides as well as to estimate displacement rates and/or to estimate the volumes of removed or accumulated material [23,27]. Digital photogrammetric techniques to investigate and monitor landslides have been widely employed by means of different tools in the optical spectrum such as Differential Synthetic Aperture Radar Interferometry (DinSAR), Light Detection and Ranging (LiDAR), Global Navigation Satellite Systems (GNSS), and the digitization of old topographic maps [27][28][29][30][31][32][33][34][35].
In the last few years, among the remote sensing techniques, there has been a rapid development of the application of Unmanned Aerial Vehicles (UAV) for landslide studies [36][37][38][39][40]. The advantages of UAV-based remote sensing are the following [41]: (i) realtime applicability; (ii) flexible survey planning; (iii) high resolution and low cost. Several landslide studies highlighted the idea that the comparisons of ortho-images and/or high resolution DTMs, besides the space-time changes of landslide area and of the estimate of the volume of accumulated or removed material, allows us to detect and analyze the spatial distribution and evolution of landforms related to the landslide activity, e.g., scarps, tension cracks, counter slopes, shear zones, trenches, displacement vectors, etc. [27,38,39,42].
The aim of this work was to analyze and better understanding the geological and geomorphological features and the space-time evolution of a large earthflow (i.e., Vomice landslide system) located in the north-eastern sector of the Calabria region (southern Italy) over a time span of 65 years (1954-2019). In particular, the analysis stressed the planoaltimetric differences caused by the geomorphological evolution of the landslide, providing an estimation of the mobilized volumes at different time windows. The research was carried out by means a geomorphological multi-temporal analysis, integrating different remote sensing data (e.g., historical and recent aerial photographs, Google earth images, LiDAR-derived DTM, and UAV survey) and field data, (geological, geomorphological, and GPS surveys).

Study Area
The Vomice landslide system is located on the right side of the Straface River, NE sector of the Calabria region (southern Italy) (Figure 1a,b). The study area ranges in elevation from 710 to 160 m a.s.l. Slope gradients vary between 0 to about 45 • , with an average value of 15 • , whereas the predominant slope aspect is towards north and east. Climate in the study area is of Mediterranean type, with hot and dry summers and mild and wet winters [43]. Pluviometric data between 1950 and 2019 coming from the station of Albidona (Figure 1a), located about 5.5 km SW of the study area (4,420,430 N, 625,482 E, 790 m a.s.l.), showed an average annual precipitation equal to 746.6 mm, distributed on 77 rainy days (Figure 1c). For the period analyzed, the maximum value of annual precipitation was 1419.3 mm, recorded in 1959, whereas the minimum value was 276.8 mm, observed in 2001. The major precipitation events are mainly concentrated in the period from October to March, with tion was 1419.3 mm, recorded in 1959, whereas the minimum value was 276.8 mm, observed in 2001. The major precipitation events are mainly concentrated in the period from October to March, with rainfall peaks generally in December and January, while the period between June and September is the driest period of the year [44]. The maximum value of daily precipitation was 244.4 mm, registered on 23 December 1990   The geological setting consists of Mesozoic-Cenozoic sedimentary rocks of Ligurid and Sicilid units [45,46]. These lithological units are crossed and cut by Quaternary highangle strike-slip and extensional faults [47,48].
The landscape of the study area reflects the complex interplay between the geological and structural arrangement. The areas where outcrop rocks of different composition and erodibility gave alternatively steep slopes in contrast to typically rounded and gentler slopes. Additionally, the variety of outcropping lithologies, often controlled by fault systems, determined the development of morpho-structural ridges bounded by fault scarps, saddles, straight ridges, and straight channels. This portion of Calabria is strongly affected by both landslides and water erosion processes because of its geological and climatic characters [44,[49][50][51]. In particular, mass movements assume a remarkable importance in the morphological evolution of the reliefs, and complex landslides and earthflows are the most common type of mass movements.

Material and Methods
A multi-temporal analysis, using remote sensing data coupled with field surveys, was applied in order to map the geological and geomorphological features and investigate the morpho-evolution of the Vomice landslide system in the period between 1954 and 2019.
The analysis was performed considering information arising from different cartographic databases along with data derived from remote sensing images and Digital Terrain Models (DTMs). In particular, the used data were: (i) topographic map related to 1954; (ii) aerial photographs (dated 1954, 1984, 1990), (iii) orthophotos (dated 2001 and 2007), (iv) Google Earth satellite images (dated 2012 and 2016), (v) a LiDAR-derived DTM (dated 2012), and (vi) UAV-based photographs (dated 2019). The sequence of archived images was also used for reconstructing the activity state of the landslide.
In order to perform the photogrammetric orientation of the UAV images, georeference the archived images, and evaluate the accuracy of and DTMs, a set of ground control points (GCPs) were collected. The GCP survey was performed using a GNSS Leica 1200™ receiver, in Real Time Kinematic (RTK) mode, with an average accuracy of about 1 cm. The GCPs were surveyed based on targets such as building corners and/or road intersections, which can be easily recognized in the aerial photos, and 50 × 50 cm artificial targets placed on the ground during the Drone flight.
Finally, all the collected data were managed and elaborated using QGIS 3.6 software. Through QGIS software, the aerial photographs were georeferenced using the WGS 84 projection datum and the Universal Transversal Mercator as a coordinate system, using a set of GCPs such as road intersections and house corners located outside the landslide area, which was assumed to be stable. A minimum of 10 GCPs for each image were used and GCPs were added until the root means square errors (RMSE) associated with the rectification process were <1 m. A second order polynomial (quadratic) transformation and a bilinear interpolation method were used to rectify all the aerial photographs. Quadratic transformations correct some errors associated with the curvature of the earth and the topography of the study area [52].

Collection and
Orthophotos dating from 2001 (1:10,000 scale) and 2007 (1:5000 scale), downloaded from the website of the Centro Cartografico della Regione Calabria (http://geoportale. regione.calabria.it/opendata), were collected and, to capture the recent planimetric changes, two satellite images dating from 2012 and 2016, available for viewing on Google Earth Pro, were used ( Table 1).
The DTMs available to perform the plano-altimetric and volumetric changes are summarized in Table 2. A DTM related to 1954 was derived by digitizing the contour lines and points of an official technical topographic map of the Calabria Region at 1:10,000 scale. The raster DTM, with a pixel size of 5 × 5 m, was obtained through a spatial interpolation using QGIS software. Furthermore, a DTM related to 2001, with a spatial resolution of 5 m, was downloaded from the website of the Cartographic Center of Calabria region (http://geoportale.regione.calabria.it/opendata), while a DTM of 2012, with a spatial resolution of 1 m, deriving by LiDAR scanning on an aerial platform, was acquired from the Italian Ministry for the Environment, Land, and Sea.

Acquisition and Photogrammetric Elaboration of UAV Images
The UAV survey was performed in February 2019 using a Parrot Anafi Drone equipped with a 21 Mpixel (5344 × 4016 pixel) RGB camera and an on-board GNSS system for the accurate geolocation of the acquired images. Due to the extension and morphology of the study area, it was necessary to split the survey area into three different flight missions with three different take-off points and with an appropriate overlap between flight areas; the different missions were performed with the same average flight altitude and frontal and side overlap. A total of 1122 nadir photographs, with a frontal overlap of 85% and a side overlap of 80%, were acquired ( Table 2). The average flying height of 131 m above the terrain guaranteed a ground sample distance (GSD) lower than 7 cm/pixel.
For photogrammetric orientation of the UAV images, a total of 45 GPS points were used. Thirty of these points were used as GCPs for the orientation process ( Table 2). The remaining were used as check points (CHKs) to calculate the RMSE of the UAV-threedimensional (3D) model.
As well as for the flight, the first elaboration step was performed by computing images of the three missions separately (but with the same methodology); only after a first elaboration were the three areas merged into a single model. The orientations of UAV digital images were computed by applying a global bundle block adjustment method [53][54][55] using Pix4D mapper software. The resulting 3D point cloud was composed by about 66.5 million points. The accuracy of the 3D-model is shown in Table 2. The RMSEs in the CHKs refer to the residuals calculated in these points after the bundle block adjustment, resulting in 0.19 m for X, 0.18 m for Y, and 0.21 m for Z ( Table 2).
The elaboration of the UAV 3D-model allowed us to produce an orthophoto of the entire study area (with a nominal ground resolution of 0.3 m × 0.3 m) and a DSM (with a resolution of 0.5 m × 0.5 m). Furthermore, a DTM with a ground resolution of 0.5 m × 0.5 m was obtained by the manual filtering of the point cloud, using Cloud Compare software (Version 2.6.1), in order to detect and remove points that corresponded with vegetation (e.g., trees or shrubs). Finally, the orthophoto and DTM were exported as raster files and further analyzed in QGIS.

Accuracy of the DTMs Collection
Before the space-time analysis of landslide evolution, the accuracy of every DTM must be calculated. To evaluate the errors of DTMs, a set of GCPs are needed for determining the difference in elevation (residual error) between the digital surface and the real surface of the same locations on the ground. This requires both a suitable sample of GCPs and suitable statistics from which to derive error terms. In this study, 15 GCPs for every DTM were used. They were located outside the landslide perimeter and were homogeneously distributed over the Vomice catchment. The stable terrain areas were mainly located in crest positions and in plots that did not suffer changes or levelling during the period considered.
The accuracy of the DEMs was assessed using an approach proposed by [56], which allows the determination of the mean error (ME) and the standard deviation error (SDE): where Z DTM is the measurement of elevation from the DTM and Z ref is the higher accuracy measurement of elevation for a sample of n points. ME can be either negative or positive and is used to distinguish between the unwanted systematic errors, whereas the SDE is used to distinguish between the expected and tolerable random errors [57,58]. Each GCP was used as a Z ref value to provide an estimation of error derived from the entire DTM.
The results of accuracy of the four DTMs are given in Table 3. Table 3. Accuracy of digital elevation models (DEMs), assessed by independent ground control points (GCPs) and expressed in terms of mean error (ME) and standard deviation error (SDE). Generally, the ME indicates an underestimation of the DTMs regarding GPS points, except for the UAV-DTM, which is slightly positively biased. In particular, the MEs were less than 0.5 m for all DTMs except for the year 1954 (1.12 m), whereas the SDEs ranged from 0.34 m for the DTM from 2019 to 4.05 m for the DTM from 1954 (Table 3).

Geological and Geomorphological Setting and Plano-Altimetric Changes Analysis of the Landslide
The geological and geomorphological investigations of the landslide area were performed through photo-interpretation of the high-resolution 2019 orthophoto generated by UAV survey, coupled with detailed field surveys performed between 2017 and 2019.
Geological surveys were used to preliminarily identify the rock types outcropping on the study catchment. Structural measurements of faults, fractures, and bedding planes were conducted to reveal the structural predisposition of slope to sliding. In particular, macroand mesostructural analyses were performed in the landslide area and its surroundings to reconstruct the geometry and kinematics of the fault planes.
The geomorphological surveys, coupled with the analysis of both high-resolution orthophoto produced after the UAV survey and a shaded relief map created from DSM raster data, were performed to identify and map the main surface features linked to gravity-driven processes-based on field recognition and freshness of the topographic signatures typical of gravity-related landforms [59], including scarps, conjugate scarps, uphill facing slopes, ridges and depressions, shear surfaces, ground cracks, levee, changes in the drainage network-and to infer the type of movement and state of activity according to [13]. Specifically, active landslides are those that are currently moving (including firsttime movements and reactivations). They are characterized by the freshness of their topographic signatures and were identified by multi-temporal geomorphological surveys carried out between February and April 2019.
Suspended landslides [11] are those that have moved within the last annual cycle of seasons (i.e., during 2018) but are not moving in 2019 (time of the last geomorphological survey). Indeed, the topographic signatures of the suspended landslides are modified by erosion and partially masked by vegetation.
Dormant landslides are those that last moved more than one annual cycle of seasons ago and whose causes of movement remain apparent (i.e., landslide bodies moved before 2018), their topographic signatures are obliterated by water erosion and masked by vegetation.
The photo-interpretation and comparison of multi-temporal ortho-images allowed us to reconstruct the planimetric changes of the landslide system in the last 65 years (1954-2019). Therefore, by means of interpretation of each ortho-images, eight maps of the complex landslide system were produced. These maps showed the landslide boundary and displayed how the landslide geometry changed with time and showed the orientation of longitudinal structures indicating earth-flow movement along the flow path. Moreover, the multi-temporal interpretation of ortho-images gives information about the landslide activity (active or dormant areas) occurring in the period considered. Specifically, the activity of the landslide portions was evaluated using several geomorphological criteria based on the morphological freshness of landslides [60], such as: presence of bare scarps or with poor vegetation cover, depletion and deposition areas well developed, irregular slope profiles, ground cracks, slope ruptures.
A GIS overlaying of these maps allowed us to quantify the evolution of landslide area and identify zones of reactivation, enlargement, or new movements. From this analysis, the retreat rate of the landslide crown, which occurred from 1954 to 2019, was computed.
In order to detect the topographic changes and related mass distribution due to the landslide activity during the period between 1954 and 2019, the DTMs related to 1954, 2001, 2012, and 2019 were compared. All DTMs used for these analyses were resampled with a cell size of 2 × 2 m to facilitate the quantitative comparisons.
Resampling was done using cubic-convolution methods [61]. However, it is important to highlight the fact that, because the morphological changes induced by the landslide activity were in the order of meters, biases between DTMs had a limited negative effect on the quantitative assessment of mass transfer along the slope. A cell-by-cell subtraction of the different DTMs provides a straightforward way of obtaining positive and negative elevation differences corresponding to accumulation or depletion, respectively. Raster

Geological Features of Landslide Area
The geological setting of the study catchment ( Figure 2), consists of Mesozoic-Cenozoic structurally complex [62,63] sedimentary rocks: the "Liguride Unit" [45,64] and the "Sicilide Unit" [65].  The Liguride Unit, outcropping in the SW portion of the study area ( Figure 2), is represented by a continuous sedimentary succession composed from bottom to top by the turbidite sequences of the Saraceno and Albidona Formation. The Saraceno Formation (Upper Oligocene-Aquitanian) is mainly made up of calcarenites and calcilutites intercalated by clays and sandstones. The Albidona Formation (Burdigalian) is mainly composed of alternating sandstone, marl, clayey marl, and silty clays.
The Sicilide Unit, which extensively outcrops in the central and north-eastern portion of the study area (Figure 2), includes both the Varicolored clay-shale Formation, composed by pelagic sediments of the Cretaceous age that shows characteristics of a strongly and pervasively deformed broken formation [62,66,67], and the Numidian Flysch Formation (Lower Miocene), outcropping only along a narrow NW-SE oriented strip, composed of silty clays with intercalated levels of quartz arenites.
The Varicolored clay-shale Formation (the so-called 'Argille Variegate' or 'Argille Varicolori') [65] that represents the lithological unit within which the Vomice landslide develops (Figure 2), from a lithological point of view, can be divided into two facies: a 'chaotic clayey' facies and a 'schistose' facies. The 'chaotic clayey' facies ( Figure 3a), widely outcropping in the catchment from the bottom to the top of the slopes, consist of a melange of pervasively deformed rocks showing "blocks-in-pelitic matrix" fabrics. The pelitic matrix is mainly represented by highly tectonized red/reddish brown and green/greenish grey shales with disrupted bedding, within which competent blocks and fragments of various sizes and shapes occur.
The latter consists of dismembered chert, limestone, and rare sandstone beds, as well as blocks of whitish marly limestones and greyish calcarenites (locally containing Nummulites). The 'schistose' facies ( Figure 3b,c), only locally outcropping at the bottom of the Vomice stream (in the central sector of the catchment) due to the deepening action of the run-off waters, includes green, red, and brown shales and marls (including centimetrethick grey marly limestone beds), with local intercalations of red shales and marls with interbedded calcirudites, calcarenites, calcilutites, marly limestones, and bedded cherts.
Along the drainage channels of the catchment and along the flanks of the valley, where the strongly deformed pelitic sediments of the varicolored clays outcrop, Holocene heterogeneous and heterometric chaotic gravitational accumulations linked to the Vomice complex landslide system, of variable thickness (from a few meters up to tens of meters), are found ( Figure 2). Moreover, in the upper part of the slope, immediately upstream of these landslide deposits, photo-interpretation analyses and field surveys highlight the presence of some slope ruptures and a trench linked to the presence of a dormant deep-seated rock slide that mainly affects the Saraceno Formation ( Figure 2).
From a tectonic point of view, in the study catchment, the pristine Late Oligocene-Pliocene low-angle tectonic contacts between the Sicilide and Liguride units (outcropping along the lower valley of the Straface River) have been sliced by Pliocene-Quaternary high-angle strike-slip and extensional faults, also clearly recognizable on a morphological basis and which controls the morphology of the catchment ( Figure 2).
Particularly, the southwest side of the study area is bordered by WNW-ESE striking left lateral and transpressive strike-slip faults, linked to the Saraceno shear zone [48], which acted as left-lateral during the Pliocene-Quaternary (Figure 2). Evidence of left strike slip motion has been found on some fault planes (Figure 3d). N-S and E-W trending dip-slip to slightly oblique Quaternary high-angle extensional faults, locally disposed in stair-step fashion, have been mapped on the central and north-eastern side of the Vomice catchment ( Figure 2). Kinematic analysis of slip lineations preserved on striated fault surfaces, outcropping within the schistose facies of the Varicolored clay-shale Formation, indicates a predominant extensional displacement (Figure 3e). In some cases, linked to the Quaternary high-angle faults, thick damage zones have been observed (Figure 3f).

Landslide Map and Main Features of the Vomice Landslide
Multitemporal geomorphological field surveys performed in the study catchment, coupled with the analysis of both high-resolution orthophotos produced after the UAV survey and DSM hillshade (Figure 4a,b), allowed us to map, in extreme detail, landslide bodies and related morphological features.
The obtained results are shown in Figure 4c. In the map, three states of activity [13] were distinguished for the landslide bodies: active, suspended, and dormant ( Figure 4c).
The Vomice landslide is characterized by a long and complex history described by [44]. This long history is characterized by a recent main re-activation that occurred in 2010 when the landslide toe reached the bottom of the valley. The landslide develops from an elevation of about 480 m a.s.l., which corresponds to the crown of the source area, to an elevation of 165 m a.s.l. within the Straface river valley, where the landslide toe is located ( Figure 4).

Landslide Map and Main Features of the Vomice Landslide
Multitemporal geomorphological field surveys performed in the study catchment, coupled with the analysis of both high-resolution orthophotos produced after the UAV survey and DSM hillshade (Figure 4a,b), allowed us to map, in extreme detail, landslide bodies and related morphological features.
The obtained results are shown in Figure 4c. In the map, three states of activity [13] were distinguished for the landslide bodies: active, suspended, and dormant ( Figure 4c).
The Vomice landslide is characterized by a long and complex history described by [44]. This long history is characterized by a recent main re-activation that occurred in 2010 when the landslide toe reached the bottom of the valley. The landslide develops from an elevation of about 480 m a.s.l., which corresponds to the crown of the source area, to an elevation of 165 m a.s.l. within the Straface river valley, where the landslide toe is located (Figure 4).
The unstable area forms entirely in the Varicolored clay-shale Formation (involving mainly clayey soils), and starts immediately downstream of the WNW-ESE striking left lateral and traspressive strike-slip faults linked to the Saraceno shear zone, where a perennial spring has been observed ( Figure 2).  The unstable area forms entirely in the Varicolored clay-shale Formation (involving mainly clayey soils), and starts immediately downstream of the WNW-ESE striking left lateral and traspressive strike-slip faults linked to the Saraceno shear zone, where a perennial spring has been observed ( Figure 2) The Vomice earthflow (Figures 4c and 5a) is a large and complex landslide system with a maximum length of 1.85 × 10 3 m, a width ranging between 40 and 300 m, and an area of about 4.21 × 10 5 m 2 (~42 ha).
The landslide is characterized by a slow, intermittent flow-like movement of prevalently plastic, clayey soil, facilitated by a combination of sliding along multiple discrete shear surfaces and internal shear strains. Following Hungr et al. [10], the landslide can be classified as earth flow. The Vomice earthflow (Figures 4c and 5a) is a large and complex landslide system with a maximum length of 1.85 × 10 3 m, a width ranging between 40 and 300 m, and an area of about 4.21 × 10 5 m 2 (~42 ha).
The landslide is characterized by a slow, intermittent flow-like movement of prevalently plastic, clayey soil, facilitated by a combination of sliding along multiple discrete shear surfaces and internal shear strains. Following Hungr et al. [10], the landslide can be classified as earth flow. The depletion zone consists of two strongly deformed source zones that feed the main earthflow, respectively named the southern and northern sources (Figure 5a), where a series of rotational or compound (roto-translational) slides mainly occur (Figure 5b,c). Source zones are characterized by active arc shaped main scarps (with an average height of some meters) overprinted by evident striae (Figure 5d), retrogressive evolution of the The depletion zone consists of two strongly deformed source zones that feed the main earthflow, respectively named the southern and northern sources (Figure 5a), where a series of rotational or compound (roto-translational) slides mainly occur (Figure 5b,c). Source zones are characterized by active arc shaped main scarps (with an average height of some meters) overprinted by evident striae (Figure 5d), retrogressive evolution of the crowns, a large occurrence of open fractures, trenches, upslope-tilted blocks, active erosion, and water-filled depressions; in addition, the morphology of the source areas is made Remote Sens. 2021, 13, 120 13 of 25 even more complicated by the presence of independent small and shallow earth-slides and earthslides-earthflow, which often occupy the landslide terraces or trigger at the foot of the main slide. The average terrain gradient of the source zones is about 20 • .
The transport zone of the landslide (where landslide style is supposed to change from pure sliding into a proper earthflow) consists of an elongated channel through which the material, produced in the depletion area, is transferred downward (with strong internal deformation)-through a run-out channelized flow-to the deposition zone (Figure 5e).
The channel is confined between marked incisions, the activity of which is testified by the presence of well-preserved striated shear planes (Figure 5f,g). Along this tract, the landslide is characterized by an evident hummocky morphology, produced by extensional cracks and contractional corrugation. In this portion, the terrain gradient may reach 15 • . The deposition zone (or accumulation area) (Figure 5h), in correspondence with the confluence with the Straface river valley bottom, consists of a fan-shaped toe that slightly deflects the watercourse. Here, diffuse cracking, terrain corrugation, and water ponds located within ephemeral topographic depressions are the most easily recognized features within the landslide body, while deep erosive channels, linked to the water's runoff, are found in the coincidence of the right and left side of the flow (Figure 5h,i). In this sector, the averages terrain gradient is around 5 • .
Data concerning the depth of the main sliding surfaces in the southern source area of the earthflow were obtained based on geomorphological field surveys and also thanks to the presence of a drainage channel that cut into bedrock. Particularly, the toe of the main sliding surface is found inside the contact zone between the chaotic and schistose facies of the Varicolored clays-shale Formation (Figure 6a), where slickensides were observed (Figure 6b). crowns, a large occurrence of open fractures, trenches, upslope-tilted blocks, active erosion, and water-filled depressions; in addition, the morphology of the source areas is made even more complicated by the presence of independent small and shallow earthslides and earthslides-earthflow, which often occupy the landslide terraces or trigger at the foot of the main slide. The average terrain gradient of the source zones is about 20°. The transport zone of the landslide (where landslide style is supposed to change from pure sliding into a proper earthflow) consists of an elongated channel through which the material, produced in the depletion area, is transferred downward (with strong internal deformation)-through a run-out channelized flow-to the deposition zone (Figure 5e).
The channel is confined between marked incisions, the activity of which is testified by the presence of well-preserved striated shear planes (Figure 5f,g). Along this tract, the landslide is characterized by an evident hummocky morphology, produced by extensional cracks and contractional corrugation. In this portion, the terrain gradient may reach 15°. The deposition zone (or accumulation area) (Figure 5h), in correspondence with the confluence with the Straface river valley bottom, consists of a fan-shaped toe that slightly deflects the watercourse. Here, diffuse cracking, terrain corrugation, and water ponds located within ephemeral topographic depressions are the most easily recognized features within the landslide body, while deep erosive channels, linked to the water's runoff, are found in the coincidence of the right and left side of the flow (Figure 5h,i). In this sector the averages terrain gradient is around 5°.
Data concerning the depth of the main sliding surfaces in the southern source area of the earthflow were obtained based on geomorphological field surveys and also thanks to the presence of a drainage channel that cut into bedrock. Particularly, the toe of the main sliding surface is found inside the contact zone between the chaotic and schistose facies of the Varicolored clays-shale Formation (Figure 6a), where slickensides were observed (Figure 6b).  In the transport and deposition zones, varicolored clays, stretched and deformed, were observed immediately below the earth-flow deposit (Figure 6c,d). The thickness of the displaced material ranges from about 20 to 25 m in the central body of the source zones (respectively northern and southern), decreasing to 5-10 m in the channel zone and increasing up to 15 m in the accumulation zone, near the landslide toe.
In addition, in the catchment area, clayey terrains exhibit high structural dynamism characterized by cracks, due to the shrinkage of clays, at the surface in the dry season, subsequently undergoing water infiltration, with consequent swelling of the clays in the wet season. This dynamism produces widespread phenomena of intense erosion, particularly along steep slopes and morphological scarps. Moreover, rills and ephemeral gullies, in some cases, affect landslide bodies and significantly contribute to sediment production.

Geomorphic Changes Detection of the Vomice Earthflow between 1954 and 2019
The analysis of ortho-images allowed the outlining of the spatio-temporal evolution of the Vomice earthflow from 1954 to 2019. The results, summarized in Figure 7, show that the landslide during this period (1954-2019) has undergone considerable modification. In the transport and deposition zones, varicolored clays, stretched and deformed, were observed immediately below the earth-flow deposit (Figure 6c,d). The thickness of the displaced material ranges from about 20 to 25 m in the central body of the source zones (respectively northern and southern), decreasing to 5-10 m in the channel zone and increasing up to 15 m in the accumulation zone, near the landslide toe.
In addition, in the catchment area, clayey terrains exhibit high structural dynamism characterized by cracks, due to the shrinkage of clays, at the surface in the dry season, subsequently undergoing water infiltration, with consequent swelling of the clays in the wet season. This dynamism produces widespread phenomena of intense erosion, particularly along steep slopes and morphological scarps. Moreover, rills and ephemeral gullies, in some cases, affect landslide bodies and significantly contribute to sediment production.

Geomorphic Changes Detection of the Vomice Earthflow between 1954 and 2019
The analysis of ortho-images allowed the outlining of the spatio-temporal evolution of the Vomice earthflow from 1954 to 2019. The results, summarized in Figure 7, show that the landslide during this period (1954-2019) has undergone considerable modification.  The earthflow has multiple sources that were most likely active at different periods. In particular, the visual interpretation of eight ortho-images highlighted the fact that the main scarps and the flanks of the source areas were affected by frequent failures, which caused a remarkable retract and an enlargement in the crown zones; while, in the accumulation zone, was observed both an enlargement and a progressive advancement of the earthflow body that involved the Vomice channel.
The retreat rates of the source areas varied significantly in space and time (Figure 8). For the time span analyzed, the maximum earthflow source area retreat (more than 200 m) occurred in the crown direction, corresponding to a mean annual retreat rate of approx. 3 m for the 65 years; however, sideward expansion was also noticed, which, at several points, reached 100 m. In addition, the advancement of the earth flow front, from 1954 to 2019, can be estimated at being in the order of 200-300 m (Figure 8). The earthflow has multiple sources that were most likely active at different periods. In particular, the visual interpretation of eight ortho-images highlighted the fact that the main scarps and the flanks of the source areas were affected by frequent failures, which caused a remarkable retract and an enlargement in the crown zones; while, in the accumulation zone, was observed both an enlargement and a progressive advancement of the earthflow body that involved the Vomice channel.
The retreat rates of the source areas varied significantly in space and time (Figure 8). For the time span analyzed, the maximum earthflow source area retreat (more than 200 m) occurred in the crown direction, corresponding to a mean annual retreat rate of approx. 3 m for the 65 years; however, sideward expansion was also noticed, which, at several points, reached 100 m. In addition, the advancement of the earth flow front, from 1954 to 2019, can be estimated at being in the order of 200-300 m (Figure 8).  (Table 4). This represents an annual extension rate of 2.78 × 10 3 m 2 /yr. During the same period, the results of the The comparison of the eight maps and the spatial analysis, performed in the GIS environment, showed that, from 1954 to 2019, the total area covered by the earthflow gradually increased from 2.41 × 10 5 m 2 in 1954 to 4.21 × 10 5 m 2 in 2019 (Table 4). This represents an annual extension rate of 2.78 × 10 3 m 2 /yr. During the same period, the results of the above analyses indicate a progressive increase of the earthflow source zones, ranging from 1.67 × 10 5 m 2 in 1954 to 2.78 × 10 5 m 2 in 2019, which represents an increase of more than 66% (Table 4). Moreover, the results revealed that the highest rate of earthflow extension (approximately 7.0 × 10 3 m 2 /yr) occurred from 2007 to 2012 (Table 4), which indicates a high activity of the earthflow during this period, both in the source areas and accumulation zone. The source areas of earthflow undergo a progressive upslope migration of the crowns because of continued mass movements, while the accumulation zone displays a toe advancement due to the fresh earth/debris flows accumulated in the lower part of the channel (Figure 7). Contrary to this, the comparison of the ortho-images of 2016 and 2019 showed the smallest areal changes, with a rate of planimetric variations of about 0.8 × 10 3 m 2 /yr (Table 4), where several rotational slide and earth flow phenomena affected the source zones, remobilizing old earthflow deposits (Figure 7). The various rates of planimetric variation of the earthflow may be related to the different intensities and durations of meteorological events that occurred in these two periods. The time span between 2007 and 2012 was characterized by two intense rainfall events that involved the whole Calabria region, respectively, in the winter of 2008/2009 and the winter of 2009/2010, causing several slope movements [68]. Whereas, in the period 2016-2019, the rainfall events were generally characterized by low intensity and short duration and, consequently, the slope failures were very low. The maps of the elevation differences (Figure 9), obtained by subtracting one DTM from the other, allowed us to quantify the surface changes and the displacement rates within the earthflow area (i.e., positive or negative values corresponding to gains and losses of material).
The results showed clearly that since 1954 to 2019 (the whole period considered), due to large mass movements and erosion processes, the earthflow dimension and morphology changed significantly both in terms of the depletion and the accumulation zones. Therefore, elevation changes exhibited high spatial variability between 1954 and 2019; the maximum variations in elevation connected to material lost and material gained within the earthflow area were −24 and +15 m, respectively (Figure 9). In addition, Table 5 summarizes the depletion and accumulation rates.
The DTM difference maps showed that the pattern of topographic changes within the earthflow area was very different in the three considered periods. A map of the elevation differences was created by subtraction of the 1954 and 2001 DTMs. The results evidenced significant lowering of the topographic surface in the source areas and in the upper portion of the deposit area of the earthflow (down to −23 m). A significant accumulation (up to 11 m) occurred in the middle sector of the earthflow body and in the ancient toe zone (Figure 9a and Table 5). Moreover, along the accumulation zone were observed local lowering to −5 m, which can be associated with the remobilization of old earthflow deposits.   Figure 9b) revealed that, in earthflow source areas, the elevation differences ranged from −0.01 m to −11 m, with a mean value of about −3 m; while in the the deposit area the mean value of elevation differences was +2.7 m (Table 5). Moreover, Figure 9b shows that the accumulation lobe was depleted by about −3 m in its upper sector. This depletion was compensated by an advancement of about 100 m of the landslide front (Figure 7), which caused an accumulation of about +10 m in the landslide foot zone (Table 5).
By exploiting the map of topographic surface changes, computed by means of comparison between the DTM of 2012 and the DTM of 2019 (Figure 9c), it was possible to reveal that, in the landslide source zones, the mean value of lowering of the topographic surface was −1.3 m (Table 5), with a local lowering exceeding −10 m in the crown area ( Figure 9c). Furthermore, in several areas were recorded elevation increases caused by mass movements that were triggered along the earthflow head scarps whose accumulations tend to fill up the internal depression of the source area. Conversely, in the deposit zone, the thickness of accumulated terrains reached about +8 m, with a mean value of elevation differences of about +1.4 m ( Table 5).
For the entire analyzed period (1954-2019), the map of vertical displacements showed that the upper part of the landslide system is characterized by widespread depletion. In particular, the mean value of depletion in the landslide source areas was about −6 m ( Table 5). The head area is affected by sub-vertical collapse caused by rotational sliding  The comparison of the DTMs of 2001 and 2012 ( Figure 9b) revealed that, in earthflow source areas, the elevation differences ranged from −0.01 m to −11 m, with a mean value of about −3 m; while in the the deposit area the mean value of elevation differences was +2.7 m (Table 5). Moreover, Figure 9b shows that the accumulation lobe was depleted by about −3 m in its upper sector. This depletion was compensated by an advancement of about 100 m of the landslide front (Figure 7), which caused an accumulation of about +10 m in the landslide foot zone (Table 5).
By exploiting the map of topographic surface changes, computed by means of comparison between the DTM of 2012 and the DTM of 2019 (Figure 9c), it was possible to reveal that, in the landslide source zones, the mean value of lowering of the topographic surface was −1.3 m (Table 5), with a local lowering exceeding −10 m in the crown area ( Figure 9c). Furthermore, in several areas were recorded elevation increases caused by mass movements that were triggered along the earthflow head scarps whose accumulations tend to fill up the internal depression of the source area. Conversely, in the deposit zone, the thickness of accumulated terrains reached about +8 m, with a mean value of elevation differences of about +1.4 m ( Table 5).
For the entire analyzed period (1954-2019), the map of vertical displacements showed that the upper part of the landslide system is characterized by widespread depletion. In particular, the mean value of depletion in the landslide source areas was about −6 m ( Table 5). The head area is affected by sub-vertical collapse caused by rotational sliding (maximum ground lowering is about −24 m) (Figure 9d). The transport of the collapsed material towards the lower parts of deposition zone caused a mean sediment accumulation of about +6.1 m ( Table 5). Deposition of material is most distinctive in the middle and lower parts of the Vomice channel and reaches maximum surface changes of +15 m (Figure 9d).
The topographic changes were also studied by means of the comparison of five representative cross sections, extracted, respectively, from 1954, 2001, 2012, and 2019 DTMs ( Figure 10). The analysis of these profiles confirmed that the topographic surface changed significantly between 1954 and 2019. The longitudinal profiles (a-a' and c-c') highlighted clearly that the most evident topographic changes appear upstream and downstream of the earthflow due the frequent occurrence of slope failures causing the retreat of the earthflow heads and sediment accumulation downslope. The two traverse profiles (b-b' and d-d'), tracked within the earthflow source area, showed negative topographic changes (depletion of terrain), whereas the profile e-e', located in the lower part of the accumulation zone, showed positive changes that prove the filling of the Vomice channel. Moreover, this profile highlights two lateral deep incisions along the accumulation zone due to gully erosion processes (Figure 5i). (maximum ground lowering is about −24 m) (Figure 9d). The transport of the collapsed material towards the lower parts of deposition zone caused a mean sediment accumulation of about +6.1 m ( Table 5). Deposition of material is most distinctive in the middle and lower parts of the Vomice channel and reaches maximum surface changes of +15 m (Figure 9d).
The topographic changes were also studied by means of the comparison of five representative cross sections, extracted, respectively, from 1954, 2001, 2012, and 2019 DTMs ( Figure 10). The analysis of these profiles confirmed that the topographic surface changed significantly between 1954 and 2019. The longitudinal profiles (a-a' and c-c') highlighted clearly that the most evident topographic changes appear upstream and downstream of the earthflow due the frequent occurrence of slope failures causing the retreat of the earthflow heads and sediment accumulation downslope. The two traverse profiles (b-b' and dd'), tracked within the earthflow source area, showed negative topographic changes (depletion of terrain), whereas the profile e-e', located in the lower part of the accumulation zone, showed positive changes that prove the filling of the Vomice channel. Moreover, this profile highlights two lateral deep incisions along the accumulation zone due to gully erosion processes (Figure 5i).   Table 6 summarizes the volumetric changes based on the vertical displacement computations in the different periods. The volume assessment, during the whole period analyzed (1954-2019), showed a depletion of more than 1.56 × 10 6 m 3 of material mobi-lized in the landslide source area, but the deposited mass in the accumulated zone is about 0.87 × 10 6 m 3 . The mass balance (the difference between depleted and accumulated volumes) shows a deficit of 0.69 × 10 6 m 3 , probably caused by intense gully erosion processes, which particularly act along the flanks of the deposit zone (see Figure 5h), and by river erosion that erodes the landslide toe.  (Table 6).

Discussion
The integration of multi-temporal ortho-images and DTMs obtained by different data sources is generally a useful tool for the characterization and evaluation of earthflow evolution [18,20,27,39,69]. In the literature, the comparisons of ortho-images and DTMs, acquired by diverse sensors and/or instruments at different times, were used to investigate landslide activity and volume estimation and their variations [8,37,[70][71][72]. Therefore, for the comparison of data acquired using different methods, with different precisions and resolutions, which were collected in various periods, the choice of spatial resolution for the data analysis and the subsequent error estimation becomes crucial [73]. In our case study, ortho-images and DTMs with variable spatial resolutions were generated, mainly using different methods of acquisition and various spatial scales of data. For example, the ground resolution of historical aerial photos, ranging from 1 m/pixel to 0.64 m/pixel, is not enough to be accurately co-registered with the recent UAV dataset (0.07 m/pixel of resolution). Nevertheless, when the entity of displacement of mass movements is particularly considerable and the morphological changes are highly visible, data with different ground resolutions can be employed in order to reconstruct the geomorphological evolution [8,40]. For the landslide investigated in this work, the entity of displacements of mass movements is of the same order of magnitude as the accuracy of historical aerial photos. For this reason, it was possible to make a direct comparison between archive image results and new UAV data.
The UAV data coupled to field observations were fundamental for defining the geological and structural setting of the study catchment and for recognizing the geomorphological features of the Vomice earthflow (Figures 4-6). The complex tectonic history experienced by Meso-Cenozoic sedimentary rocks [48] accounts for both the occurrence of a close net of discontinuities (faults and joints) and the chaotic arrangement of sediments. Additionally, the hydrographic network is strongly controlled by tectonic patterns, and there are several examples where stream flow direction has been controlled by structural features (Figure 2).
Particularly, the southern source area of the Vomice earthflow-where the depth of the main sliding surface is located inside the contact zone between the chaotic and schistose facies of the Varicolored clays-shale Formation (Figure 6a)-is controlled by various highangled fault segments belonging to the N-S and NE-SW systems, to which extensive fault zones are associated (Figures 2 and 3). The recognized faults greatly contributed to the weakening of the substrate (Figure 3b,c,e), presenting a preferential route for groundwater infiltration and migration.
The geomorphological study provided data on landslide deformation pattern and style; the movement is retrogressive on the source areas, advancing in the mid-lower sector and partially widening at the toe. The scale of deformations and displacement velocities decreases, moving from the source sector in the direction of the toe, which is generally the last section of the landslide to reactivate. Moreover, we can assume that, during the paroxysmal phases of the landslide, the usual displacement rate is "slow", accordingly to [13], but in limited sectors of the landslide and for a limited period of time this rate may be classified as "moderate" or "rapid". The response of the Vomice earthflow to rainfall is typically delayed, with long periods of accumulated precipitation required to cause the earthflow to move [20,44,71,82,83].
With regard to the geomorphic changes detection, both long-and short-term multitemporal analyses, performed by means of the archive aerial photographs, orthophotos, and Google Earth images, indicate that the Vomice earthflow, between 1954 and 2019, was periodically affected by slope failures, which caused a progressive expansion of both the source areas and the lower sector of the depositional area. The results highlighted the fact that, during the last 65 years, the landslide area has undergone deep morphological changes (Figures 7 and 8) and a significant increasing of the landslide surface (Table 4). Large terrain displacements were recorded in the upper parts of the source areas, characterized by multiple landslide scarps, as highlighted by the detailed geomorphological survey; in these areas, the thickness of collapsed material is more than of 20 m. Moreover, within the landslide source areas, where many single shallow landslides occur, counter slopes, trenches, and depressions of a small size, often masked by the mostly pelitic nature of lithology, were observed. In addition, several tension cracks affect the main scarps, which are subjected to significant retrogressive and enlarging evolutions.
Thanks to the overall analysis of the results, we discerned that the landslide activity is spatially and temporally complex. Since 1954, earthflow landslide alternates between long periods of relatively slow movement and moderately rapid accelerations. During the last 15 years, its activity increased considerably. The time span from 2007 to 2012 can be considered as one of the most intense active periods, causing an important extension of the source area, because of its retrogressive evolution. Consequently, a prominent advancement of the accumulation area and a wide fan-shaped toe was detected (Figure 7). The different rates of movement of the landslide can be related to the severity of rainfall events [84]. According to Gullà et al. [68] and Ferrari et al. [85], this sector of Calabria has been hit by a series of intense meteorological events during the last 65 years, which have caused the activation or reactivation of several landslides. The most severe rainfall events occurred in hydrological years 1953-1954, 1959-1960, 1972-1973, 1976-1977, 1989-1990, 2008-2009, and 2009-2010. Therefore, we can assume that the most important reactivations of the Vomice earthflow occurred during these rainfall events.
Moreover, the analyses of ortho-images and field observations showed that the landslide, in the upper portion of the slope, involved and damaged several sectors of a municipal road (Figure 8).
The multi-temporal analysis further revealed that, within the Vomice earthflow area, different slope sectors have being active at different times. This behavior is typical of earthflows and has been observed in similar various geological and physiographical settings [5,8,16,17,70,[86][87][88].
The comparison of the four different DTMs (Figure 9) demonstrated the complex distribution, both spatially and temporally, of the depletion and accumulation zones of the landslide. This quantitative assessment of elevation changes allowed us to determine the rates of surface modification for different landslide sectors; moreover, the analysis provided information on the volume of material mobilized from the source area, both that of accumulation and that of material lost due to erosion. From 1954 to 2019, the landslide morphodynamic changed considerably in the source zones, in the transition zone, and even in the accumulation zone, where large vertical displacements of the terrain surface were observed, both in terms of the depletion and the uplift zones ( Figure 9). The volumetric analyses clearly showed that different parts of the Vomice earthflow were active at different times, and with different rates of topographic change. During the entire period analyzed, the volume of the material depleted was about 1.56 × 10 6 m 3 , whereas the deposited material along the transition zone and in the accumulation zone was about of 0.87 × 10 6 m 3 . The mass balance shows a deficit of 0.69 × 10 6 m 3 (about 44% volume loss). Negative mass balances were generally attributed to the removal of material by surface runoff or river erosion, or, in some case, is the result of DTM errors [38,88,89]. In our case, the multi-temporal analysis and field observations showed that part of the mobilized material reached the Straface River and was eroded by the water flow, which has also limited the development of the accumulation area. The negative balance can be also attributed to the intense erosive processes, mainly rills and gullies, that act along the accumulation zone (Figure 5h,i). However, in agreement with [90], even the accuracy of DTMs may have influenced the elevation difference results.
Overall, the approach used in this study, based on the integration of multi-source and multi-temporal images and DTMs, coupled with field surveys, strongly improves the knowledge of the geological and geomorphological features of the study area and identifies the topographic changes that have affected the Vomice earthflow over the last 65 years. The data collected provide a basis for an interpretation of the long-term morphodynamic and kinematic evolution of the landslide. In addition, these data could give useful support for the monitoring, numerical modelling, and hazard mitigation of the landslide.

Conclusions
This paper emphasizes the importance and great potentiality of the integration of multi-source and multi-temporal data for the diachronic reconstruction of morphological changes that have occurred in the last 65 years due to the Vomice earthflow, a typical and representative landslide phenomenon of wide sectors of the Apennine chain of southern Italy.
The Vomice earthflow is a 1.85 × 10 3 m long earthflow, located in northern Calabria, that has affected an area of 4.21 × 10 5 m 2 (~42 ha). The landslide exhibits two source zones: a narrow and elongated transport zone and a lobate accumulation zone. The landslide movement is retrogressive on the source areas, where a series of rotational or composite (roto-translational) slides occur, advancing in the mid-lower sector and partially widening at the toe. The landslide is characterized by a slow, intermittent flow-like movement of prevalently plastic, clayey soil (varicolored clays). The scale of deformations and displacement velocities decreases, moving from the source sector in the direction of the toe, which is generally the last section of the landslide to reactivate. The response of the Vomice earthflow to rainfall is typically delayed, with long periods of accumulated precipitation required to cause the earthflow to move.
The results of the comparison data allowed us to evaluate the space and time evolution of earthflow from 1954 to 2019. Geomorphic changes in the landslide, quantified using DTM differencing, has detected complex patterns between the source areas and the accumulation zone during the considered time span. This analysis highlighted a widespread lowering of the topographic surface in the source area and a significant uplift at the landslide toe. In addition, the performed multi-temporal analysis showed a high increase of landslide surface (more than 66%) over the last 65 years and a parallel increase of the volume of mass moved.
However, a very important aspect of this quantitative analysis is the assessment of the data accuracy. The accuracy of the applied method and of the results obtained is highly dependent on the quality of the data used. For instance, the ground resolution of aerial photos and the accuracy of the instruments used for data acquisition of course influence the quality of the output data (e.g., orthophotos and DTMs). As our study demonstrated, plano-altimetric changes computed from multi-source DTM analysis may be used for areas affected by large landslides. In fact, in the study case, the values of errors of DTMs were acceptable compared to the topographic changes obtained for the considered time span.
In conclusion, the outputs of this work highlight the fact that the proposed approach can be considered a useful tool for investigating and monitoring the space-time morphological changes of similar landslides in other geo-environmental contexts.

Data Availability Statement:
The data presented in this study are available on request.