Mapping and evaluating kinematics and the stress and strain field at active faults and fissures: a comparison between field and drone data at the NE rift, Mt Etna (Italy)

We collected drone data to quantify the kinematics at extensional fractures and normal faults, integrated this information with seismological data to reconstruct the stress field, and critically compared the results with previous fieldwork to assess the best practice. As a key site, we analyzed a sector of the northeast rift of Mt Etna, an area affected by continuous ground deformation linked to gravity sliding of the volcano's eastern flank and dike injections. The studied sector is characterized also by the existence of eruptive craters and fissures and lava flows. This work shows that this rift segment is affected by a series of NNE- to NE-striking, parallel extensional fractures characterized by an opening mode along an average N105.7^∘ vector. The stress field is characterized by a σ_Hmin trending northwest–southeast. Normal faults strike parallel to the extensional fractures. The extensional strain obtained by cumulating the net offset at extensional fractures with the fault heave gives a stretching ratio of 1.003 in the northeastern part of the study area and 1.005 in the southwestern part. Given a maximum age of 1614 CE for the offset lavas, we obtained an extension rate of 1.9 cm yr⁻¹ for the last 406 years. This value is consistent with the slip along the Pernicana Fault system, confirming that the NE rift structures accommodate the sliding of the eastern flank of the volcano.

1 Introduction

Quantifying offset at recent and active faults and other structures is fundamental to deciphering the kinematics and stress and strain of deformation zones. This is a basic step for the assessment of seismic hazard (Lyakhovsky et al., 2012); in the case of volcanic zones it also contributes to the understanding of the crustal conditions that may facilitate magma ascent and thus the evaluation of volcanic hazard (Keir et al., 2006). At rift zones, the precise definition of the spreading direction and extensional rate requires the collection of a huge amount of data that must cover the whole rift extension. Since rift zones are composed of swarms of tens to hundreds of normal faults and extension fractures, the collection of a statistically robust amount of data requires a great fieldwork effort. Moreover, logistically complex conditions can affect the performance of fieldwork, as for example in the East Africa Rift system, where crustal extension rates have been frequently evaluated indirectly from plate tectonic models (Jestin et al., 1994; Chu and Gordon, 1999).

On active volcanoes, the presence of rough terrains and the possible exposition of researchers to explosive products frequently prevent optimal field surveys. These complex logistic conditions, in fact, do not permit having a detailed evaluation of the strain field due to the difficulties in obtaining a sufficiently large number of measurements along an extension fracture or a fault. In fact, only the collection of a large amount of horizontal dilation values can allow the precise reconstruction of the strain field. At faults, the reconstruction of heave and throw values requires the precise measurement of offset. In any case, the measurement of fault slip profiles is very time-consuming and can be very difficult in the case of faults with offsets in the order of tens of meters.

In the last few years, the abovementioned difficulties have been overcome by the use of structure-from-motion (SfM) photogrammetry applied to images collected by unmanned aerial vehicles (UAVs or drones) in active volcano-tectonics studies (Bonali et al., 2019a, 2020; Trippanera et al., 2019; Weismüller et al., 2019) and to assess volcanic hazard (Müller et al., 2017; Darmawan et al., 2018; De Beni et al., 2019). Therefore, in this work we use this technique, which allows us to reconstruct very detailed orthomosaics and digital surface models (DSMs) of the surveyed areas. The resulting images, which can attain a resolution as precise as 1 cm, allow us to collect several pieces of high-resolution structural data also in 3D and take direct measurements of structures and morphostructures, like dilation values along faults and fractures, even using immersive virtual reality tools (Tibaldi et al., 2020).

The present paper has a double focus: on the one hand, it describes new data useful for the interpretation of the activity of the NE rift, and on the other hand, it aims to present a methodology useful for similar studies. We show that the UAV-supported methodology can attain a precision comparable to field surveys in areas affected by active deformation. We also wish to show that UAV surveys have sufficient precision that may allow us to quantify the increment of extensional deformation by successive, repeated surveys. For this, we selected a sector of the NE rift, located on the northern summit part of Mt Etna (Italy) (Fig. 1), which is characterized by ongoing extensional fracturing, eruptive fissuring, and normal faulting. This also contributes to improving our knowledge of this important volcano-tectonic structure of Mt Etna, where only a few structural surveys were conducted several years ago by Garduño et al. (1997) and Tibaldi and Groppelli (2002). The 2002–2003 eruption took place here, accompanied by the development of new fractures and deposits, and thus new mapping is necessary.

Figure 1 (a) Map showing the geodynamic context where the volcano Etna is located. (b) Map showing the main structures of Mt Etna, with the eastern flank characterized by instability (a and b modified after Villani et al., 2020). (c) Digital elevation model of the upper part of Mt Etna with the location of the 118 well-constrained earthquake epicenters (yellow dots) that occurred during the 2002–2003 eruption (data credits: Tiziana Tuvè, INGV) and available focal mechanisms between 2008 and 2019 (blue dots) (source: http://sismoweb.ct.ingv.it/maps/eq_maps/focals/index.php, last access: 15 January 2021). The location errors of the dataset are on average 0.62 km for the epicenter and 0.41 km for the depth. The corresponding values for the median are 0.3 km (epicenter) and 0.1 km (depth). For more details on the location uncertainties, see Figs. 4 and 5 of Mostaccio et al. (2013). The black rectangle locates the area of the NE rift represented in Figs. 2 and 3. PF: Piano Provenzana Fault; PFS: Pernicana Fault system; NEC: Northeast Crater; VOR: Voragine Crater; BN: Bocca Nuova Crater; SEC: Southeast Crater.

The area turned out to be very suitable for such studies because (i) it is not covered by vegetation due to the high altitude (2000–2500 m a.s.l.); (ii) it is characterized by high deformation rates (in the order of 2 cm yr⁻¹; Tibaldi and Groppelli, 2002); (iii) the deposits affected by faulting and fracturing are historic, and as a consequence the effects of erosion are negligible and structures are perfectly preserved.

2 Geological background

Mt Etna, one of the most active volcanoes on Earth, is located in a compressional environment (Lanzafame et al., 1997; Cocina et al., 1997, 1998), on the border between the African and the European plates (Fig. 1a).

The volcano Etna has a constantly opened central conduit feeding four summit craters named Voragine (VOR; formed in 1945), the Northeast Crater (NEC; 1911), Bocca Nuova (BN; 1968), and the Southeast Crater (SEC; 1971) (Fig. 1c). After 2007, a new large summit cone grew on the southeast flank of the SEC: the New Southeast Crater (Del Negro et al., 2013; Behncke et al., 2014; Acocella et al., 2016). Eruptions from these summit craters are classified as summit eruptions (Acocella and Neri, 2003). Flank eruptions occur along radial fissures mostly on three main “rift zones”: the W rift, the S rift, and the NE rift (Fig. 1b) (Cappello et al., 2012). These flank (or lateral) eruptions are usually fed by shallow (1–3 km) dikes that propagate laterally from the central conduit (Acocella and Neri, 2009).

The NE rift is a network of N- to NE-striking eruptive fissures, 0.5 km wide and about 7 km long, extending from the NEC (∼ 3320 m a.s.l.) to ∼1400 m of altitude (Garduño et al., 1997) (Fig. 2). The upper portion of the rift strikes north from the summit down to 2500 m a.s.l., whereas the lower section strikes northeast down to the Monte Nero area (Fig. 1c). The rift is bordered to the southeast by a 200 m high tectonic scarp (Piano Provenzana Fault, PF in Figs. 1c and 2) partially covered by recent volcanic products. To the west, the rift is limited by a small scarp, crossed by recent eruptive fissures and largely concealed by historic lava flows and cinder cones. Faults and non-eruptive fractures cut the central portion of the NE rift; they strike between 0 and 60^∘ and have different kinematics, including pure extension or right-lateral and left-lateral transtension (Tibaldi and Groppelli, 2002). About 35 % of the fractures are associated with extrusive volcanic activity, which affected the lower portion of the same fractures: along these, hornitos, craters, and small and large cinder cones are common. Moving downslope, cones and craters take on a more and more pronounced elliptical shape, with the main axis striking 30–60^∘. Within them, the feeding magmatic dike often crops out, striking 10–20^∘ in the central part of the rift and 40–50^∘ in the northeastern part (Geshi and Neri, 2014).

Figure 2 Structural map of the NE rift superimposed on a shaded relief lidar-derived image of the area. Numbers show date of historical eruptive fissures (dotted red lines), modified after Neri et al. (2004). Yellow circles mark the epicentral location of earthquakes recorded during the 2002–2003 eruption (for location error, see caption of Fig. 1). In white, major faults; PF: Piano Provenzana Fault; PFS: Pernicana Fault system. The area outlined by the black line represents the area of Figs. 5–6 surveyed with the drones. The NE rift is located in Fig. 1.

At ∼1800 m a.s.l., the NE rift meets the 18 km long Pernicana Fault system (PFS, Figs. 1c and 2), an active left-lateral transtensional structure bounding the unstable flank of the volcano (Groppelli and Tibaldi, 1999; Acocella and Neri, 2005). Both the NE rift and the PFS are the NW margins of a wide sector of Etna involved in seaward displacement (Fig. 1b) (Borgia et al., 1992; Solaro et al., 2010; Ruch et al., 2010, 2013; Acocella et al., 2013; Apuani et al., 2013; Mattia et al., 2015), affecting an onshore area >700 km² (Neri et al., 2004) and with a thickness of 1–4 km (Ruch et al., 2010; Siniscalchi et al., 2012; Ruch et al., 2012). This corresponds to the unstable flank delimited by the upper slip surface of Guardo et al. (2020), since this upper surface tends to emerge in correspondence to the Etna summit–NE rift zone and goes from above sea level down to almost 4 km b.s.l., from west to east. The unstable area also continues below sea level, until it reaches the abyssal plain at a depth of over 2000 m (Urulaub et al., 2018). Several authors have recently highlighted the possible relationship between eruptive activity and flank deformation, showing that the acceleration of flank deformation may trigger flank eruptions and vice versa. In some cases, it was demonstrated that tectonic activity along the PFS triggers eruptions from the NE rift (Neri et al., 2004, 2005; Walter et al., 2005; Bonforte et al., 2011; Ruch et al., 2012; De Novellis et al., 2019).

Applying the lithostratigraphic units following the standards suggested by Salvador (1994), two main groups of volcanic deposits are detectable in the NE rift area (Fig. 3): the products belonging to the Il Piano Synthem (Mongibello Volcano lithosomatic unit; 15 420 ± 60–0 yr BP) and those belonging to the Concazze Synthem (Ellittico Volcano lithosomatic unit; 56.6 ± 15.4 ka–15 420 ± 60 yr BP) (Coltelli et al., 1994; Garduño at al., 1997; Coltelli et al., 2000). Both represent volcanic units made up of products erupted during the last ∼ 57 ka and belonging to the Mongibello Supersynthem (Branca et al., 2011).

Figure 3 Geological map of the NE rift showing lithostratigraphic units, the 2002–2003 earthquake epicenters (yellow dots), and the main faults (in white), superimposed on a shaded relief DSM derived from the area (Gwinner et al., 2006). The NE rift is located in Fig. 1.

With reference to the stratigraphy in the bottom part of Fig. 3, the Concazze Synthem coincides with Ellittico Volcano, a large stratovolcano with a main, summit eruptive vent approximately coinciding with the current summit of Etna but higher (3600–3800 m a.s.l.). The stratigraphic succession consists of alternating lavas and pyroclastic deposits. It ends with plagioclase-rich porphyritic lava flows and reddish subaphyric lavas and scorias (Pizzi Deneri Formation and Portella Giumenta Formation, respectively; Coltelli et al., 1994; Branca et al., 2011). The deposits of the final explosive activity at Ellittico date to 15 420 ± 60 yr BP (Condomines et al., 1982; Cortesi et al., 1988; Gillot et al., 1994; Coltelli et al., 2000; De Beni et al., 2011); this explosive activity also generated the formation of a large and deep summit caldera, whose remains today crop out at the edge of the highest portion of the NE rift, namely at Punta Lucia and Pizzi Deneri. In the NE rift area, the products of the Portella Giumenta Formation overlap deeply eroded cinder cones and porphyritic lavas belonging to Piano Provenzana Formation.

The Il Piano Synthem constitutes the present active volcano. The lower boundary coincides with the Ellittico caldera, while the upper boundary is the current topographic surface. In the NE rift, the volcanics belonging to the Pietracannone and Torre del Filosofo formations largely crop out. During the last 120 years, the NE rift eruptions lasted 21 d on average, with 7 m³ s⁻¹ eruption rates. In ∼ 1 d, the eruptive fissures reached the maximum length (3825 m) by propagating at an average speed of 0.053 m s⁻¹ (Neri et al., 2011).

3 Seismological data

The number of stations and geometry of the seismic network for the permanent monitoring of Etna has significantly changed in the last 40 years. The boost to the seismic network (from five short-period, vertical-component stations of the 1970s to the 30 broadband, three-component sensors of 2020) has also involved signal transmission (from analogue to digital) and acquisition systems (Patanè et al., 2004). In the framework of seismic activity at Etna, the NE rift, along with the PFS and the PF, plays a key role. Indeed, they form a discontinuity that accommodates stress changes related to the magma intrusion and tectonic loading, interpreted as the main sources causing the eastward sliding of the eastern sector of the volcano (Alparone et al., 2013a, b).

The area of interest in this case study is situated at the westernmost part of the PFS–PF fault system. The structure has a bow-like geometry, striking northeast at its western tip, and bending along an E–W strike direction towards east. Here the earthquakes are usually shallow (depths mostly between 0 and 3 km b.s.l.) and with small to moderate magnitude (< M_L 5) (Fig. 1c). Despite their magnitude, these superficial earthquakes can be damaging as documented by macroseismic studies, which highlight the high seismic hazard of this sector of the volcano (Azzaro et al., 1998; Azzaro, 2004). Alparone et al. (2013a, b) report that seismic activity at the PFS–PF fault system increased from September 2002 on, starting shortly (a month) before the 2002–2003 eruption. Overall, 874 earthquakes with M_max 4.1 heralded the onset of that eruptive episode. Focusing on seismic activity during the 2002–2003 eruption, Mostaccio et al. (2013) tested NonLinLoc (Lomax et al., 2000), a nonlinear probabilistic earthquake location method, using a 3D velocity model. From the 328 well-constrained locations obtained by Mostaccio et al. (2013), we extracted a subset of 118 shallow earthquakes, which are located in an area encompassing the zone of our case study, part of the Valle del Bove and of the PFS–PF fault system. We used this 2002–2003 dataset because it contains the best-located earthquakes and because it is representative of the typical seismic activity of this sector of Mt Etna. Figure 1c highlights a bow-shaped distribution of 90 epicenters starting from the summit caters. The striking correlation between epicenter location and structural elements is visible comparing Figs. 1 and 2, since both have an arcuate shape and seismicity recalls the bow shape of the faults and fracture distribution. Along the NE rift, the earthquake distribution trends north-northeast–south–southwest and tends to bend to northeast–southwest and east–west moving in an eastern direction. This group of earthquakes is well separated from a second, smaller group (26) located more to the south, whose position clearly marks the northern rim of the Valle del Bove (Fig. 1). It is worth noting that the peculiar distribution of earthquakes in Fig. 1 is not only typical of the 2002–2003 eruptive period. Indeed, analyzing the distribution of seismic foci during the years from 2000 to 2009 (with the exclusion of time spans with volcanic activity), Alparone et al. (2013a, b) identified a cluster of earthquakes, which closely marked the NE rift and which are aligned along the PFS–PF fault system. The magnitude of these earthquakes was small, with only two of them with M_L>2. A second cluster belonged to the easternmost sector of the PFS–PF; they had stronger magnitude, reaching values M_L 4.1.

Figure 1c depicts the fault plane solution calculated for the earthquakes reported in Table S1 (Supplement), which occurred from 2008 to 2019 (source: http://sismoweb.ct.ingv.it/maps/eq_maps/focals/index.php, last access: 15 January 2021). Even though the mechanisms show some scatter, common elements may be identified. The earthquakes located in the southwestern corner of the picture have T axes striking mostly north, northeast, and east–northeast. Most of these earthquakes have strike-slip or normal faulting mechanisms, with almost vertical P axes. The earthquakes located in the northern part of the figure express the general trend of deformation along the PFS–PF fault system, with either normal faulting or a horizontal strike-slip mechanism. All mechanisms of this group have T axes striking in a SE direction.

4 Methods

4.1 Drone survey, photo, and ground control point (GCP) collection

In the present work, we applied the structure-from-motion (SfM) photogrammetry techniques, through drone surveys, to build up the high-resolution digital surface model (DSM) and orthomosaic for the target area. We followed the overall workflow that has been successfully tested in volcanic terrains and in challenging logistic conditions (Bonali et al., 2019a, 2020). Such a workflow has been designed to work with commercial quadcopters over large areas in volcanic terrain; that is exactly the situation we tackled in the present work where we used the DJI Phantom 4 PRO. This device is supplied with an incorporated chip set to work with a space-based satellite navigation/referencing system (GPS/GLONASS) and a high-resolution camera sensor (20 megapixels) in order to enhance the quality of the surveying as well as to obtain georeferenced pictures (geographic coordinates/WGS84). The overall area has been surveyed by several different flight missions where each of them has been planned to consider the presence of natural obstacles – mainly identified as topographic highs – known from a previous field survey, topographic maps, and satellite images. Flight height has been set up to 80–95 m above the ground, reaching the excellent pixel size of 2–3 cm for the resulting orthomosaic. We set the flight path considering the wind speed and direction and chose an overlap ratio of 85 % and 80 %, along the flight path and in a lateral direction, respectively (Gerloni et al., 2018; Antoniou et al., 2019; Bonali et al., 2019a, 2020; Fallati et al., 2020). The constant speed velocity was set up by the app considering all the above settings, and pictures have been captured using equal time interval modality.

As a parallel and complementary activity to the UAV survey, we collected several GCPs, distributed all over the area, essential to scale and reference the SfM-derived models as well as to avoid any bulging effect (James and Robson, 2012; Turner et al., 2012; Westoby et al., 2012; Smith et al., 2016; Vollgger and Cruden, 2016; James et al., 2017; Esposito et al., 2017). We targeted 34 natural targets, as already successfully performed by Bonali et al. (2020), to speed up the GCP collection avoiding the deployment and recovery of artificial targets. This method allowed us to save a day (8 h) of fieldwork. All GCPs have been collected with the GPS/GNSS Stonex S850A multi-frequency receiver in RTK configuration (with sub-centimeter accuracy). Depending on the 3G network availability, the GPS was linked in real time with the sicili@net network, a real-time correction service based on the caster NTRIP tool (http://193.206.223.39:5099/spiderweb/frmIndex.aspx, last access: 10 January 2021), or, in base-rover configuration, data were post-processed thanks to the Stonex Cube manager using Monte Conca and Pizzi Deneri permanent station correction data of the Istituto Nazionale di Geofisica e Vulcanologia (INGV) network.

4.2 Photogrammetry processing

For the photogrammetry processing of the 4018 captured pictures, we used a commercial SfM software – Agisoft METASHAPE (http://www.agisoft.com/, last access: 20 January 2021), which is commonly used due to the outstanding quality of the resulting output models (Cook, 2017; Burns and Delparte, 2017; Benassi et al., 2017). We also used the Agisoft Cloud beta service for data processing. The processing passed through some key steps (workflow), described in detail by Verhoeven (2011) and Brunier et al. (2016), which led to the realization of the DSM and the orthomosaic as final products (see Fig. 4 and Tables 1–2), which are based on the SfM-derived sparse and dense clouds.

Figure 4 Resulting DSM (a) and orthomosaic (b) of the surveyed area; blue dots represent the location of surveyed GCPs. Ref. system: UTM33/WGS84.

Table 1 Summary of setting and results related to the photogrammetry processing.

Table 2 Outline of time production for DSM and orthomosaic, including the time (in hh:mm:ss) for UAV survey and image collection.

We divided the overall workflow in four principal steps: (i) aligning of pictures, (ii) georeferencing, (iii) dense cloud generation, and (iv) DSM and orthomosaic production. A summary of all details regarding the photogrammetry processing and results are reported in Tables 1 and 2.

4.3 Mapping and measurement gathering on SfM-derived models

The DSM and the orthomosaic have been uploaded in a GIS environment, where we were able to trace all normal faults, extension fractures, and eruptive fissures we recognized in the area (Fig. 5a). Structures with a continuous vertical offset >20 cm, as measured on the DSM, have been classified as normal faults, whereas structures with a vertical offset <20 cm have been classified as extension fractures. Regarding eruptive fissures, they have been traced considering morphometric parameters of the eruptive centers, as explained in Tibaldi (1995), Bonali et al. (2011), and Tibaldi and Bonali (2017): the strike of the feeding fracture is directly related to the elongation of the cone base and the crater, to the direction of the line connecting the depressions on the crater rim, and to the alignment of cones. In some eruptive centers, the outcropping dike was visible, giving information about the direction and component of opening of the eruptive fissure.

Figure 5 (a) Geological–structural map of the studied area; structures are traced on high-resolution SfM-derived models; reference system: UTM33N-WGS84. Location in Fig. 2. Rose diagrams show the strike of the extension fractures (b), normal faults (c), and eruptive fissures (d).

On the basis of the derived models, we collected a total of 574 structural data (432 at extension fractures and 142 at faults): on the DSM, we measured strike and vertical offset along the normal faults by calculating the difference in elevation along topographic profiles traced every 10 m, orthogonally to the fault scarp. With regard to extension fractures, we determined the local strike, the opening direction vector, and the amount of net dilation. The latter two values were obtained by tracing a line that connects the two piercing points located on the opposite sides of an extension fracture, whenever they were undoubtedly recognizable. The length of the line gives the net dilation and related opening vector.

Finally, along two NW–SE-trending transects traced parallel to the resulting overall direction of opening, the total amount of horizontal dilation has been calculated, cumulating each single value, as well as the extension rate and stretch in the area. These values were compared with the data collected in the field by Tibaldi and Groppelli (2002).

5 Results

We surveyed an area of 2.2 km² through the collection of a total of 4018 photos. Thanks to the above-described workflow, a high-quality orthomosaic and a DSM were reconstructed (Fig. 4), with a resulting ground resolution of 2.97 and 11.86 cm pixel⁻¹, respectively. In these models, we recognized the presence of 20 normal fault segments, 250 extension fractures, and 54 eruptive fissures (Fig. 5a). The extension fractures strike mainly N20–50^∘, as highlighted in Fig. 5b, with an average strike of N24.6^∘. Normal faults strike mainly N10–40^∘ (Fig. 5c) with an average strike of N25.3^∘, whereas eruptive fissures strike mainly N20–40^∘ (Fig. 5d) with an average strike of N29.8^∘. Regarding normal faults, they depict a graben running at the foothill of the westernmost pyroclastic cone, known as Monte Pizzillo (Fig. 5a). These faults have a maximum height of the SE-dipping scarp of 2.3 m (Fig. 6a), giving a heave of 0.6 m and assuming a fault plane dip of 75^∘. The fault scarp facing northwest reaches a height of 3 m giving a heave of 0.8 m. The fault scarp located towards the southeastern part of the studied area faces southeast and reaches a maximum height of 4.2 m, resulting in a heave of 1.1 m.

Figure 6 Geological–structural map showing (a) vertical offset amounts along normal faults and values of horizontal dilation measured along the two transects (black lines: traces of the transects) and (b) amount of opening and opening direction in the 144 detected structural stations. Rose diagram shows the distribution of net opening direction (grey) and of the respective local fracture azimuth (blue).

To better understanding the active deformation processes affecting the area, we also collected a series of quantitative data at 144 sites along the extension fractures, totalling 432 structural measurements (Fig. 6). The latter include (i) local fracture strike; (ii) the amount of fracture dilation; and (iii) the opening direction. Opening direction values are in the range N72–163^∘, with a mean value of N105.7^∘, and most values are between N90–100^∘ (rose diagram in Figs. 6b and 7a). The fracture strike compared with the fracture opening direction highlights a clockwise rotation of fracture strike with the increase in opening directions (Fig. 7a). We have also quantified the local extension fracture azimuth, obtaining values between N329.8^∘ W and N78.8^∘, with a peak between N0–10^∘ and a mean value of N19.1^∘ (rose diagram in Fig. 6b), suggesting a slight overall left-lateral component of 3.4^∘. In more detail, the fractures with a lateral component <5^∘ were here classified as pure extensional fractures, whereas the remaining fractures have a left-lateral or a right-lateral component (Fig. 7b), counting 65 fractures with a left-lateral component, 40 fractures with a right-lateral component, and 39 pure extensional fractures, out of our 144 total data. Moreover, the component of left-lateral motions (up to 52^∘) is larger than the right-lateral component (up to 36^∘). The graph of Fig. 7b also shows the relation between lateral components of motions and fracture azimuth: with an increase in the fracture strike, the lateral component tends to change from the right-lateral component to the left-lateral component.

Figure 7 (a) Graph comparing the local strike of extension fractures with the net opening direction at each site. (b) Graph comparing the local strike of extension fractures with the lateral component of motions at each site. (c) Histogram showing the frequency of the net dilation amount values measured along extension fractures. (d) Histogram showing the frequency of vertical offset values measured along normal faults. (e) Graph comparing the local strike of extension fractures with the net dilation amount at each site. (f) Graph comparing the strike of all the 250 extension fractures with their length. (g) Graph comparing the length of each structure with its maximum offset/dilation values, distinguishing between SE-dipping faults, NW-dipping faults, and extension fractures.

The dilation values measured along extension fractures are in the range of 0.07–4.14 m (Figs. 6b and 7c), the average value is 0.4 m, and almost all values <1 m. These values have been related to the local strike, showing that the greatest dilation values are associated with strike values of about N20^∘E; moving away from these strike values, dilation decreases gradually, especially if strike rotates in a NW–SE direction (Fig. 7e). Regarding structure length, extension fractures reach a maximum length of 93.8 m, with an average value of 13.8 m. Fractures with greater lengths show strike values between N40–60^∘ E, decreasing gradually if strike rotates in a counterclockwise direction and more abruptly if it rotates in a clockwise direction, as shown in Fig. 7f.

Regarding normal faults, we measured vertical offset every 10 m along all the fault segments, obtaining values included in the range 0.1–7 m (Figs. 6a and 7d); the average value is 1.6 m and about half of the values <1.5 m. Faults present greater lengths, reaching a maximum of 299.6 m, with an average value of 92.2 m (Fig. 7g). In particular, in Fig. 7g it is evident that greater lengths characterize SE-dipping faults, with respect to NW-dipping ones, which are all <85 m. Also, regarding the amount of displacement, we can observe that SE-dipping faults are characterized by greater values of offset than NW-dipping ones: the former reach a maximum offset of 7 m, whereas the latter reach a maximum value of 3 m (Fig. 7g). Regarding length / displacement ratios for normal faults, these are between 11.3 and 284.7, with an average value of 67.9.

Considering a dip of 75^∘ for normal faults, we were able to calculate the dilatational component at both extension fractures and normal faults. Thus, we determined the total extensional component along two transects (traces in Fig. 6a), both with a length of 1.43 km, oriented in the given overall spreading direction, obtaining a total value of 5.3 m in the northern part of the area and 7.6 m in the southern part, which correspond to a stretching ratio of 1.003 in the north and 1.005 in the south.

6 Discussion

6.1 Rift geometry, structuring, and kinematics

The part of the rift where we focused our study belongs to the 4.5 km long NE rift of Mt Etna, a volcano-tectonic feature that showed important volcanic and tectonic activity in historic times. The orientation of the structures surveyed by UAVs is coherent with the remaining structures of the central part of the NE rift, which generally strike northeast–southwest. The remaining northeastern and southwestern portions of the NE rift show a slight clockwise and counterclockwise rotation with respect to the central part, giving the rift a gentle concavity towards the southeast (Fig. 2). Most of the rift, and similarly the studied area, is dominated by extension fractures and a few normal faults. Among faults, the largest one is represented by the PF normal fault, facing southeast, which constitutes the westernmost termination of the PFS (Groppelli and Tibaldi, 1999; Tibaldi and Groppelli, 2002; Acocella and Neri, 2005). Although the PFS has left-lateral strike-slip to transtensional motions, as shown by the focal mechanism solutions of Fig. 1, in the studied sector the PF produces an escarpment, facing southeast, about 200 m high that separates the ridge of the NE rift from the flat plain of Piano Provenzana. This dominant downdip motion is linked to the rotation of the PFS–PF fault system that turns from an E–W orientation in the eastern part towards a NE–SW strike direction in the studied area. The other faults form small grabens, one of which is present in the studied sector: this graben is 35 m wide at the southwestern foothill of the Mt Pizzillo pyroclastic cone and widens up to 80 m in correspondence to the upper portions of the cone. We suggest that this geometry is mainly due to the interference between the fault dip and the shape of the conical edifice, as observed also at Mount Laki in Iceland (Trippanera et al., 2015) or along the Harrat Lunayyir Fault in Saudi Arabia (Trippanera et al., 2019).

Regarding normal faults in the area, the calculated length / displacement ratios (11.3–284.7) are smaller than the ones obtained by Gudmundsson et al. (2013), which are in the range 42–362 (average about 130). In any case, our work considered a much smaller area than the one studied by these authors, which studied all the principal faults of the eastern flank of the volcano, with lengths up to 12 950 m and displacements up to 190 m.

The 144 opening directions measured along the extension fractures on the images obtained by drone surveys and SfM, indicate a clear homogenous dominant extension vector trending N105.7^∘. This vector is perpendicular to the largest slope of the area that is represented by the scarp of the PF and coincides with the direction of the maximum slope gradient of this part of Mt Etna (Favalli et al., 1999). As a consequence, we retain that the opening vector of the NE rift is strongly influenced by gravity effects linked to the shape of the Mt Etna edifice. A comparison with seismicity indicates that here the focal mechanism solutions have T axes trending northwest–southeast (Fig. 1), consistent with the opening directions measured by the drone surveys.

The fact that the faults dipping to the southeast are longer, have larger offsets, and are more frequent than those dipping to the northwest, can be linked to the fact that the NE rift has a strongly asymmetric profile measured in a NW–SE direction. Towards the northwest, in fact, there is a smooth topographic decrease, whereas in the opposite direction there is the steep slope gradient produced by the PF scarp. Moreover, the northwestern side of the NE rift is buttressed by the stable northern volcano slope, whereas the southeastern side is involved in the gravity sliding of the eastern volcano slope. This means that there is a larger gravity force acting in the SE direction.

The fact that most opening took place at the fractures with a strike ranging N0–50^∘ and that the longest fractures strike N10–55^∘ can be explained by assuming that this range corresponds to the local orientation of the opening vector linked with the instability of the eastern volcano flank. Finally, the generally larger values of the length of faults with respect to fractures is compatible with the concept that extensional fractures represent an immature stage of the evolution towards faults (Gudmundsson, 1987; Acocella et al., 2003; Tibaldi et al., 2019); faults can, in fact, derive from the linkage between different fracture segments during the evolutionary process.

Notwithstanding the general orthogonality of the opening direction with respect to the average trend of the NE rift, in greater detail we observed that 45 % of all extension fractures present a left-lateral component, whereas 27 % are characterized by pure extension. Considering also that the component of left-lateral motions is larger than the right-lateral component, we conclude that the NE rift has a left-lateral transtensional behavior linked to the large strike-slip component of motions along the PFS.

Finally, we used the software “Lissage” (Lee and Angelier, 1994) and the unpublished software ATMO-Stress, prepared in the framework of the NEANIAS project (https://www.neanias.eu/, last access: 10 February 2021) of the EU to calculate the stress field. Lissage is a C-based software designed to reconstruct paleo-stress trajectories in a given area (Lee and Angelier, 1994), using multiple local stress determinations as input data, including P and T axes derived from seismological data, the direction of principal stress axes from stress inversion, and any other data that describe the azimuth of σ_Hmax or σ_Hmin. Such software can be used to reconstruct stress trajectories both using local field data and regional paleo-stress database (e.g., Hu et al., 1996; Munoz-Martín et al., 1998; Maestro et al., 2007; Bonali et al., 2019b). The Lissage software was here used to quantify the stress field at a broader scale by the interpolation of σ_Hmin resulting from T axes of focal mechanism solutions and from the single points of measurement of net dilation direction at extension fractures in our study area (Fig. 8a), assuming that net dilation is parallel to σ_Hmin. The ATMO-Stress software is the online version of such software and was here used to calculate in detail, at a more local scale, the stress field based only on net dilation direction at extension fractures (Fig. 8b). Both results indicate a clear NW–SE trend of σ_Hmin, although in detail there appears a slight counterclockwise rotation from east to west. We maintain that this rotation is linked to the transition from the strike-slip dominion of the PFS, which strikes east–west, to the more extensional dominion of the rift-PF, which strike northeast–southwest.

Figure 8 Stress field trajectories obtained by the interpolation of σ_Hmin in the whole NE rift (a) and in our study area (b). For stress computation, we used the program “Lissage” (Lee and Angelier, 1994) in (a), and the unpublished program ATMO-Stress in (b).

6.2 Extensional rate

The total extension measured along the faults and fractures that crop out in the studied area along the two transects is 5.3 m in the northern part of the area and 7.6 m in the southern part. Assuming the age of 1614–1624 CE for the oldest lavas affected by the brittle structures, we obtain an extension rate of 1.87 cm yr⁻¹ at the southern transect for the last 406 years. Our measured extension rates are based on a wider dataset than that published in the Tibaldi and Groppelli (2002) paper, who indicated extension rates measured at single fractures with values of 1.8, 1.3, and 0.6 cm yr⁻¹ along the rift moving from southwest to northeast. Our studied area corresponds to the southwestern portion of the NE rift that, based on the data of Tibaldi and Groppelli (2002), opens at a higher rate. Our data thus indicate that this part of the rift opens at a slightly higher rate than previously suggested (1.87 cm yr⁻¹ instead of 1.3–1.8 cm yr⁻¹) and we retain our result is based on a more robust statistical dataset.

This result is also coherent with the gross general rate of deformation measured along the nearest main structure that is represented by the PFS. This fault is linked to the NE rift, and this justifies its comparison with our data. The PFS slip-rate has been quantified in its central-western part at 0.4 ± 0.1 to 2.2 ± 0.1 cm yr⁻¹ and in the eastern part at 0.2 ± 0.1 to 0.8 ± 0.4 cm yr⁻¹, where the PFS splays in the Vena–Presa zone (Tibaldi and Groppelli, 2002). The PFS short-term slip-rate has been assessed to be 1 cm yr⁻¹ in the vertical component and 2.8 cm yr⁻¹ in the left-lateral component by Azzaro et al. (2001). A lower fault slip-rate of 0.6–1.5 cm yr⁻¹ since 3 ka has been more recently calculated by D'Amato et al. (2017), although this refers only to the throw rate and does not consider the strike-slip component that is especially high in the western and eastern part of the PFS. This consistency between the extension rate at the NE rift and the slip along the PFS confirms that they accommodate the seaward sliding of the eastern volcano flank.

6.3 Methodological aspects

Results from the present work support the utility of using UAV-based SfM as a complementary tool to increase quality data collection, in addition to classical fieldwork, here aimed at defining the architecture and active processes working in a rift zone in volcanic areas. If we consider the time necessary to carry out the drone surveys, plus the time necessary to process the data and interpret them, we reach a total of 10 d of work. This yielded the collection of a huge amount of structural data, comprising 432 structural measurements (opening directions, amount of opening, and local azimuth) at extension fractures. In comparison, Tibaldi and Groppelli (2002) collected just 22 opening directions in the field for the same area and the related structural map over a total of 1 week of work.

The use of UAVs in geoscience is quickly increasing for a series of reasons: (i) UAV data acquisition is cheaper with respect to other methods, such as airborne and terrestrial laser scanning and lidar (Cawood et al., 2017; Lizarazo et al., 2017); (ii) the proposed approach reduces work time compared to field data collection, especially in the case of the study of long structures as those of Tibaldi and Ferrari (1992), Kozhurin et al. (2006), and Trippanera et al. (2019); (iii) UAVs provide the possibility of reaching sites that can be inaccessible for logistic conditions or can be dangerous, such as an active volcano. Moreover, the excellent accuracy of the SfM-derived models allows us to carry out observations and measures at a level of detail in the order of centimeters and thus a scale that is comparable to field surveys. This scale, together with the velocity of the workflow and related processing, can even improve performance.

Furthermore, the rapid development of UAV technology (furnished with a real-time kinematic positioning system) will increase the accuracy of the SfM outputs in the future, flying also at a higher elevation from the ground. Similarly, the UAV flight stability, camera quality and battery capacity can also improve thanks to the continuous delivery of new UAVs and related items.

With regard to limitations arising from the use of the UAV-based SfM, we mention the dependence on the flight time that is linked to the battery life, weather conditions, flight rules, and the fundamental step of placing and surveying the GCPs. The latter are needed to precisely scale and reference the model, but their placement is time-consuming, slowing down the total time needed for the UAV missions, particularly over broad areas. Smith et al. (2016) suggested considering a minimum of three GCPs, whereas other authors considered a higher number of GCPs but decreased the overlap ratio among the pictures to 60 %–70 % (Javernick et al., 2014). Finally, by increasing overlap ratio, it is possible to decrease the number of GCPs considered for scaling and referencing (Esposito et al., 2017).

7 Conclusions

We applied UAV-based structure-from-motion (SfM) photogrammetry to analyze a high-altitude area characterized by rough terrains in the northeastern part of Mt Etna. The area is affected by the presence of the NE rift, a volcano-tectonic feature composed of NE-striking historic eruptive fissures, extension fractures, and normal faults. The stratigraphic deposits span an age from 56.6 ± 15.4 ka to nowadays, whereas in the area surveyed by the drones, the deposits are mostly historic in age.

The highly detailed drone survey, in the order of 2.8 cm of resolution, showed that the studied sector of the NE rift is affected by 250 extension fractures, 20 normal fault segments, and 54 eruptive fissures. The extension fractures strike mainly N20–50^∘, with an average strike of N24.6^∘. The normal faults strike mainly N10–40^∘ with an average strike of N25.3^∘. The eruptive fissures strike mainly N20–40^∘ with an average strike of N29.8^∘.

We measured 432 structural data, comprising local fracture strike, dilation amount, and opening direction. Opening direction values are in the range N72–163^∘, with a mean value of N105.7^∘. A comparison with respect to fracture strike indicates the presence of a slight overall left-lateral component of 3.4^∘. Moreover, 65 fractures have a left-lateral component, 40 fractures a right-lateral component, and 39 pure extension. The component of left-lateral motions (up to 52^∘) is larger than the right-lateral component (up to 36^∘). These data suggest the effect of the left-lateral normal Pernicana Fault system on the NE rift kinematics.

Seismicity seems to be concentrated along the Pernicana–Piano Provenzana faults, whereas earthquake focal mechanism solutions show T axes trending perpendicularly to the NE rift, consistent with the results of our measurements of dilation orientations.

Assuming the age of 1614–1624 CE for the oldest lavas affected by the studied fractures, we obtain an extension rate of 1.87 cm yr⁻¹ for the last 406 years. This rate is consistent with the data measured by other authors at the Pernicana Fault system.

Results indicate that SfM photogrammetry coupled with drone surveys allows us to collect large datasets with a detail comparable to field surveys. Drone survey has the advantage of collecting more data in the same time period with respect to classical fieldwork and also allows data collection in difficult terrains where logistics can represent an insurmountable obstacle.