Reconstruction of Late Pleistocene-Holocene Deformation through Massive Data Collection at Kraﬂa Rift (NE Iceland) Owing to Drone-Based Structure-from-Motion Photogrammetry

: In the present work, we demonstrate how drone surveys coupled with structure-from-motion (SfM) photogrammetry can help to collect huge amounts of very detailed data even in rough terrains where logistics can a ﬀ ect classical ﬁeld surveys. The area of study is located in the NW part of the Kraﬂa Fissure Swarm (NE Iceland), a volcanotectonic rift composed of eruptive centres, extension fractures, and normal faults. The surveyed sector is characterized by the presence of a hyaloclastite ridge composed of deposits dated, on a stratigraphic basis, to the Weichselian High Glacial (29.1–12.1 ka BP), and a series of lava ﬂows mostly dating back to 11–12 ka BP. The integration of remotely sensed surveys and ﬁeld inspections enabled us to recognize that this segment of the Kraﬂa rift is made of grabens arranged en- é chelon with a left-stepping geometry. A major graben increases in width in correspondence of the hyaloclastite cone; we interpret this geometry as resulting from the mechanical contrast between the sti ﬀ er lava succession and the softer hyaloclastites, which favours the development of concentric faults. We also measured a total extension of 16.6 m and 11.2 m along the fractures a ﬀ ecting the lava units, and a total extension in the hyaloclastites of 29.3 m. This produces an extension rate of 1.4 mm / yr in the Holocene lavas and 1.7 ± 0.7 mm / yr in the Weichselian hyaloclastite deposits. The spreading direction we obtained for this area is N97.7 ◦ E, resulting from the av. of 568 opening direction values.


Introduction
Quantifying the rate of crustal extension and the spreading direction across a rift zone is of paramount importance for a series of practical applications, which span from the assessment of seismic hazard [1], to the evaluation of factors that can contribute to facilitating magma uprising and thus to the assessment of volcanic hazard [2]. Anyway, the precise definition of the extensional rate, as well as the spreading direction, requires the knowledge of a series of parameters that include the architecture of the rift, the geometry of each single fault and extension fracture, the timing of rifting, and the direction of extension, and much more importantly, the data must be as precise as possible. Since it is difficult to obtain all these data with enough accuracy, usually rates of crustal extension are just estimated or approximated, and this is especially true in offshore rift zones such as oceanic ridges.

Methods
This section is devoted to illustrating the overall methodology we adopted in this work, starting from the first steps focused on planning of UAV surveys and related photo collection, to photogrammetry processing with the aim of building the 3-D model, the orthomosaic and the DSM, and, finally, to the collection of quantitative data on the above mentioned outputs.

Drone Surveys and Data Collection
Our approach, based on the application of SfM photogrammetry techniques through drone surveys (Figure 2b), has been made possible thanks to the DJI Phantom 4 Pro (Figure 2a), a versatile quadcopter, which, according to previous studies [17,19,63], has been proved to be suitable for surveying volcanic terrains with difficult logistic conditions (e.g., outcropping lavas). This UAV is equipped with a 20 Megapixels camera and an integrated Satellite Positioning System GPS/GLONASS (referred to the WGS84 datum). Metadata of the collected images are stored in EXIF files (Exchangeable Image File Format) which include information such as shutter speed, apertures, ISO and GPS coordinates (Datum WGS84).
Before starting all UAV surveys, a reconnaissance of the selected area has been done to locate any potential terrain obstacle and to evaluate the potential flight path accordingly. After this first

Methods
This section is devoted to illustrating the overall methodology we adopted in this work, starting from the first steps focused on planning of UAV surveys and related photo collection, to photogrammetry processing with the aim of building the 3-D model, the orthomosaic and the DSM, and, finally, to the collection of quantitative data on the above mentioned outputs.

Drone Surveys and Data Collection
Our approach, based on the application of SfM photogrammetry techniques through drone surveys (Figure 2b), has been made possible thanks to the DJI Phantom 4 Pro (Figure 2a), a versatile quadcopter, which, according to previous studies [17,19,63], has been proved to be suitable for surveying volcanic terrains with difficult logistic conditions (e.g., outcropping lavas). This UAV is equipped with a 20 Megapixels camera and an integrated Satellite Positioning System GPS/GLONASS (referred to the WGS84 datum). Metadata of the collected images are stored in EXIF files (Exchangeable Image File Format) which include information such as shutter speed, apertures, ISO and GPS coordinates (Datum WGS84). been taken into account: for example, the sun orientation, which is the same for all surveys (sun to the zenith) since they have been performed in time sequence, and the light conditions, which were always optimal and suitable for the camera ISO range . Also, we considered the wind speed and direction, being always parallel to the main and longer paths during each drone survey, and we maintained a constant path direction in each subarea in order to be sure to obtain the entire coverage of the whole big area.  Before starting all UAV surveys, a reconnaissance of the selected area has been done to locate any potential terrain obstacle and to evaluate the potential flight path accordingly. After this first visual evaluation, we created a polygonal shapefile in a GIS environment covering the target area, tailored on a satellite image, and we divided it into nine subareas that covered all the features that we wanted to map: extension fractures and normal faults affecting the hyaloclastites and post-LGM lava units (Figures 1 and 2c). The extension of each subarea has been defined to use a maximum of two UAV batteries for each area, considering a flight height between 80 and 95 m to reach a ground resolution of 2-3 cm for the orthomosaic. Shapefiles were converted into KML format and imported in DJI Ground Station PRO (https://www.dji.com/ground-station-pro), a free iPad app for DJI drones operations. This app allowed to generate nine different missions and to define and change for each of them the flight path direction, the flight height, as well as to fix camera parameters and directly check data acquisition on the iPad display. In fact, for all the surveys, fixed acquisition parameters were set (Table 1) to obtain proper images alignment and to reduce the distortions on the resulting models: (i) The overlap has been slightly decreased, with respect to previous authors [17,[64][65][66][67], to 85% along the path and 80% in lateral direction, in order to save time during each flight mission; (ii) photos have been captured using equal time interval mode and using constant velocity [17,68]: these two parameters are automatically managed by DJI Ground Station PRO in order to guarantee the set-up photos overlap at the given flight height. Additionally, a series of environmental conditions have been taken into account: for example, the sun orientation, which is the same for all surveys (sun to the zenith) since they have been performed in time sequence, and the light conditions, which were always optimal and suitable for the camera ISO range . Also, we considered the wind speed and direction, being always parallel to the main and longer paths during each drone survey, and we maintained a constant path direction in each subarea in order to be sure to obtain the entire coverage of the whole big area.
Furthermore, a crucial step in order to georeference the resulting models (orthomosaics and DSMs) and thus to allow the collection of high-precision quantitative measurements, was to collect World Geodetic System (WGS84) coordinates of ground control points (GCPs) in each drone-surveyed area. The placement and acquisition of GCPs was planned considering the specificities of the study area. We collected 37 GCPs (Figure 2c) distributed both at the corners and in the central part of each subarea to reduce the doming effect resulting from SfM processing and reduce systematic error [69][70][71][72]. Such GCPs can be artificial or natural objects, both must be easily recognizable on drone-captured photos ( Figure 3). Artificial GCPs deployment and recovery, on an extensive area, is really time-consuming since targets must be placed before starting the survey and retrieved afterward. For this reason, we used both artificial and natural targets as GCPs: artificial targets include objects with bright colours (e.g., Figure 3b), the corner of the backpack case of our drone (Figure 3c), and series of slices of white bread arranged to form an X on the ground (Figure 3d), whereas natural targets include massive rocks with defined geometric shapes (e.g., Figure 3e) and geometric lava flow margins (e.g., Figure 3f). Slices of white bread and natural targets were easily recognisable from UAV images and they have been used since the unnecessity to retrieve them could significantly speed up the survey.
The collection of GCPs has been performed with the Emlid Reach RS© single frequency receivers in RTK configuration (with centimetre-level accuracy), linked in real time by NTRIP to the Icelandic CORS network (ICECORS-https://www.lmi.is/en/icecors-network/), used as reference station. Because of the use of RTK/GNSS GPS information, we expect to obtain a model with centimetric accuracy, as already highlighted by several authors; for a more in-depth reading about the resulting accuracy we suggest the work proposed by the following authors: [73][74][75][76].

Photogrammetry Processing
UAV images have been processed using Agisoft METASHAPE (http://www.agisoft.com/), a commercial SfM software, and Agisoft Cloud beta service for data processing. The software was chosen for its user-friendly interface, spontaneous workflow, the good quality of the outputs [77][78][79] and the reasonable price of the educational license. These qualities have made it a broadly used application by the earth science community. SfM photogrammetric workflow, described in detail by [80,81], consists in a sequence of steps that from a set of overlapping images lead to the creation of a sparse and dense cloud, a 3-D tiled model, an orthomosaic and a DSM as final products (Table 2).

Photogrammetry Processing
UAV images have been processed using Agisoft METASHAPE (http://www.agisoft.com/), a commercial SfM software, and Agisoft Cloud beta service for data processing. The software was chosen for its user-friendly interface, spontaneous workflow, the good quality of the outputs [77][78][79] and the reasonable price of the educational license. These qualities have made it a broadly used application by the earth science community. SfM photogrammetric workflow, described in detail by [80,81], consists in a sequence of steps that from a set of overlapping images lead to the creation of a sparse and dense cloud, a 3-D tiled model, an orthomosaic and a DSM as final products (Table 2). The processing approach can be divided into four main stages: 1.
Photo alignment: an initial low-quality photo alignment using common key points in each image generated a sparse 3-D point cloud ( Figure 4a). The focal length and photo dimensions were automatically computed by the software and then used for the subsequent calibration of the core parameters of the camera (principal point coordinates, lens distortion coefficients).

2.
Georeferencing: GCPs were identified into the photos and markers have been assigned with the surveyed coordinates. This process allowed to scale and georeference the point cloud and thus to improve the accuracy of the final models. Then, the images were realigned using high accuracy settings and default values for Tie and Key points settings, 40,000 and 4000 respectively, and using Generic and Reference preselection settings.

3.
Dense Point Cloud Building: a 3-D dense point cloud was generated, using a mild depth filtering and medium quality settings, from the sparse point cloud (Figure 4b). 4.
DSM and orthomosaic generation: a DSM was generated from the dense point cloud, and successively the orthomosaic was produced using the corresponding DSM (Figure 5a,b).

5.
The software also allows to generate a high-quality 3-D tiled model, both including the mesh and the texture ( Figure 5c). The latter can be used as a scenario for immersive virtual reality [31]; here we suggest a Medium quality set-up with a tile size of 4096 × 4096 pixels.
All outputs refer to the geographic coordinate system (WGS84).

Data Collection on 2-D and 3-D Outputs
Once the resulting models had been processed, we imported the DSM and the orthomosaic in a GIS environment and mapped all fractures in the area, and remapped all post-LGM and LGM stratigraphic units (Figure 6a): these steps have been made possible thanks to the very high-resolution of the orthomosaic. Fractures have been further classified as normal faults when associated with a continuous and appreciable, and persistent, vertical offset > 0.5 m, calculated on the DSM, whereas they have been classified as extension fractures when associated with vertical offset less than 0.5 m.

Data Collection on 2-D and 3-D Outputs
Once the resulting models had been processed, we imported the DSM and the orthomosaic in a GIS environment and mapped all fractures in the area, and remapped all post-LGM and LGM stratigraphic units ( Figure 6a): these steps have been made possible thanks to the very high-resolution of the orthomosaic. Fractures have been further classified as normal faults when associated with a continuous and appreciable, and persistent, vertical offset > 0.5 m, calculated on the DSM, whereas they have been classified as extension fractures when associated with vertical offset less than 0.5 m. Apart from the DSM, the use of the 3-D tiled model, that can be defined as a virtual outcrop [82,83], has been crucial in order to perform visual checks in key points of each examined structure. Such models can be used for observation using an immersive virtual reality [31,64,65,84] and non-immersive virtual reality [85] approach as well as 3-D visualisation on a computer workstation [82]. For each set the av. value and standard deviation is represented as continuous and dashed lines respectively, the regional spreading direction range is reported as dark grey lines. (e) Fracture strike vs. length (for both extension fractures and normal faults), av. value is reported as a continuous line, the SD range is represented as a dashed line (box), and the regional spreading direction range is represented as a dark grey rectangle.
Different types of quantitative data have been taken: along normal faults, we quantified the strike and the vertical offset, by measuring differences in altitude along 50-m-spaced topographic profiles traced on the DSM, orthogonally to the fault scarp, whereas along extension fractures we collected the strike, the amount of dilation and opening direction values. More in detail, wherever piercing points on each side of the extension fractures were clearly visible, we traced a line connecting them on the orthomosaic, and we calculated its strike and length to obtain the direction and the amount of opening, respectively. All factors that could alter the accuracy of the resulting values have been avoided, by making sure to collect data only at fractures not affected by erosion and sedimentation, without falling rocks on their sides and far away from fault rims to avoid local gravity effects. For these quantitative measurements, we also used the SfM-derived 3-D model for some For each set the av. value and standard deviation is represented as continuous and dashed lines respectively, the regional spreading direction range is reported as dark grey lines. (e) Fracture strike vs. length (for both extension fractures and normal faults), av. value is reported as a continuous line, the SD range is represented as a dashed line (box), and the regional spreading direction range is represented as a dark grey rectangle. Different types of quantitative data have been taken: along normal faults, we quantified the strike and the vertical offset, by measuring differences in altitude along 50-m-spaced topographic profiles traced on the DSM, orthogonally to the fault scarp, whereas along extension fractures we collected the strike, the amount of dilation and opening direction values. More in detail, wherever piercing points on each side of the extension fractures were clearly visible, we traced a line connecting them on the orthomosaic, and we calculated its strike and length to obtain the direction and the amount of opening, respectively. All factors that could alter the accuracy of the resulting values have been avoided, by making sure to collect data only at fractures not affected by erosion and sedimentation, without falling rocks on their sides and far away from fault rims to avoid local gravity effects. For these quantitative measurements, we also used the SfM-derived 3-D model for some checks to better recognize all fault scarps and visualize piercing points. Finally, all collected data have been used for two aims: defining the spreading direction in the area and quantifying the overall extension rate during late Pleistocene and Holocene times. We calculated all horizontal dilation values along N105 • E-trending transects, in order to obtain the cumulative dilation, stretch, and extension rate in the two main outcropping successions (LGM and post-LGM). These transects have been located in order to cover both the LGM unit and the post-LGM units. The N105 • E is the mean value of the spreading direction given by DeMets et al. [54] and Hjartardóttir et al. [60], both considering geological data for a long-term time window.

High Resolution Outputs from Photogrammetry Processing
The areal coverage of our study area is about 2.7 km 2 subdivided into nine different surveys, performed at 80-95 m of height from the ground; the surveys were carried out during the summer 2019 campaign, during which we collected a total of 6068 photos during a total flight time of 4.68 hrs. Because of the above described workflow, such photos were used to reconstruct a high quality orthomosaic and a DSM (Figure 5a,b), respectively characterized by a resolution of 2.59 and 10.40 cm/pixel; the resulting GeoTIFFs have a storage dimension of 15.8 and 1.11 GBs respectively, the orthomosaic has been converted to JPEG 2000 (JP2) for data collection in GIS. The 3-D tiled model is characterised by a texture resolution that is identical to the orthomosaic and a storage dimension of 9.87 GBs. Regarding the processing time, as detailed in Table 2, all outputs have been generated in 46.8 h, where the most time-consuming processing have been the dense cloud and 3-D tiled model generation.

New Geological-Structural Data
During the present work, we mapped a total of 1355 fractures classified as normal faults (86) and extension fractures (1269): these structures strike from NNW-SSE to NE-SW, but are mostly in the range N0-20 • E, whereby the av. strike is N8.2 • E, as shown in the rose diagram in Figure 6b. The longest fracture reaches a length of 1156.1 m, whereas the overall av. length is 28.3 m (Figure 6e). Although the whole Krafla rift has an about N-S elongation, once you define the structures in detail it emerges that several single faults have a bimodal strike given mostly by N-S and NNE-SSW segments, with a curvilinear trace in plan view.
Owing to the high resolution of the new DSM and orthomosaic, the number of mapped fractures is dramatically greater with respect to the 45 main fractures documented on the 1:100,000 scale geological map [39] (Figure 1b). In addition to the newly mapped fractures mentioned above, we have also redrawn the geological map of the area, tracing with very high detail the geological limits, as shown in Figures 6a and 7.
In regard to normal faults, they mainly strike from N-S to NE-SW, dipping towards WNW and ESE (Figure 6c), and show an overall av. and maximum length of 167.9 m and 1156.1 m respectively (Figure 6e). On the other hand, extension fractures mainly strike from NNW to NNE, with an av. value of N8.1 • E (Figure 6d), and an overall av. and maximum length of 18.8 m and 240.7 m respectively (Figure 6e). In order to better define the deformation processes that have been occurring since late Pleistocene times, we collected a series of quantitative data at a total of 852 sites, 284 along normal faults and 568 along extension fractures. Regarding normal faults, we collected 284 vertical offsets, thus obtaining a maximum value of 42 m along faults affecting LGM units and of 19 m along faults affecting post-LGM units (Figure 7). Regarding extension fractures, we collected a total of 1704 structural data, consisting of strike, opening direction and fracture dilation, at 568 different sites  Figure 6a. The rose diagram shows both extension fracture strikes and opening directions. For each set of data, the av. value and standard deviation is represented as continuous and dashed lines respectively, the regional spreading direction range is reported as dark grey lines.  Figure 6a. The rose diagram shows both extension fracture strikes and opening directions. For each set of data, the av. value and standard deviation is represented as continuous and dashed lines respectively, the regional spreading direction range is reported as dark grey lines.
Opening direction values are in the range N46.3-142.3 • E, with a mean value of N97.7 • E and a SD of 13.3 • : only 109 values fall within the spreading direction range (N104-115 • E), which is based on spreading direction estimates given in the literature for the NVZ [17,54,[57][58][59][60][61]. We also measured the extension fracture azimuth that ranges N48.1 • W-N61.5 • E, with a mean value of N7.9 • E and a SD of 18.5 • : only 119 values show an azimuth that is compatible with the spreading direction range. If we consider both the extension fractures opening direction and strike, only 43 out of 568 are compatible with the spreading direction range N104-115 • E (Figure 8a).

Rift Architecture and Kinematics
The studied rift segment is located along the north-western part of the KFS and is characterised by the presence of grabens arranged en-échelon with a left-stepping geometry. In the study area, two main grabens are present: one, to the north, affects the Holocene lavas only, and disappears in correspondence of the other graben. This second graben is located further south and affects mostly the hyaloclastic ridge and, to a lesser extent, prolongs into the surrounding Holocene lava succession.
The southern graben is a few tens of meters wide at its northern and southern terminations, whereas its width increases up to 300 m in its central part. This geometry in plan view cannot be related to the presence of the volcanic ridge respect to the surrounding flat area for two reasons: (i) Classical studies showed that the presence of a volcanic edifice strongly influences fault orientation by deflecting and converging faults that tend to become radial at the cone, assuming an hourglass pattern [87,88]; in our case instead, faults tend to diverge in correspondence of the ridge. (ii) In plan On the other hand, by increasing the extension fractures strike, we observe an increase in opening directions with a sort of clockwise rotation of values (Figure 8a). In fact, considering N-S to NE-SW fractures, opening direction values range from N73.7 • E to N142.3 • E, whereas NW-SE fractures are characterized by opening direction values ranging from N46.4 • E to N125.5 • E. On the other hand, by increasing the fracture strike, we notice a decrease in the lateral component, that changes from right-lateral to left-lateral (Figure 8e). Finally, values of the amount of opening slightly increase for fractures characterized by a lower component of motion (Figure 8f).
We also measured the extensional strain in the area by cumulating the dilation values obtained at both extension fractures and normal faults, along three N105 • E-striking transects, each with a length of 0.62 km (Figure 7). To calculate the dilatational component at normal faults, we used the measured vertical offset, considering a dip of the fault plane of 70 • [86]. Along the northernmost transect, which cuts structures affecting 11-12 ka BP lavas, we measured a cumulated dilation of 16.6 m, corresponding to a stretch of 1.027; along the central one, which cuts only faults affecting the hyaloclastites, we obtained a cumulated dilation of 29.3 m, corresponding to a stretch value of 1.049; along the southernmost transect, which cuts structures affecting 12 ka lavas, the cumulated dilation has a value of 11.2 m, with a stretch value of 1.018. The overall dilation, calculated in the central part of the area along a 1.23 km-long transect and considering all fractures affecting LGM and post-LGM units, is 36.4 m, corresponding to a stretch of 1.0304.

Rift Architecture and Kinematics
The studied rift segment is located along the north-western part of the KFS and is characterised by the presence of grabens arranged en-échelon with a left-stepping geometry. In the study area, two main grabens are present: one, to the north, affects the Holocene lavas only, and disappears in correspondence of the other graben. This second graben is located further south and affects mostly the hyaloclastic ridge and, to a lesser extent, prolongs into the surrounding Holocene lava succession.
The southern graben is a few tens of meters wide at its northern and southern terminations, whereas its width increases up to 300 m in its central part. This geometry in plan view cannot be related to the presence of the volcanic ridge respect to the surrounding flat area for two reasons: (i) Classical studies showed that the presence of a volcanic edifice strongly influences fault orientation by deflecting and converging faults that tend to become radial at the cone, assuming an hourglass pattern [87,88]; in our case instead, faults tend to diverge in correspondence of the ridge. (ii) In plan view, inclined normal faults tend to diverge in correspondence of a topographic high, like a volcanic cone, but in our case it is possible to observe that faults diverge and graben width increases to the north exactly in correspondence of the contact between the hyaloclastites and the lava flows, and to the south in the flat area occupied by the lavas. Moreover, the difference in altitude between the ridge and the surrounding flat area is modest and cannot account for the bulge of fault traces as simply because of fault outcrop following topography. We propose that the increase of graben width in correspondence of the hyaloclastite ridge could instead be related to the mechanical contact between the stiffer lava succession and the softer hyaloclastites: the rheological contrast between the two rock units leads to faults that tend to follow the outline of the contact, favouring concentric faults. This is similar to what has been noticed after simulating a magma chamber and diffuse rift extension [88]: in this case, the presence of a magma chamber affects the faults orientation that close to the intrusion edge change their direction, mimicking the intrusion shape by curving around it.
Focusing on the kinematics of extension fractures, the collected data suggest an av. spreading direction of N97.7 • E, with a range of N46.3-142.3 • E for the Holocene times (Figure 7). This av. value is lower than those given in the literature that range between N104 • E and N115 • E [54,[57][58][59][60]. The av.
spreading value calculated in the present research is closer to those that take into account geological data that cover a larger time window [17,54,60,61]) with respect to only GPS data [57][58][59], suggesting a slight overall anticlockwise rotation of the spreading direction in the NW part of the Krafla rift.
The overall fracture azimuth is suggested to be orthogonal to the spreading vector range discussed in the above section: in fact, longer and more developed fractures show an azimuth N-S to NNE-SSW (Figure 6a,b,e), as suggested for the active Krafla fissure swarm [60].
Finally, we paid attention to the presence of strike-slip components along the extension fractures, where the av. value suggests an overall very limited left-lateral component, even though the majority of data are characterized by values with right-lateral components (Figure 8e,f). We believe that this component is the result of local perturbations by dyking at shallow depths, as shown for the 2014-2015 Bárðarbunga dyking event [89] and as suggested in the ThFS [16,17]. In support of this hypothesis, the area has been subjected to recent historical rifting events characterized by N-propagating dykes [48], and Tibaldi et al. [37] described textbook examples of dyke-induced deformation at the surface in the eastern part of the area presented in this study.

Rift Extension Rate
The total extension measured along the fractures cropping out in the Holocene lava succession is, along two transects, 16.6 m and 11.2 m, whereas the total extension in the hyaloclastites is 29.3 m, coherently with the older stratigraphic age of the latter. Assuming the age of 11.5 ± 0.5 ka BP for the post-LGM lava units, and of 20.6 ± 8.5 ka BP for the hyaloclastite deposits, we obtain an extension rate of 1.4 ± 0.7 mm/yr and 1.7 ± 0.7 mm/yr, respectively. Our measured extension rates of 1.4 mm/yr and 1.7 ± 0.7 mm/yr for the shorter and longer periods are coherent with each other, but much slower than the GPS velocity field measured after the 1975-1984 rifting/dyking period (3-4.5 cm/yr). We thus consider that our measures do represent the real long-term deformation field of the western part of the Krafla rift. This represents only a fraction of the total extension that characterizes this rift, which in turn is a fraction of the extension of the whole NVZ of 1.8-2.3 cm/yr [54,55]. From another point of view, focusing on the amount of horizontal dilation, the stretch in the area is 1.018-1.027 for post-LGM units, 1.049 for LGM unit and is 1.030 if we consider both. In the Krafla area, a stretch value of 1.009 was previously calculated [90] along a 30-m-long transect in an area located just north of the 1984 Krafla lava flow, totalling a horizontal dilation of 30 m for the last 10,000 yrs. Even though, owing to the high resolution of our SfM-derived models, we collected a larger number of data, also at fractures with centimetric dilation, we think that both values describe the contribution of both tectonic and magmatic forces (dyke intrusion) in dictating deformation at the surface in terms of extension fractures opening and normal faulting. Differently, Paquet et al. [91] calculated a stretch ratio of 1.036-1.046 across the Tertiary Alftafjordur dyke swarm, describing the deformation that lasted 1 Ma or more. The stretch value that we calculated on the LGM units resulted like Paquet et al. [91], even though they considered a much longer time span and a greater depth of analysis since they studied an ancient rift zone exposed by glacial erosion. We finally suggest that a stretch value <1.030 can be considered as representative of the influence of both tectonic and magmatic forces; on the contrary, based on data from the present work and from [91], any classification for stretch values >1.030 is controversial and must be subject to further research.

Methodological Aspects
Even though the application of the method can have some limitations (e.g., the age of structures should be kind of 'recent', the research area should be well exposed with a low vegetation coverage, basic geological researches should be completed ahead of the survey), our work confirms the usefulness of the UAV-based SfM as a support for classical fieldwork to collect massive amounts of high resolution quantitative data, here aimed at reconstructing late Pleistocene-Holocene deformation in an active volcano-tectonic region, being here a key area within the active Krafla rift. At a general level, the use of UAVs in geoscience has increased over the past decade [92][93][94][95][96][97][98][99][100][101]. The affordability of such method is one of the main reasons why the geosciences community broadly adopts it: UAV data acquisition represents, in fact, a less expensive alternative to other techniques, such as Terrestrial Laser Scanning (TLS), Airborne Laser Scanning (ALS) and LiDAR [102,103]. Another reason for which this methodology has been broadly used is the possibility to save time with respect to fieldwork, especially in vast areas and along structures that can be tens of km long [15,104], and to enable to reach key geological sites that are generally inaccessible because of hard logistic conditions, like volcano edifices affected by caldera or lateral collapse [105][106][107]. Furthermore, the rapidity of data capturing, the quick processing and the high accuracy of the output model generated through the UAV-SfM technique allow to observe at a centimetric detail the characteristics of the ground surface and outcrops, for example by catching also micro-topographic variations. As an example of the great amount of measurements that can be collected through this method, Bonali et al. [17] were able to collect 1098 measurements at 387 sites in a volcano-tectonic field, comprising azimuth, dilation and opening direction of extension fractures by using a commercial-grade UAV and applying SfM algorithms, thus providing a contribution to the reconstruction of the tectonic evolution in such an active rift zone in NE Iceland.
Another advantage of this methodology is represented by the continuous improvements over the years of the consumer-grade UAVs, in terms of flight stability, battery capacity and camera quality, owing to the constant release of new products: for example, in this work we used the DJI Phantom 4 Pro (20 Megapixel and ISO range 100-12800), which has been released less than a year later with respect to the DJI Phantom 4 (12 Megapixel and ISO range 100-1600). This implies that in the next years, owing to such fast improvements of the UAV technology, geological studies will achieve better quality by flying at a higher altitude and covering a wider area. Finally, the SfM technique could be applied with a low-cost workflow, particularly suitable for field research, which allows any user to process the data [108].
Through a methodological point of view, we tested two novelties in our approach: (i) we slightly decreased the image overlap to save time during each single flight mission. This overlap decrease resulted in 30% of saved time. Considering our overall survey, we saved about 2 h, which is plenty of time in areas usually affected by wind and rain. (ii) We mainly used natural, ecological, or already available artificial targets for GCPs collection, this action halved the time usually necessary for this fundamental task. We obtained excellent models as output where we collected very high-resolution structural data, especially regarding piercing points along extension fractures. Of course, both novelties can be applied to further researches and areas. In addition, the area of this research is characterized by an overall morphology that is more complicated than in Bonali et al. [17]: in fact normal faults have a different morphological expression due to the age and characteristics of volcanic units. Finally, the presence of the Hituholar volcanic edifice, with an elevation of more than 100 m from the surrounding area, has made it more challenging to plan and execute UAV missions. Even so, we reached excellent resolution and accuracy in all SfM outputs, allowing us to collect all data necessary to define spreading direction and stretch in the studied area, enhancing the importance of this methodology in this kind of studies.
The UAV-based SfM also presents a few limitations, such as the dependence on the battery life/flight time, weather conditions and legal flight restrictions. Also the necessity of positioning some GCPs, in order to reference the model accurately, could be time spending and thus could slow down the survey, especially over vast areas. Smith et al. [109] recommended the use of at least 3 GCPs; other authors used a higher number of GCPs (up to 95 across a 100 m grid), but reducing the image overlap to 60-70% [110]. With a greater overlap, on the other hand, a lesser number of GCPs can be used [68,111]. The number of GCPs can also be reduced if the geometry network results sufficiently strong [72] or by using an UAV incorporating onboard RTK-GNSS sensors that would allow to generate models with centimetric spatial accuracy without placing GCPs.
In our work, we used 37 GCPs, both along the corners and at the centre of every area, in order to reduce the 'doming' effect and reduce systematic errors [69][70][71][72]. Another strategy to reduce the 'doming' effect on the resulting DTM/DSM [73], consists of collecting aerial photos with an oblique inclination with respect to the ground [112]. Nevertheless, this strategy results in greater consumption of the battery power and so it is useless for commercial UAVs equipped with low-cost GPS chipset, where GCPs collection is crucial to georeference the model; it is thus more efficient to design a GCPs survey as proposed in the present work.
We thus selected our flight strategy to reach a good compromise between battery consumption, survey time, area extension, data resolution and accuracy of the models. In just two days of fieldwork and almost 46 h of processing of about 7000 UAV pictures, we mapped an area of 232 Ha (2.3 km 2 ). Moreover, the achieved resolution of the orthomosaic (2.6 cm/px) and the DSM (10 cm/px) resulted sufficiently detailed for surveying extension fractures and collecting plenty of data of opening vectors, down to centimetric values of dilation. This high accuracy opens the possibility to carry out repeated UAV surveys across the same area at regular time intervals, in order to study the present-day deformation and observe its evolution over time. This could be applied, for example in the KFS, by comparing the resulting models (orthomosaics and DSMs) before and after a rifting event, if the superficial deformation is larger than the model resolution. Regardless, additional improvements are suggested to define the most efficient survey approach for DSMs reconstruction, mainly devoted to measuring vertical offset at normal faults. The fast evolution of UAVs, equipped with RTK positioning system, can represent a solution for the above-presented limitations, increasing the accuracy of the model, flying also at greater heights. Finally, the resulting 3-D models can be navigated in the immersive virtual reality, useful both for scientific research and teaching and divulgation purposes [17,31,64,65].

Conclusions
We used drone surveys coupled with structure-from-motion photogrammetry, introducing two novelties to increase the efficiency of drone surveys, to analyse a remote area characterised by rough terrains in the NW part of the KFS (NE Iceland). This is a volcanotectonic rift composed of eruptive centres and N-S to NNE-SSW fissures and normal faults. The central part of the surveyed sector has a hyaloclastite ridge composed of deposits dated, on a stratigraphic basis, to the Weichselian High Glacial (LGM here dated 29.1-12.1 ka BP). These are surrounded by a series of flat lava flows dated at 11-12 ka BP, and locally by 1975-1984 AD lavas.
The integration of the data enabled us to recognize that this segment of the KFS is made of grabens arranged en-échelon with a left-stepping geometry. A major graben increases in width in correspondence of the hyaloclastite ridge; we interpret this geometry as resulting from the mechanical contrast between the stiffer lava succession and the softer hyaloclastites, which favours the development of concentric faults.
The direction of extension has also been measured at hundreds of sites, resulting in an E-W to WNW-ESE trend, normal to the fissure and fault strikes. The total extension results to be 16.6 and 11.2 m along the faults and fissures cropping out in the Holocene lava succession, whereas the total extension measured in the hyaloclastites is 29.3 m. This indicates an extension rate of 1.4 mm/yr in the Holocene lavas and 1.7 ± 0.7 mm/yr in the latest Pleistocene hyaloclastite deposits. We consider these data to represent the real values of long-term extension in this part of the Northern Volcanic Zone, much slower than the GPS velocity field measured around the Krafla system that resulted from the 1975-1984 rifting and dyking events. The av. spreading direction here quantified is N97.7 • E.