Structure of massively dilatant faults in Iceland: lessons learned from high resolution UAV data

15 Normal faults in basalts develop massive dilatancy in the upper few hundred meters below the Earth's surface with corresponding interactions with groundwater and lava flow. These massively dilatant faults (MDF) are widespread in Iceland and the East African Rift, but the details of their geometry are not well documented, despite their importance for fluid flow in the subsurface, geo-hazard assessment, and geothermal energy. We present a large set of digital elevation models (DEM) of the surface geometries of MDF with 5-15 cm resolution, acquired along the 20 Icelandic Rift zone using unmanned aerial vehicles (UAV). Our data present a representative set of outcrops of MDF in Iceland, formed in basaltic sequences linked to the Mid Ocean Ridge. UAV provide a much higher resolution than aerial/satellite imagery and a much better overview than ground-based fieldwork, bridging the gap between outcrop scale observations and remote sensing. We acquired photosets of overlapping images along about 20 km of MDF and processed these using photogrammetry to create high 25 resolution DEMs and ortho-rectified images. We use this dataset to map the faults and their damage zones to measure length, opening-width and vertical offset of the faults and identify surface tilt in the damage zones. Ground truthing of the data was done by field observations. Mapped vertical offsets show typical trends of normal fault growth by segment coalescence. However, openingwidths in map-view show variations at much higher frequency, caused by segmentation, collapsed relays and tilted 30

MDF in Iceland, formed in basaltic sequences linked to the Mid Ocean Ridge.
UAV provide a much higher resolution than aerial/satellite imagery and a much better overview than ground-based fieldwork, bridging the gap between outcrop scale observations and remote sensing. We acquired photosets of overlapping images along about 20 km of MDF and processed these using photogrammetry to create high 25 resolution DEMs and ortho-rectified images. We use this dataset to map the faults and their damage zones to measure length, opening-width and vertical offset of the faults and identify surface tilt in the damage zones. Ground truthing of the data was done by field observations.
Mapped vertical offsets show typical trends of normal fault growth by segment coalescence. However, openingwidths in map-view show variations at much higher frequency, caused by segmentation, collapsed relays and tilted 30 blocks. These effects commonly cause a higher than expected ratio of vertical offset and opening-width for a steep normal fault at depth.
Based on field observations and the relationships of opening-width and vertical offset, we define three endmember morphologies of MDF: (i) dilatant faults with opening-width and vertical offset, (ii) tilted blocks (TB), and (iii) opening mode (mode I) fissures. Field observation of normal faults without visible opening invariably shows that 5 these have an opening filled with recent sediment. TB dominated normal faults tend to have a largest ratio of opening-width and vertical offset. Fissures have opening-widths up to 15 m with throw below a 2 m threshold.
Plotting opening-width versus vertical offset shows that there is a continuous transition between the endmembers.
We conclude that for these endmembers, the ratio between opening-width and vertical offset R can be reliably used to predict fault structures at depth. However, fractures associated with MDF belong to one larger continuum 10 and consequently where different endmembers coexist, a clear identification of structures solely via the determination of R is impossible.
These massively dilatant faults (MDF) are common in rift zones such as the Icelandic Rift and the East African Rift or in active volcanic systems such as Hawaii (Acocella et al., 2003;Gudmundsson, 1987b;Martel and Langley, 2006;Rowland et al., 2007). MDF guide the flux of water, magma or hydrocarbons and are therefore of interest for applications such as geo-hazard assessment, hydrocarbon exploration and geothermal energy 20 production (Crider and Peacock, 2004;Faulkner et al., 2010;Ferrill and Morris, 2003;Grant and Kattenhorn, 2004;Gudmundsson, 1987a;Kettermann et al., 2015Kettermann et al., , 2016Rowland et al., 2007).
The surface geometries have been described including dilatancy, tilted blocks (i.e. rigid, detached blocks tilted towards the hanging wall of a fault) and extension fractures. However, many of the observations are based on local measurements considered representative for the regional structures.
In this study we investigate the structure and evolution of massively dilatant faults in Iceland (Fig. 1, 2) by identifying and characterizing surface geometries at the regional scale at centimeter resolution. We achieve this by extracting data of MDF from high resolution maps generated from unmanned aerial vehicle (UAV) based photogrammetry. Mapping faults in centimeter resolution over kilometers length allows for bridging the gap between outcrop scale and regional observations. This enables us to quantify the geometry of the studied faults at 5 high detail over the entire fault lengths. The ultimate goal of this work is to introduce a new classification scheme that correlates a ratio of measured fault aperture and fault throw with actual underlying fault structures that are often overprinted by sedimentation or erosion. To describe different types of discontinuities we use the terminology of Peacock et al. (2016).

Massively Dilatant Faults 10
Massively Dilatant Faults (MDF) ( Fig. 1) with several meters to tens of meters aperture can form close to the surface in cohesive rocks. They are characterized by distinct geometries such as sub-vertical fault scarps, rotating blocks and fractures that remain open (or often filled with sediment or rubble) up to hundreds of meters deep (Acocella et al., 2000;Grant and Kattenhorn, 2004). Existing studies of MDF in outcrops (Holland et al., 2006) and analogue models show that MDF form in layers where the ratio of rock strength to effective stress is 15 sufficiently high (van Gent et al., 2010). With depth, MDF transition to shear faults due to the increase of effective stress.
Geometries of MDF have been described for several sites in Iceland with respect to fracture length, opening, throw, obliquity and segment linkage (Acocella et al., 2000;Bonali et al., 2019a;Bubeck et al., 2017;Gudmundsson, 1987bGudmundsson, , 1987avon Hagke et al., 2019;Tentler and Mazzoli, 2005;Trippanera et al., 2015). Relationships of length 20 and opening are complex (Hatton et al., 1994), with largest openings at the fault center (Gudmundsson, 1987a) and smaller opening and throw closer to the fault tips (Tentler and Mazzoli, 2005). Correlations of opening and throw are weak, but show different distributions at fault tips and centers, possibly reflecting different stages of fault growth (Gudmundsson, 1987a;Tentler and Mazzoli, 2005). However, a direct relation of larger vertical offsets and larger dilatancy has been suggested by Acocella et al. (2003) and Trippanera et al. (2015). 25 Several models explaining the growth of MDF exist, ranging from (1) upwards propagation of the fault (Grant and Kattenhorn, 2004), perhaps linked to dike intrusion (Trippanera et al., 2015) over (2) linkage of shear faults at depth with tensile fractures at the surface (Abe et al., 2011;Hardy, 2013;Vitale and Isaia, 2014) to (3) nucleation at the surface and downward propagation of the fault until tensile failure is no longer possible (Acocella et al., 2003;van Gent et al., 2010). The transition from tensile to shear fractures can be envisioned as a broad zone 30 reaching down to depths of more than 1000 m, depending on the mechanical stratigraphy of the fractured units, as well as on fault kinematics Kettermann et al., 2019, Gudmundsson andBäckström, 1991).
In the following we quantify the surface geometry of MDF along representative set of faults. We develop a method to distinguish different characteristic surface structures, and show that it is possible to make inferences on fault processes and geometries in the shallow subsurface based on analysis of surface data only. 5

Geological background and study sites
Iceland is a volcanic island in the Atlantic Ocean on the Mid-Atlantic Ridge separating the Eurasian and North American plates. It is linked to a deep mantle plume (Einarsson, 1991;Lawver and Muller, 1994;Vink, 1984;Wolfe et al., 1997), with associated melt production forming the Icelandic shelf  with a local crustal thickness of at least 25 km (Allen et al., 2002). It is located between the Reykjanes Ridge segment in 10 the SW and the Kolbeinsey Ridge in the N (Fig. 2). Minimum horizontal stresses in Iceland are oriented orthogonal to the rift axes (Ziegler et al., 2016), with a NW-SE trend on the Reykjanes Peninsula Ridge and E-W trend in the North Volcanic Zone. The orientation of the faults, fissures and dikes follow the trend of the rift axes ( Fig. 2; e.g. Grant and Kattenhorn, 2004;Gudmundsson, 1983Gudmundsson, , 1987aGudmundsson, , 1987bHjartardóttir et al., 2012;Opheim and Gudmundsson, 1989). The bedrock is mostly of volcanic origin with ages increasing with distance from the 15 rift (Fig. 2). The succession of lava flows with cooling joints, paleo-soils and hyaloclastite result in a complex mechanical stratigraphy. Faults and fissures crosscutting this mechanically heterogeneous sections reactivate the pre-existing cooling joints close to the surface, leading to complex geometries (Forslund and Gudmundsson, 1991;Gudmundsson, 1987a;Gudmundsson and Bäckström, 1991;Hatton et al., 1994).
The sites chosen for this project include the West Volcanic Zone (WVZ) and North Volcanic Zone (NVZ) (Fig.  20 2). The less than 8 Ma old NVZ (Saemundsson, 1974) is composed of seven volcanic systems, each with a central volcano and associated N-NNE striking fault and fissure swarms (Gudmundsson, 1995;Hjartardóttir et al., 2015;Saemundsson, 1974). We selected these sites because they are representative for the variability of faults on Iceland, as they include faults in purely extensional, but also in oblique rift kinematics. They offer the best outcrop conditions with well-defined structures with only minimal vegetation, soil cover or erosion. Associated with the 25 Krafla volcano in the NVZ, the Krafla fissure swarm stretches ca. 40 km south and 50 km north, mostly in postglacial lava flows (Hjartardóttir et al., 2012), with faults and fissures reaching a maximum opening-width of 40 m and vertical offsets up to 42 m (Opheim and Gudmundsson, 1989). Two recent rifting episodes are documented, with eruptions in 1724-1729, (the Mývatn Fires) and 1975 -84, (Krafla Fires), both accompanied by strong earthquakes and movement along active faults. The horizontal displacement of the rift system during the 30 Krafla Fires was about 8 m, equal to ca. 500 years of plate divergence (Hollingsworth et al., 2012;Tryggvason, 1994).
The Theistareykir fissure swarm is located in the rift zone within the NVZ and composed of N-S-striking Holocene fissures and normal faults (Gudmundsson et al., 1993). The westernmost normal fault of the Theistareykir fissure swarm, known as the Gudfinnugja Fault, connects with the Húsavik-Flatley fault which also offsets the Kolbeinsey 5 Ridge in the North (Fig. 2) (Gudmundsson et al., 1993;Pasquarè Mariotto et al., 2015;Tibaldi et al., 2016). Further interactions of the Theistareykir fissure swarm with the Tjörnes fracture zone have been suggested by Bonali et al. (2019a) due to variations in strike direction.
The Thingvellir fissure swarm is linked to the Pleistocene Hengill volcanic system by the continuity of the faults and the documented ground movement during the last rifting episode in 1789 (Saemundsson and Saemundsson, 10 1992). The volcanic system includes ca. 100 associated fissures and faults, of which some reach opening-widths > 60 m and vertical offsets of 40 m, representing the largest postglacial structures of Icelandic rift zones (Gudmundsson, 1987b). The fissure swarm consists of two Holocene main faults, known as Hrafnagjá and Almannagjá, that envelop lake Thingvallavatn. Almannagjá, which is subject to this study, has been described by Gudmundsson (1987b) as 7.7 km long and locally up to 64 m wide. The southern tip consists of several en-échelon 15 extension fractures, which are interpreted to be related to an older weakness, e.g. a Pleistocene fault underneath (Gudmundsson, 1987b), as further north the fault consists of several parallel fractures and parallel extension fractures at the northern tip (Gudmundsson, 1987b).
The Vogar fissure swarm is located on the NE of the Reykjanes Peninsula in postglacial lava flows (Fig. 2 a,b).
The geometry of the fractures of the fissure swarm has been characterized as anastomosing or sinuous (Clifton and Kattenhorn, 2006;Clifton and Schlische, 2003;Gudmundsson, 1987a).

Methods 25
We combined satellite and airborne imagery and used drone imagery to create 3D surface models and digital elevation models (DEMs). These different methods focus on different scales to acquire data, thus making it possible to bridge the gap between cm-and km-scale observations.

Satellite-borne data
To identify areas of interest, we used published datasets of satellite and airborne imagery. Google Earth was used in combination with the aerial photographs from Loftmyndir Inc., that are freely accessible via https://www.map.is/base/. As Iceland has been subject to several landscape shaping volcanic eruptions during the time span that is covered by remote sensing data, we included aerial imagery dating back to the 1950s provided 5 by the National Land Survey of Iceland in our preliminary remote sensing. Some examples of young volcanic eruptions include the Krafla Fires (Hjartardóttir et al., 2012;Tryggvason, 1986) with nine eruptions between 1975-1984(Thordarson and Larsen, 2007, the 1991 and 2000 eruptions of Hekla Rose et al., 2003), or the Holuhraun eruption during 2014 -2015, the largest Icelandic eruption since more than two centuries (Geiger et al., 2016;Müller et al., 2017;Schmidt et al., 2015). Published mappings of faults and fissures 10 of the NVZ (Hjartardóttir et al., 2012) and the Reykjanes Peninsula (Clifton and Kattenhorn, 2006) were further used to identify possible sites. TanDEM-X WorldDEM™ tiles with a resolution of 12 m were used at a later stage to complement our own DEMs, which are highly resolved but only able to cover respectively smaller areas. With TanDEM-X WorldDEM™, general surface slopes were identified to quality check our own models in order to avoid typical error sources such as doming (James and Robson, 2014) known to possibly accompany Structure 15 from Motion-based DEMs. Furthermore, the elevation models were used to aid and back the interpretations in areas where we do not have coverage with own spatial data.

Unmanned aerial vehicle photogrammetry
Unmanned aerial vehicles (UAVs) have been used increasingly in geosciences within the last years, using commercial ready to fly sets or self-made products. They have been applied successfully to e.g. aid the mapping 20 of faults and joints (Bonali et al., 2019b;Trippanera et al., 2019;Vasuki et al., 2014), landslides (Niethammer et al., 2012), or moss beds (Lucieer et al., 2014). The data for this study was acquired using three different UAVs: the DJI Phantom 4 and Mavic pro with 12 MP sensors and the Phantom 4 Advanced, with a 20 MP sensor. We took front-and sideways overlapping photographs of the faults with the cameras facing 90° downwards, covering the area to be mapped, and added oblique photographs to reduce possible 'doming' effects as suggested by (James 25 and Robson, 2014). Due to the large distances we covered, most flights to acquire the photographs were undertaken manually. For later processing in Agisoft Photoscan, we aimed for a frontal overlap of at least 70 % and a sidelap of at least 50 %, being much higher in practice, especially at the fractures where > 9 perspectives are achieved for most cases. We varied the altitude from which the photographs were taken between 30 m and 70 m above ground according to the estimated accuracy of the resulting DEM and the general dimension of the area to be mapped. We 30 thereby kept the flight altitudes as low as possible, depending on the dimension of the survey area. The focus of our photographs was set on the fractures and the adjacent areas to capture structures linked to the fracture geometry, such as TB or the damage zone and identify possible surface variations such as topographic slopes.
Our method is similar the one used by Bonali et al. (2019b) who collected data from 50 and 100 m altitude and have shown that with these setting one obtains sufficient resolution for detailed analysis of faults and fissures. We 5 did not place further ground control points (GCPs) and relied on the integrated GPS receiver of the UAV in favor of time efficiency, since absolute elevations are not relevant to our research goals. We used Agisoft Photoscan as processing tool to align the photographs on high settings and create a sparse 3D point cloud. As second step, a dense point cloud was created, that also served as basis for the DEM, in medium or high resolution, in case the medium setting resulted in DEMs with resolutions > 15 cm/pixel. The resulting DEMs were further used as basis 10 for ortho-mosaics with half the spatial resolution, thus representing the lower threshold of the accuracy of our measurements. The combination of our UAV images and the processing software allowed us to create digital elevation models and ortho-rectified mosaics with a resolution varying between 5 and 15 cm per pixel, depending on the i) the altitude the photographs were taken from, ii) the count of perspectives achieved per object and iii) the quality settings used during the processing. 15 GPS tags by the UAV onboard receiver lead to an absolute horizontal accuracy of ~ 5 m per photograph and a vertical error in the magnitude of several tens of meters. To account for the vertical error of the UAV onboard GPS, we used relative elevations in our DEMs, assigning 0 m elevation to the respectively lowest point in the model. With the large amounts of photographs taken, the horizontal error is further reduced during processing, so that we cannot identify an offset in the horizontal orientation of our models compared to e.g. the lower resolved 20 TanDEM-X WorldDEM™. Therefor we assume a horizontal error in the range of our onboard GPS receiver, as previously explained. To avoid misinterpretations stemming from model tilt, possible slopes and artefacts have been ruled out by comparison with TanDEM-X WorldDEM™ data. To more carefully review the point clouds created during the process and to take first test measurements, we used the Compass plugin (Thiele et al., 2017) for CloudCompare. 25 Ground truthing and quality control of the resulting DEMs was acquired by including reference objects of known dimensions in our models. Our model accuracy lies within the same error range as in Bonali et al., (2019b), who achieved accuracies between 0.04 m -0.07 m for 50 m -100 m flight altitude horizontally and up to ca. 20 cm vertically, thus below and within the same order of magnitude as our accuracy.

Data extraction
The DEMs and ortho-mosaics were imported in Esri ArcMap to map the fractures manually as polylines on a scale of 1:100. The chosen mapping scale represents the middle ground between a more time-consuming, high-accuracy or a fast, but low-accuracy mapping scale, still enabling us to accurately map details of few tens of centimeters such as the edges of cooling columns. With a manual mapping accuracy of a few mm using GIS for the interpretation 5 in the DEM and ortho-mosaic at a 1:100 scale, the mapping error of the interpretation is in the same order of magnitude as our spatial resolution. We traced all observable fractures along their surface expressions, including the surface traces of faults, the foot-and hanging wall of dilatant faults as well as the adjacent fracture traces of dilatant joints. In a later step, the polylines were merged to polygons, representing the opening of the fracture on the surface. The strike directions of the fractures were measured as represented by a straight line from tip to tip. 10

Opening-width and vertical offset
To extract the opening-widths, parallel scanlines with 1 m spacing were created orthogonally to the average strike direction. Small deviations in strike were found to have no significant impact on the results, as the relative error is below 1% for deviations < 60° (Supplement S2). The scanlines were subsequently clipped with the polygons of the faults to measure opening-width (e.g. S6). The overall dilatancy of a fault is defined as the summed opening-15 width of subparallel faults and fissures.
Vertical offset was measured on the same scanlines along which the opening-width was determined. One measurement point was created along the scanline on the hanging wall and on the footwall (S6). The vertical offset was calculated from the elevation difference and linked to the opening-width measurement and the associated scanline. To prevent errors caused by local variations of surface elevation associated with the damage zone, the 20 elevation data were extracted a few meters away from the fracture, if possible, without crossing another fracture.
In this process, we also tried to avoid areas of rough surface with wavelengths of a few meters (large vegetation or lava blocks). To extract geometries of fractures in single models, these methods of opening-width and vertical offset extraction are viable. To compare absolute elevations in several DEMs, however, the use of GCP and the correction of the elevations is suggested. 25

Ground truthing and field observations
Ground based observations were used to complement the airborne datasets and interpretations. Vertical fracture walls, which are not as well resolved as horizontal surfaces in the top-down UAV photographs, mostly consist of successive lava flows (Fig. 3, supplement S5). The thickness of the lava flows varies between few cm to several meters. Contacts of the lava flows can be smooth transitions or sharp edges, locally with remains of volcanic glass, 5 mm to a few cm thick. The different morphology of the contacts is interpreted to influence the cohesion between the layers. Patterns of cooling joints may vary strongly between the layers, resulting in columns of few cm to several meters in length and diameter.
The correlation of piercing points of adjacent fracture walls as presented in Bonali et al., (2019b) is more accurate for fissures with small opening-width with respect to larger opening-widths. More correlation points can be found 10 directly on the fracture walls. Thus, we experienced the identification of reliable piercing points to be increasingly challenging with increasing opening-width and vertical offset. We attribute this to the observation of rubble and disintegrated columns within larger fractures, indicating the advanced erosion of the fracture wall. This is supported by the observations that small fissures may remain unfilled to large depths, while faults are filled with broken columns, typically a few meters under the hanging wall cutoff (Fig. 4b, d). Maximum depths for accessible 15 (larger) fractures and faults show peaks at 30 m up to 50 m in Thingvellir. The true depth of the opening is not observable, as the voids are filled with rubble, disintegrated columns and sediment. The filling of the voids can reach up to the fracture edges at the surface, also covering these ( Fig. 3a). This leads to the observation of a single scarp with no distinguishable opening from the elevated perspective of an UAV.
The apparent opening-width of fissures at the surface can be strongly misestimated when not directly measured at 20 the rock surface: erosion of soil-cover into the fracture leads to funnel shape on the surface, causing overestimation of the opening as also described in Bonali et al., (2019b). We further observed the reduction of visible openingwidth by vegetation that i) is large enough to prevent a clear line of sight on the fracture when seen from above and ii) moss and lichen patches that grow over the edges, thus reducing the visible opening (Fig. 5). As the thickness of the moss patches observed can vary between centimeters to decimeters, the relative error becomes 25 larger for overall small opening-widths, eventually covering the whole opening on the surface. For opening-widths summed over several small fractures, e.g. en-échelon arrangements of early stages of oblique faults, this error will sum up, leading to a large underestimation of the true opening-width at the surface.

Digital elevation models
Here, we focus on several representative DEMs of our dataset that were used to map the faults in detail and to quantify their geometries. A direct comparison between the resolution of our high-resolution DEMs and the ones provided by TanDEM-X WorldDEM™ is provided in Figure 6. To simplify and visualize the fracture geometries, 5 measurements of opening-width and vertical offset were sorted according to position along strike and plotted as x-y scatterplot, with the position along strike in meters on the x-axis and the extracted value, (opening-width, vertical offset) on the y-axis, scaling with the DEM above ( Fig. 7-11).

Asbyrgi
The field area 'Asbyrgi' is a fracture accumulation in the center of a large graben in the Kelduhverfi area ( fractures with vertical offsets and opening-widths below 2 m. They are underlapping a structure that can be traced from ca. 450 m to 1400 m, consisting of at least 3 segments of similar size, from i) 450 -700 m, ii) left stepping from 700 -1100 m and iii) right stepping from 1100 m -1400 m. The opening-width of the larger structure reaches its maximum values around 10 -12 m close to its center between 700 and 900 m along strike (al.st.), and decreases towards the tips, with a steeper gradient towards the south than towards the north. The vertical offsets 20 are < 2 m in i) and increase to 3.5 m in ii), where they reach a local minimum in the center, but again increase towards segment iii), in which the vertical offset remains between 2 and 4 m.

Theistareykir South
The DEM of Theistareykir South (

Thingvellir (Almannagjá) 25
With almost 7 km length, on average 200 m width and a relative elevation difference up to 53 m, the DEM of the NNE-SSW striking Almannagjá fault in Thingvellir is one of our largest high-resolution data sets (Fig. 10).
Bounded by Lake Thingvallavatn in the east, the DEM includes the western main fault of the postglacial graben. Measurements of opening-width, however, could be performed based on the ortho-mosaic.

Vogar
This DEM covers the adjacent shoulders of a graben associated with the Vogar fissure swarm (Fig. 11). Our focus while capturing the drone photographs was on the two NE-SW striking main faults of the graben. We were able to connect the two photosets by including several traverses orthogonal to the fault strike. Thus, we were also able to 15 cover several smaller fractures in the graben center. The maximum extent of the resulting DEM is more than 2 km

General observations
Comparing the four field areas and the observed fracture geometries therein, no major differences in structure are observable. Fractures and faults in all areas can be described according to opening-width, filled or covered 5 openings, vertical offset and associated structures such as tilted blocks (Fig. 12). In each of the field areas we can define the following endmember structures associated to MDF from field observations and insights from our DEMs: 10 1) Dilatant Faults Dilatant faults are opening mode fractures with vertical walls and measurable opening-width and vertical offset (Fig. 12a). The opening can be filled with tilted columns, broken columns, sediment or younger lava flows and vegetation. Basalt columns can be jammed between the fracture walls or form part of the highly porous aggregate between the walls (Fig. 4d). Filling with sediments on top of the rubble can cover the gap completely (Fig. 12c), 15 so that the opening-width can no longer be determined at the surface (Fig. 3a) and DEM (Fig. 14). This has also been described by Trippanera et al. (2015) as their type A structure. The complete filling of the opening causes the faults to appear as a single scarp on the surface. In basalts, the exposed fracture is sub-vertical and follows the geometry of the basalt columns.
The faults may develop dilatancy of up to 15 m accompanied by vertical offsets in the same order of magnitude. 20 The void between the walls may remain open down to depths of 20 -30 m, but is usually shallower (Fig. 1). The faces of the fracture walls expose successive layers of lava flows with mostly vertical cooling joints (cf. supplement S5). Field examination of the walls on both sides usually do not allow to match patterns of columns, indicating that material between the walls is missing due to erosion. In the DEMs (Fig. 15) MDF have measurable openingwidths, (fracture walls are visible) with clear vertical offsets (Fig. 12a). All faults, including the ones that appear 25 non-dilatant, can be shown to have an opening. At apparently non-dilatant faults the opening is covered by sediment so that it is not easily visible in the field or DEM (Fig. 14).

2) Fissures
Fissures are opening mode fractures with a prominent surface aperture and no significant vertical offset (Fig. 12b).
Their local geometry is governed by pre-existing cooling joints in the basaltic lava flows (Grant and Kattenhorn, 2004;Holland et al., 2006) and thus develop in dm-scale sawtooth patterns along the boundaries of the basalt columns (Fig. 3a). Fissures can be early stages of MDF or represent the lateral ends of MDF, which show enéchelon fracturing when formed at oblique slip (Acocella et al., 2003;Grant and Kattenhorn, 2004;Gudmundsson, 5 1987a;von Hagke et al., 2019). Earlier studies defined fissures as fractures with vertical offsets < 1 m (Grant and Kattenhorn, 2004;Gudmundsson, 1987b). However, our data (Fig. 15) show a more substantial surface roughness, leading us to define a maximum allowed vertical offset of 2 m (cf. supplement S1). Fissures tend to cluster around fault tips and occur in zones parallel to faults, and are rarely longer than 100 m.

3) Tilted Blocks
Tilted blocks (TB), also referred to as monoclines (Grant and Kattenhorn, 2004;Martel and Langley, 2006;Smart and Ferrill, 2018;Sonnette et al., 2010) can develop in the hanging wall of dilatant faults, creating a surface dipping away from the footwall. The length of single tilted blocks ranges from several meters up to several hundred meters; widths range between several meters to tens of meters, and depending on their subsurface geometry and kinematics 15 three different types of TB can be distinguished . An example of a Type III (Kettermann, 2019) TB is provided in Figure 3d. A sketch showing the expected behavior of opening-width and vertical offset relationships of a Type I TB is provided in Figure 12d. In the DEMs we identify TB quantitatively based on the slope on the hanging wall, dipping away from the footwall (Fig. 16), following Kettermann et al., (2019). Dips of TB are commonly few degrees; exact measurements may be perturbed by vegetation cover and surface roughness 20 of the lava flows. For a more detailed analysis of TB including kinematic models the reader is referred to Kettermann et al., (2019), here we focus on the presence or absence of TB along the faults.

Interpretation of the mapped DEM data
We mapped the fractures as previously described and measured opening-width and vertical offset along their strike 25 in the DEMs. Surface structures were identified either in the high-resolution models or, for larger structures, with TanDEM-X WorldDEM™. Based on the definitions above, each measurement in the database was assigned to one of the proposed endmember types, i.e. Fissures, Tilted Blocks (TB) or Dilatant Faults, which are subdivided in dilatant faults with discernible opening (DF) and no discernible opening. 30 Asbyrgi (Fig. 17a): Due to the small vertical offsets and en-échelon arrangement, the fractures in the south are identified as mode I fissures. The larger segmented structure towards the north shows different endmembers: i) is a fissure, followed by two segments of TB, as their aspect ratio is biased towards larger opening-widths in relation to the vertical offset, further aided by the surface gradient of the hanging wall, which is dipping away from the footwall. As the 5 high-resolution DEM coverage is restricted in these directions, TanDEM-X WorldDEM™ elevations were used to complement the data.
Krafla N (Fig. 17b): 10 The section of the MDF shows DF and several TBs. Their direct influence on the measured opening-width is apparent in the plot, as the opening-width strongly increases along the TB, and decreases again in sections without Thingvellir, Almannagjá (Fig. 19) : 30 As described by Gudmundsson (1987b), the Thingvellir fissure swarm consists of two types of fractures: extension fractures with vertical offsets < 0.5 m, in which category most fissures fall, and dilatant normal faults, that often turn into fissures at their tips. Reviewing our DEM, the fault parallel structures along the hanging wall, we interpret the eastwards dipping slope as several TBs, in accordance with the interpretation of Kettermann et al. (2019). In combination with the large opening-width and vertical offset, the proposed endmember type is TB. The en-échelon fractures in the southern part of the DEM show vertical offsets < 1 m, when measured on the adjacent sides of single fractures, thus qualify as fissures. However, the eastwards sloping surface is very prominent at this location, 5 resulting in a combined vertical offset of the faults locally exceeding 10 m.
The en-échelon fractures have been interpreted as surface expression of an old weakness, possibly a Pleistocene fault that has been covered by postglacial lava flows (Gudmundsson, 1987b). Thus, the opening-width and vertical offsets are interpreted to represent the geometry of a larger, underlying structure. We classify this structure as a very large TB. The overall trends of opening-width as well as vertical offset are maximal at the center of the fault 10 and decrease towards the tips. This is interpreted as the typical behavior of an ellipsoidal fault. The dataset of the Almannagjá fault is not only the largest of our datasets with almost 7 km length, but also includes the largest values of opening-width, vertical offset and the most prominent TBs. We infer that these large offsets are possible due to the relatively larger dilation rate of Thingvellir fissure swarm as compared to e.g. the proximate Vogar fissure swarm (Gudmundsson, 1987b). 15

Discussion
In this study, we used UAV based photogrammetry to create high resolution DEMs of representative faults and fractures of the Iceland Rift. We show that these DEMs can be used to map faults and fractures in much more detail when compared to aerial photography or satellite imagery. Furthermore, these DEMs can be used to extract 20 geometries much faster than taking measurements in field. Resolutions of 5-15 cm/pixel are appropriate to map fractures in volcanic settings, corroborating the findings of Bonali et al., (2019b). We derived three endmember types of surface expressions of fractures linked to the massively dilatant faults in Iceland, based on the ratio of opening-width and vertical offset, similar to earlier studies on dilatant fractures (Tentler and Mazzoli, 2005;Trippanera et al., 2015, and references therein). 25 We visualized this ratio for different field areas, color coded for different endmembers (Fig. 20 and Fig. 21). The straight lines plotted in the figure represent the relationship between vertical offset and horizontal opening corresponding to a simple dilatant normal fault at depth with dips ranging between 60°and 70°, being within the dip range as commonly inferred for Iceland (Angelier et al., 1997;Grant and Kattenhorn, 2004;Gudmundsson, 30 1992;von Hagke et al., 2019;Trippanera et al., 2015). We further define the parameter R as ratio of opening-width and vertical offset as Eq. (1): where O is the opening-width and V the vertical offset measured on one scanline. We calculated R for all 5 endmember types with measured opening-width sorted by area (Fig. 22). Expected R values for fault dips between 60° and 70° are within the interval of R(60°) = 0.58 and R(70°) = 0.36.
The fractures in Asbyrgi have been classified as fissures and tilted blocks. The proximity of the fissures to the fault and the prominent surface slope may possibly be interpreted as doming of the surface related to a subsurface dike 10 intrusion, as described for similar structures by Tentler and Temperley (2007). However, in this model, a slope on both sides of the fractures is expected. Since in Asbyrgi the local slope can only be identified at the hanging wall, we interpret the structure as TB. The plot of our data from Asbyrgi does not show a clear separation between the point clouds for fissures and TB (Fig. 20). The different endmember classification in the plot is a result of our definition for the cutoff of vertical offset, which automatically defines fractures with vertical offsets < 2 m as 15 fissures, thus creating an artificial threshold in a smooth transition within a continuum. The data show that TB and fissures are endmembers of a continuum, without a gap or separation in the data. Most measurements from Asbyrgi show much higher R values as would be expected for slip along a basement fault dipping steeper than 60° (Fig.   12a, d), in agreement with the geometry of TBs  and opening mode fissures.
In the Asbyrgi area we noted TB data points with very low opening-widths. These measurements are in areas 20 where the visible opening is reduced by vegetation and soil cover. This effect can be seen mostly in relay zones where the opening-width is summed over smaller fractures or at the fault tip in the north, thus indicating an underestimation of opening-width by mapping errors, depending on the type and coverage of vegetation, as explained in 2.3.2 and Figure 5. A general trend of the surface dipping at approximately 5% from W to E can be identified in our high resolution and the TanDEM-X WorldDEM™ data. Since the vertical offsets have been 25 measured at a distance to the fracture, an additional error is made, leading to the overestimation of vertical offsets and thus underestimation of R.
Krafla N includes the two endmember types of MDF: DF and TB, which form two clusters (Fig. 20). DF with the lowest R ≈ 0.49 are interpreted to be the closest to "classic MDF" (Fig. 12a), but most of the data show much 30 larger R with TB often having R > 1 (Fig. 20). However, the two clusters of R in Krafla overlap. R in data of Theistareykir S show a similar trend, where DF are within the commonly assumed interval of R with R ≈ 0.47, while TB tend towards larger R with R ≈ 1.07. However, in Theistareykir S, the vertical offsets are underestimated at the TBs due to the missing horizontal part of the hanging wall in the DEM, leading to a too high R-value.
Assuming the same vertical offset as in the DF part of the fault, the vertical offset of the TBs might be underestimated by ca. 5 m. Overlaps of DF and TB can be explained by the resolution of surface tilt; TB with 5 surface dips << 5° may be unrecognized and misinterpreted as DF. Particularly the transition from TB to DF is prone to interpretation errors when no clear boundary is visible, e.g. due to a smooth transition from tilted to horizontal hanging wall.
The tendency of TB towards large opening-width with respectively smaller vertical offset becomes even more 10 apparent in Thingvellir, where TB are well developed and the majority of the measurements have R > 1 (median R ≈ 1.59). The smaller, detached cluster between 0 -5 m vertical offset and 5 -20 m opening-width (Fig. 20) results from the measurements taken on the en-échelon fractures in the southern part of the main fault. These classify as fissures when viewed as single fractures (cf. supplement S3), but when counted cumulatively, reflect the underlying fault structure as a TB. This strain partitioning has been described as separate fault structure by 15 Trippanera et al. (2015).
The fractures mapped in Vogar are of particular value, as all proposed endmember types associated to MDF are present: fissures form a continuum with DF in the north (Fig. 20). Faults with no distinguishable opening-width were all confirmed in the field to be filled with sediment at the northern and southern fracture. R-values of DF in the north (R ≈ 0.65) and south (R ≈ 0.62) are similar to R values of TB (R ≈ 0.67), while TB trend towards larger 20 vertical offsets (Fig. 20). The high value of R for DF is most likely caused by strong erosion of the fracture walls, leading to overestimation of the opening-width, and accumulation of material in the opening, up to completely covering it. With our selected areas all types of TB as defined in Kettermann et al. (2019) are covered: Type I in Vogar and Krafla, Type II in Theistareykir and Type III in Thingvellir, all resulting in R values larger than expected for the "classic" MDF. Solely relying on measurements of opening-width and vertical offset, expressed as ratio R, 25 different types of TBs cannot be distinguished. This is because vertical offset is less influenced by surface structures, but deep processes and fault geometries, when compared to opening-width. Furthermore, vertical offset has been measured outside the influence area of TB, while opening-width depends on TB geometry and erosion. Therefore, simple numeric measurements are not sufficient to identify the type of TB. Additionally, opening-width can vary strongly over short distances. Generally, in all areas the data of the endmember classifications form one 30 continuous cloud.
We added further measurements from Asbyrgi and Krafla (cf. supplement S4) in combined plots (Fig. 21) including data from Iceland and Ethiopia (Tentler and Mazzoli, 2005;Trippanera et al., 2015, and references therein). The distribution of the endmember types remains consistent. However the data of Tentler and Mazzoli, (2005) Acocella et al. (2003) and Gudmundsson (1992). Measurements classified as DF that show higher values of R (associated with shallower basement fault dips) are the result of erosion and the disintegration of the fracture walls, leading to an overestimation of the opening-width.

iii)
TB plot in clusters that also overlap with DF, but mainly have R > 0.5, depending on fault geometry.
Faults with large vertical offsets in relation to the opening-width can produce TB clusters similar to those 20 of DF, as observed in Vogar. However, when compared to their non-TB counterparts, the TB clusters trend towards larger opening-widths (cf. supplement S4&5). This trend can be explained by the rotation of the hanging wall away from the footwall and a resulting increased aperture .
Therefore, measurements of opening-width on TBs lead to an overestimation of the overall dilatancy of the fractures. 25 Vertically elongated clusters in the plots (Fig. 21) are the result of relatively stable vertical offsets with smooth gradients over long distances along strike at most faults, while the opening-width shows much more local variations. Piercing point correlations of adjacent fracture walls are more reliable for small openings and openingwidth measurements can vary strongly over the fault length. Consequently, errors increase towards the fault center where displacement is largest, as shown for Asbyrgi, Krafla and Theistareykir. Maximum opening-widths are 30 larger than 60 m in the Thingvellir dataset with a maximum vertical offset of 40 m. Our data is consistent with the measurements of Tentler and Mazzoli (2005) and Trippanera et al. (2015), who have however not studied faults with large values as present in Thingvellir. Figure 22 shows the distributions of R for the different endmembers, but the overlapping distributions indicate that all mapped structures belong to a larger continuum with smooth transitions between endmembers. Furthermore, vertical offsets with small opening-widths are rare in our and the data of Tentler and Mazzoli (2005) and Trippanera et al. (2015), resulting in a gap between MDF with and without TB, and DF with no discernible opening (Fig. 21). The reasons for this gap in our data is i) the mapping procedure 5 of the fracture traces (cf. supplement S6) at the transition between DF and no opening, and ii) the interval length between scanlines (cf. supplement S6). The transition between opening and non-opening occurs over shorter lengths, therefore the transition between opening and the fully filled state with decreasing opening-widths is not fully covered by the scanlines, since a smooth transition between the filled and open state will be missed when it appears over a distance smaller than our scanline interval of 1 m. 10

Data availability
A table containing the measurements of opening-width and vertical offset and additional figures are provided in the supplement.

Author contributions
Christopher Weismüller acquired the data, did the processing and data extraction and prepared the manuscript with 5 contribution from the co-authors. Janos L. Urai contributed to the conceptualization, discussion and interpretation of the data and acquired funding. Michael Kettermann acquired data and contributed to the processing, conceptualization and discussion. Christoph von Hagke and Klaus Reicherter contributed to the conceptualization and discussion of the manuscript and acquired funding.

Competing interests 10
The authors declare that they have no conflict of interest.

Acknowledgements
We would like to thank the park administration of the Thingvellir National Park for giving us the permission to use our drone in the park outside of business hours and Daniele Trippanera for sharing his data with us. Further thanks go to the Deutsche Forschungsgemeinschaft for funding this project. Project number 316167043. Project 15 name: MDF: The structure and evolution of near-surface massively dilatant faults. TanDEM-X WorldDEM™ data is provided by a DLR Science grant, 2017. We further thank two anonymous reviewers who helped to improve this manuscript with their comments.