Rifting Kinematics Produced by Magmatic and Tectonic Stresses in the North Volcanic Zone of Iceland

In the North Volcanic Zone of Iceland, we studied with the greatest possible detail the complete structural architecture and kinematics of the whole Theistareykir Fissure Swarm (ThFS), an N-S-trending, 70 km long active rift. We made about 7500 measurements along 6124 post-Late Glacial Maximum (LGM) extension fractures and faults, and 685 pre-LGM structures. We have collected the data over the last 6 years, through extensive field surveys and with the aid of drone mapping with centimetric resolution. In the southern sector of the study area, extension fractures and faults strike mainly N10°-20°, the opening direction is about N110°, and the dilation amount is in the range 0.1–10 m. In the central sector, faults and extension fractures strike mainly N00-10°, the opening direction is N90-100°, and the dilation amount is 0.1–9 m. In the northern sector, extension fractures and faults strike N30-40°, the opening direction is about N125°, and the dilation amount is 0.1–8 m. The variations in strike are attributable to two processes: the interaction with the WNW-ESE-striking Husavik-Flatey transform fault and Grímsey Oblique Rift (Grímsey lineament), and the structural inheritance of older NNE- to NE-striking normal faults. Most extension fractures show a minor strike-slip component: a systematic right-lateral component can be accounted for by the interaction with the WNW-ESE-striking fault zones and the regional, oblique opening of the rift. We regard dyke propagation as a possible cause for the more complex strike-slip components measured at several other fractures. Cumulated dilation and fracture frequency decrease along the rift with distance away from the Theistareykir volcano, situated in the central sector of the ThFS. This is interpreted as a decrease in the number of dykes that are capable of reaching great distances after being injected from the magma chamber.


INTRODUCTION
Rifting is a widespread process in both submarine and terrestrial environments, the study of which enables enhancing knowledge about how plate separation works and about related geohazards such as volcanism and seismicity. Although classical rifting is assumed to take place as plate separation along a direction perpendicular to plate boundary, frequently this process is complicated by an oblique component that results in shear along the rift zone (Withjack and Jamison, 1986;Tron and Brun, 1991). Oblique opening brings about a complex network of normal as well as strike-slip faults and transtensional structures, whose kinematics and orientation are guided by the angle between the direction of relative motion of the two diverging plates and the rift trend (Corti, 2012). This complexity is made even more so by the influence of magma-induced stress that produces deformation at the surface, especially in relation to shallow dyke emplacement (Mastin and Pollard, 1988;Rubin and Pollard, 1988;Rubin, 1992;Gudmundsson et al., 2008;Tibaldi, 2015). Recently, the development of strike-slip motions has been described as taking place along fractures induced by shallow dyking at the 2014-2015 Bárðarbunga eruption in southern Iceland (Ágústsdóttir et al., 2016;Ruch et al., 2016). Consequently, it seems that both tectonic and magma-induced stresses are capable of inducing strike-slip components along rift structures. This setting is made even more complicated by the fact that rift zones with rift-parallel shear components can result in different fault patterns: on one side, oblique rifting typically produces en-échelon fault systems with segments not perpendicular to the direction of plate separation (Fournier and Petit, 2007; Figure 1A). On the other side, Karson (2017) proposed the partitioning of motions with strike-slip faults at the sides of a propagating rift, within which normal faults do occur ( Figure 1B). In the latter case, the presence of strike-slip faults does not necessarily result from oblique plate separation, but it can be the effect of the outward rift propagation from a hotspot center.
In order to shed light on the possible causes of shear components at rifts, it is necessary to carry out detailed surveys of fracture location, strike, age, and kinematics. While at continental rifts the surface rupture processes allow for a direct observation of fracture systems, in oceanic rifting settings the processes mostly occur at great depth below sea level, thus hindering detailed data collection. The North Volcanic Zone (NVZ) of Iceland represents the surface expression of the Mid-Atlantic Ridge, where plate spreading occurs at about 2 cm/yr (Geirsson et al., 2006;DeMets et al., 2010). Along the Icelandic rift system it is possible to study ongoing processes of seismic deformation and volcanism; this geodynamic setting provides a formidable natural laboratory to investigate the complex interactions between magma-and tectonic-induced stresses (e.g., Einarsson and Brandsdóttir, 1980;Björnsson, 1985;Heki et al., 1993;Jonsson et al., 1997;Cattin et al., 2005;Wright et al., 2012;Tibaldi et al., 2013;Pasquaré Mariotto et al., 2020). In the NVZ, previous authors suggested the presence of shear components along the Theistareykir Fissure Swarm (ThFS) (Karson, 2017;Tibaldi et al., 2019;Bonali et al., 2019a,b) and interaction with transform faults (Tibaldi et al., 2016a(Tibaldi et al., ,b, 2018; this makes the ThFS an outstanding site for gaining major insights into deformation patterns at complex diverging plate boundaries. With the purpose of contributing to a better understanding of the above, we hereby provide a complete picture of the ThFS, composed of thousands of extension fractures and faults extending for about 70 km (Figure 2). The ThFS represents the westernmost rift of the NVZ. This has been active since 8-9 Ma (Saemundsson, 1974;Young et al., 1985; FIGURE 1 | Sketches of rift zones with rift-parallel shear component that takes place through different fault patterns: (A) oblique rifting with en-échelon faults, the strike of which is not perpendicular to plate separation direction (modified after Fournier and Petit, 2007). (B) Partitioning of motions with strike-slip faults at the sides of a rift that is propagating outward from a hotspot center (modified after Karson, 2017).
Frontiers in Earth Science | www.frontiersin.org FIGURE 2 | Tectonic setting of northeastern Iceland. Within the North Volcanic Zone, the Mid-Atlantic Ridge is offset by the Húsavık-Flatey Fault (HFF) and by the Grimsey Lineament (GRL), which, together with the Dalvík Lineament (DAL), compose the Tjörnes Fracture Zone. Yellow stars locate historical earthquakes with M > 6 (after Stefansson et al., 2008). Orange stripes represent volcano-tectonic rift zones, also termed "fissure swarms," where the main normal faults, extension fractures and volcanic fissures are shown (after Magnúsdóttir and Brandsdóttir, 2011;Hjartardóttir et al., 2012; Quaternary central volcanoes are highlighted as white triangles. Box locates ( Figure 4A). The inset shows the location of the area and the main volcano-tectonic rift zones of Iceland (after Einarsson and Saemundsson, 1987;Hjartardóttir et al., 2016a) Bergerat and Angelier, 2008) and is suggested to have been propagating northwards during the past million years (Gudmundsson, 2007). The goal of our work has been to define, with an unprecedented detail, the kinematics of a rift, by performing systematic field and drone surveys (Figure 3) aimed at assessing the geometry, direction and amount of opening of faults as well as extension fractures. Measuring several sites along each structure has allowed us to collect a great deal of data (about 7500 measures), instrumental for defining strike-slip components and distinguishing the different roles played by tectonic and magma-induced forces in affecting post-Late Glacial Maximum (LGM) deformation within the ThFS. These data are meant to integrate the first investigations made along the same rift by Bonali et al. (2019b), who focused only on its northern part. We are now able to present new data also from the central and southern parts of the rift; in detail, we focused on: (i) the geometric characterization of all the structures, with special reference to the strike and length of faults and fractures; (ii) the quantification of the opening directions of extension fractures and the possible presence of shear components; (iii) the estimation of the strain amount along the rift; (iv) the assessment of the origin of the structures; (v) the evaluation of the overall working of the rift.

GEOLOGICAL BACKGROUND
The NVZ is a 200 km long volcano-tectonic structure extending from the central portion of the Icelandic hotspot as far the Tjörnes Fracture Zone to the north (Figure 2). The latter is a zone of shear divided into three main, N120 • -striking segments that are, from north to south: the Grímsey Oblique Rift (also known as Grímsey Lineament), the Húsavík-Flatey Fault and the Dalvík Lineament. All these three zones are seismically active and connect with the Kolbeinsey Ridge (KOR) through a complex oceanic transtensional rifting zone (Einarsson, 1991;Magnúsdóttir and Brandsdóttir, 2011).
The NVZ is composed of five, about N-S-trending, volcanic systems: the ThFS, Krafla, Fremrinámar, Askja and Kverkfjöll . Each consists of 5-20 km-wide and 60-100 km-long fracture swarms and a central volcano (Saemundsson, 1974). Other features are normal faults, eruptive fissures and extension fractures that strike parallel to the rift (Figure 2). The ThFS is an active, N-S trending, 70 km long volcano-tectonic rift that forms the westernmost edge of the NVZ, marked by the presence of normal faults and extension fractures, parallel to sub-parallel to each other, striking mostly N00-30 • (Angelier et al., 1997;Acocella et al., 2000;Dauteuil et al., 2001;Hjartardóttir et al., 2012Hjartardóttir et al., , 2016aPasquaré Mariotto et al., 2015). The ThFS is also dotted with a few volcanic vents, but their number is too limited to enable any statistical analysis. The rift is also home to the major Theistareykir central volcano. Along this rift, some structural changes can be noticed, such as changes in strike, most of which occur in the sector where the ThFS intersects the Husavik-Flatey Fault, the latter characterized by dominant right-lateral strikeslip movements (Gudmundsson et al., 1993;Fjäder et al., 1994;Magnúsdóttir and Brandsdóttir, 2011;Pasquaré Mariotto et al., 2015;Tibaldi et al., 2016a,b).
Along the ThFS, several depositional units crop out, represented by sedimentary sequences from Miocene to Pleistocene times, as well as lava flows of Pliocene, Pleistocene and Holocene age (Saemundsson, 1974;Garcia et al., 2002). The latest eruption here took place about 2.4 ka BP, with the emission of the "Theistareykjahraun" lava flows between the central volcano and the Husavik-Flatey Fault (Saemundsson et al., 2012; Figure 4A). Active vertical deformation at Theistareykir volcano, measured through GPS-based methods (Metzger et al., 2011), was interpreted as induced by inflation of the underlying magma chamber. During historical times, the following rifting and volcanic unrest episodes are documented: (i) between 1867 and 1868, an M 6 earthquake hit in the northern part of the ThFS, accompanied by volcanic eruptions in the northern offshore sector of the rift, within the Mánáreyjar volcanic system (Figure 2; Thoroddsen, 1925;Halldórsson, 2005;Magnúsdóttir and Brandsdóttir, 2011); (ii) between 1884 and 1885, an earthquake swarm culminated in an M 6.3 earthquake along the rift, near the coast to the north (Halldórsson, 2005); (iii) between 2007 and 2009 an 8 cm uplift in correspondence of an area centered on the Theistareykir volcano was recorded, and it was interpreted as the result of an inflation of the magma chamber below (Metzger et al., 2011;Metzger and Jónsson, 2014).

MATERIALS AND METHODS
All the structures documented in the present work were first mapped on aerial photos, and then they were studied and classified during fieldwork. In order to collect data along the identified structures, both field surveys and Unmanned Aerial Vehicles (UAV) surveys were employed, the latter integrated with the aid of the Structure-from-Motion technique, enabling to produce Orthomosaics and DSMs with centimetric resolution. The structures investigated in the field were classified as: (i) normal faults, when a systematic and continuous vertical displacement > 50 cm is observed (Figures 3A,B) and (ii) extension fractures, which are dilatant fissures with a vertical component <50 cm ( Figure 3C). In light of the above, the term "fracture" is used in the present work to represent the overall population composed of faults (with a vertical component > 50 cm) and extension fractures, the latter including structures with pure extensional opening ( Figure 3D) or with a minor strike-slip component ( Figure 3E).
In the field, wherever possible, we verified that most faults have pure dip-slip motions or a very minor strike-slip component. We also observed the presence of fissures in correspondence of the footwall block at normal faults. Wherever the fissures were located close to the associated faults, they were not included in the calculation of the opening amount; in fact, they can be regarded as due to gravitation instability induced by unbuttressing along the fault scarp. These fissures are usually quite short and can be classified as belonging to the fault structure.

UAV Use and Structure From Motion
Some areas along the rift were surveyed by means of UAVs; in particular, we used three different types of quadcopter, which is the most used and versatile type: the DJI Phantom 4, the DJI Phantom 4 PRO and the DJI Spark. The DJI Phantom 4 and the DJI Phantom 4 PRO are both equipped with a high-quality camera together with GPS coordinates, provided by the integrated Satellite Positioning Systems GPS/GLONASS (referred to the WGS84 datum). Photos are taken at 12.4 and 20 Megapixels, respectively. The UAVs can fly for a long time (28-30 min approximately for each battery). Thanks to a built-in, high precision 3-axis gimbal, flights are smoother and the collection of the images is more stable.
We also made use of the smaller DJI Spark for surveys over more limited areas. This quadcopter is equipped with a camera taking photos at 12 Megapixels, and is capable of flying for as long as 16 min.
The UAV-captured photos have been taken in such a way as to attain a 90-85% overlap along the path and a 80-75% one in a lateral direction: this is crucial to obtain a good alignment of the images and to reduce the distortions on the resulting Orthomosaics. With the DJI Phantom 4 (and the PRO version), we collected pictures between 70 and 100 m of height from the ground. With the DJI Spark, pictures were collected between 20 and 30 m of height. The collected images were then processed through Agisoft Metashape, a commercial Structurefrom-Motion (SfM) software.
In order to allow for the co-registration of datasets and the calibration of the models (orthomosaics and DSMs) resulting from SfM processing, World Geodetic System (WGS84) coordinates of four artificial Ground Control Points (G) were fixed near each corner within each surveyed area; one was fixed in the central part (e.g., James and Robson, 2012;Turner et al., 2012;Westoby et al., 2012), with the purpose of reducing the "doming" effect resulting from SfM processing (e.g., James et al., 2017). The GCPs were surveyed by way of low-cost single frequency receivers in RTK configuration (with centimeter-level accuracy). The SfM technique enabled us to identify matching features in different images, collected along a defined flight path, and to combine them to create a sparse and dense cloud, a mesh, an Orthomosaic, a DTM/DSM (further details in Stal et al., 2012;Westoby et al., 2012). Particularly, our Orthomosaics and DSMs are characterized by a resolution in the following ranges: 1.0-4.5 and 2.0-18.0 cm/pixel, respectively.

Data Collection
For each structure, both normal faults and extension fractures, we assessed the overall strike and length in a GIS environment, where we plotted all the measured structures, mapped and classified as explained above. Other measurements, which will be described below, were made only for structures longer than 20 m.
As far as faults are concerned, in the field we determined the dip, dip-slip and strike-slip offsets by systematically measuring displacements of piercing points such as stream beds, crests, lava flows and deposits. As the downdip slip is the dominant one, we measured the vertical offset in the field at intervals of 50 m, or continuously through drones, and we were able to come up with a fault slip profile, following Manighetti's approach (Manighetti et al., 1997(Manighetti et al., , 2001. In the field, vertical offsets were measured by tape and/or laser rangefinder, where fault scarps were in the order of meters, and by GPS measurement where the scarps were in the order of tens of meters; furthermore, we assessed offsets on DSMs by measuring the difference in altitude along a topographic profile traced orthogonally to the fault scarp. Regarding extension fractures, both in the field and on DSMs we measured strike, opening direction and net opening amount. In the field, we measured the opening direction with a compass, where corresponding piercing points on each side of the extension fracture could be located (Figures 3D,E); the net opening amount was measured along the opening direction, by means of a tape or a laser rangefinder. In case fractures were filled by sediments, we measured the amount of opening across their lower portions and only where fracture rock walls were detectable, so as to avoid any possible overestimation of the dilation amount due to gravity, frost weathering, or a combination of both, which can create a funnel-shaped profile ( Figure 3F). On the orthomosaics obtained through UAV surveys, the opening direction and net amount of opening were quantified by assessing the direction and measuring the length of a line connecting two piercing points on each side of a fracture (in a GIS environment). All these measurements were made only along extension fractures that are at least 50 m away from a scarp, so as to avoid local gravity effects. Moreover, opening amounts were also collected along N109 • -striking transects, each spaced 1 km in a N-S direction; this was done with the purpose of evaluating along-rift variations in the extensional strain parallel to the average direction of plate spreading affecting this part of Iceland (DeMets et al., 1994(DeMets et al., , 2010Metzger et al., 2011Metzger et al., , 2013Hjartardóttir et al., 2012;Drouin et al., 2017).

Pattern of Faults and Extension Fractures
In the study area, the structures are found in pre-LGM lithostratigraphic units and in post-LGM units ( Figure 4A). The 685 structures belonging to the first group, located in the western part of the area, strike mainly from N-S to NNE-SSW, and subordinately NE-SW and NW-SE (Figures 4B,C). They are represented by 679 normal faults, mostly dipping E and W, and subordinately SE to S, and 6 extension fractures, striking NNE-SSW and NW-SE. The majority of post-LGM structures (total = 6124) strike NNE-SSW, whereas minor sets strike N-S and NE-SW (Figures 4A,D). The same orientations can be recognized both analyzing extension fractures only (total = 4181) and faults only (total = 1943) (Figures 4E,F). The structures belonging to the Husavik Flatey Fault strike mostly NW-SE, and subordinately WNW-ESE and NNW-SSE ( Figure 4B). In the northern part of the study area, most normal faults have downthrown to the east, as exemplified in Transect 1 of Figure 4G. In the central part, most normal faults have downthrown in the opposite direction (Transect 17 in Figure 4G). South of the Theistareykir volcano, most normal faults dip to the east again (e.g., Transect 30 in Figure 4G), while further south both fault dips are developed.
The length of pre-LGM structures is shown in Figure 5A, resulting to be generally longer than the structures affecting post-LGM units ( Figure 5B)

Opening Directions
As can be observed in Figure 6A, the sites of measurement for a total of 1633 are distributed homogeneously in order to have a complete areal coverage of the whole rift zone. The resulting range of measures is from N 27 • to N 174 • , with an average of N 111 • and a standard deviation of 17. By analyzsing the distribution of opening directions with respect to Easting (Figure 6B), it is possible to observe the following pattern: in the western part of the rift, directions are mostly in the N80-140 • range, whereas they trend mostly between N100 • and N160 • in the eastern part. As regards the distribution of opening directions with respect to Northing (Figure 6C), the southern and central sections of the rift show similar ranges, mostly in the order of N80-140 • , whereas the northern part shows a range mostly of N100-160 • . Thus, combining the two graphs, it is possible to highlight a clockwise rotation of opening directions from the southwestern to the northeastern part of the ThFS.
The graph of opening directions with respect to extension fracture strikes shows that, as the latter rotate clockwise from NW-SE to NE-SW, so do opening directions, which rotate also clockwise from around N60 • to N160 • (Figure 6D). Although this indicates that the two are linked together, as will be discussed in the related section, it is worth underscoring the wide data scatter, which may indicate a more complex relation.  , 2010Hjartardóttir et al., 2012). Based on GPS data from different time windows (1997-2011, 2006-2010, 2008-2014), the short-term spreading vectors are N109-115-112 • (Metzger et al., 2011(Metzger et al., , 2013Drouin et al., 2017).

Strike-Slip Components
By plotting the opening direction azimuth with respect to the strike-slip component (Figure 7B), it is possible to note a shift from left-lateral to right-lateral component, along with the clockwise rotation of opening directions. The comparison of the strike-slip component with respect to Easting shows a slightly higher frequency of right-lateral components along the westernmost part of the rift (Figure 7C). The comparison of the strike-slip component with respect to Northing indicates the more frequent occurrence of right-lateral components northwards ( Figure 7D).

Extensional Strain
We measured the extensional strain by cumulating the contribution from normal faults ( Figure 8A) and extension fractures ( Figure 8B). Opening amounts were measured in the field as well as on UAV-derived high-resolution images across the study area, totaling 2745 measurements (1698 along extension fractures and 1047 along normal faults). In order to better assess extensional deformation, we cumulated extension values also along N109 • -striking paths, each spaced 1 km in a N-S direction, parallel to the average long-term spreading direction, obtained as explained above. At extension fractures, as already mentioned, we quantified the net opening. At normal faults, we calculated the horizontal extension based on the downdip slip and the dip of the fault plane, averaged at a value of 75 • (Forslund and Gudmundsson, 1991). Extension fracture opening is in the 0.1-10 m range, with an average of 1.2 m (Figure 9A). The higher values are in the south-central part of the rift. Horizontal dilation values along normal faults are in the 0.1-55 m range, with an average of 2.4 m ( Figure 9B).
In regard to the extensional strain measured along the N109 • transects, values have been calculated both considering only fractures affecting Post-LGM units and considering also fractures affecting Pre-LGM units. In the first case, the cumulated dilation is in the range of 3-6 m north of transect 11 (Figure 9C). A minimum is observed at transect 12, and then a rather steady increase can be noticed southward as far as transect Regarding the number of fractures calculated along the various transects, it ranges mostly around 1-20 in Post-LGM units ( Figure 9D). This range is also valid in the rift section where there is the greatest value of cumulated dilation (i.e., between transect 20 and 22). This means that in this area the same amount of fractures takes up a larger deformation than in the other sections of the rift. Considering also Pre-LGM units, the number of fractures increases showing a range between 4 and 28.
The average fracture spacing is 0.40 km, with a maximum value of 1.35 km if we consider only fractures affecting Post-LGM units, whereas it presents an average value of 0.67 km and a maximum value of 1.52 km considering also Pre-LGM units ( Figure 9E).
In Post-LGM units, the area affected by dilation is wider in correspondence of the central sector of the ThFS (transects 17-19) (Figure 9F). From here, it slightly increases northward from 0.9 km to 3 km, whereas it decreases southward from 4.5 to 2 km. In the southernmost sector of the rift, each transect (36-40) intersects only one fracture, setting to 0 km the average spacing between fractures and the area of deformation. Taking into account also the Pre-LGM units, the width of the area affected by dilation stays the same between transects 17-19, whereas it increases going northward up to 11.5 km (transect 1); going southward it remains greater than 6 km up to transect 37, decreasing further south.
The stretch ratio is a measure of the extensional strain of a linear element, which can be defined as the ratio between the final length and the initial length. In Post-LGM units (Figure 9G), it has been calculated on the basis of the maximum width of the area affected by fractures (as wide as 9.8 km); it ranges up to a maximum of 1.0035, with an average of 1.00084. Being generally constant in the northern part, it reaches a peak in the central part (at transect 21), and it gradually decreases in the southern part. Considering also Pre-LGM units, the maximum width of the area affected by fractures is 11.5 km, which results in a maximum stretch ratio of 1.107, with an average of 1.033. Two peaks are observed at transects 28 and 30, with a decrease in the southern part, whereas in the northernmost part an increase can be noticed.

Interaction Between Transform Faults and the ThFS Rift
The most relevant changes of fault and extension fracture strikes are observed in correspondence of the possible prolongation of the Husavik-Flatey Fault and in the northernmost part of the rift, near the sea coast in front of the Grimsey Lineament ( Figure 4A). Approaching the Husavik-Flatey Fault, the rift structures bend from a dominant NNE-SSW orientation to a NNW-SSE orientation. This can be better understood in Figure 10A, in which the faults and extension fractures located along the prolongation of the Husavik-Flatey Fault are plotted, with the exception of those striking NNE-SSW, which are more frequent north and south of this area. In light of the above, our observations corroborate the preliminary results provided by Tibaldi et al. (2018) who suggested that this anticlockwise reorientation of the rift structures is due to the possible role played by the right-lateral strike-slip Husavik-Flatey Fault underneath lava flows of Holocene age. As attested by analog modeling as well (Figure 10B), shallow right-lateral normal faults, striking NNW-SSE to N-S, can form at the surface with an en-échelon geometry, on top of a WNW-ESE transtensional right-lateral wrench fault affecting the basement.
The interaction between the structures of the ThFS and the Grimsey Lineament, in the northernmost part of the rift, can be interpreted in a different manner. Here, the rift structures tend to rotate clockwise, so much as to attain a preferential NE-SW orientation. The right-lateral, transtensional Grimsey Lineament is not located below the studied ThFS structures (as opposed to the above outlined setting, marked by the intersection of the structures with the Husavik-Flatey fault), but it is found further north. Hence, we can suggest a process, according to which the structures of the ThFS are dragged along the southern part of the shear zone, attaining a heterogeneous simple shear with a smoothly increasing strain gradient northwards ( Figure 11A). This setting, in fact, may result in a clockwise rotation of the structures as a consequence of the shearing deformation. In theory, an alternative interpretation may call for a rigid clockwise block rotation, that may also take place in correspondence of a right-lateral shear zone (Figure 11B; e.g., Nelson and Jones, 1987). However, in the latter case, the book-shelf rotation produces secondary zones of left-lateral shear in correspondence of the sides of the rotating blocks. On the contrary, our hundreds of field data point to a dominant right-lateral component of slip, whereas no fractures or faults with a dominant left-lateral strike-slip have been observed here. Although left-lateral strike-slip faulting has been suggested by Stefansson et al. (2008), to account for the 1885 earthquake, it is worth noting that there is no indication of fault direction in the micro-earthquake pattern of this seismic event, and only its location and magnitude have been obtained from intensity observations by Halldórsson (2005). Stefansson et al. (2008) suggested left-lateral slip along a NNE-SSW plane only on the basis of the distribution of small events associated with the large M6.2 1976 earthquake (Rögnvaldsson et al., 1998). Actually, this event took place 25 km further NE than the 1885 earthquake, in an offshore environment. The focal mechanism solution of the 1976 event might suggest left-lateral strikeslip along a NE-SW plane or right-lateral strike-slip along a NW-SE plane parallel to the Grimsey Lineament (Einarsson, 1987). Through mapping, mainly with the aid of aerial photos and satellite images, lineaments and fractures in the ThFS,  suggested that the local, right-stepping en-échelon arrangement of some of the NNE-SSW fractures might indicate left-lateral strike-slip displacements. However, our detailed survey along all these fractures does not confirm this kinematics. We thus can conclude that there are no data documenting book-shelf faulting in the onshore, northern portion of the ThFS.  Figure 8A).
Finally, also at other locations along the ThFS, a few segments of faults and fractures can be found, striking preferentially NNE to NE. In agreement with Bonali et al. (2019b), we also believe that these changes in strike orientation may derive from the structural inheritance of an older fault pattern located below the Holocene lavas. This hypothesis is based on the fact that within the conterminous, older rocks of Upper Miocene-Pleistocene age that crop out immediately west of post-LGM units (Figure 4A), two main sets of faults can be noticed, striking around N-S and NE-SW. The scarps of the NE-SW striking faults, in particular, stop abruptly at the stratigraphic contact with the post-LGM lava flows, indicating that the faults must prolong below the younger lavas. Such substrate faults may have locally guided the propagation of fractures upward with a similar NE-SW strike.

Opening Directions
Our extremely detailed survey of the opening directions along the entire set of structures has made it possible to unveil a quite peculiar pattern. From Figure 6D it can be seen that, as the fracture strike rotates, their opening direction changes as FIGURE 11 | (A) Sketch of dragging of the about N-S structures located south of the Grimsey Lineament right-lateral shear zone by heterogeneous simple shear with a smoothly increasing strain gradient. (B) Sketch of rigid clockwise block rotation (book-shelf type) that also may occur in correspondence of a right-lateral shear zone (after Nelson and Jones, 1987). Model B is not reconciled with field data collected in the northern part of the study area.
well. This may be claimed to reflect a process that is linked to fracture orientation, and only subordinately to the remote tectonic stresses. If the tectonic least principal stress (σ 3 ) would dominate, the opening direction should remain quite stable at fractures with different orientations. We thus suggest that the genesis of most extension fractures in the ThFS is linked to a strong local dilatant stress that forces the fracture walls to open along a direction from perpendicular to sub-perpendicular to the fracture strike. This type of stress can be produced by magma migrating in the crust along shallow-depth dykes. In fact, according to a commonly accepted model of dyke emplacement at shallow depth, the upward advancement of a dyke tip line should induce, at the surface, the development of a narrow graben or extension fractures above it (Pollard and Holzhausen, 1979;Mastin and Pollard, 1988;Rubin, 1992;Bonafede and Olivieri, 1995;Chadwick and Embley, 1998). In Iceland, such process has been observed, for example, during dyke emplacement events at Krafla (Opheim and Gudmundsson, 1989) and Holuhraun Ruch et al., 2016). In particular, Gudmundsson (1995), thanks to field observations in Iceland, suggested that dyke emplacement does not necessarily lead to graben formation but, instead, results more often in the development of extension fractures above it. We believe that, when magma propagates horizontally (i.e., laterally) along a shallow, vertical dyke, it can produce also fractures with strikeslip components as well. This is testified to by the onset of seismic events, contemporaneous to dyke intrusion, with double-couple focal mechanism solutions showing strike-slip components of deformation, as observed in the past at Long Valley caldera (Savage and Cockerham, 1984), during the 2014 event (Shelly et al., 2016) at Yellowstone caldera (Russo et al., 2017(Russo et al., , 2020, and during the 2014 Bárðarbunga dyking episode in Central Iceland (Ágústsdóttir et al., 2016). Figures 6A,B show that the southern and central sections of the rift display similar opening directions, in the N80-140 • range, whereas the northern part shows a systematic clockwise rotation of opening directions, that attain a N100-160 • range. This is reflected in both Northing and Easting variation because the rift is oriented NNE-SSW and, therefore, the measurement sites positioned further east correspond also to those located at the northern termination of the rift. We regard this clockwise rotation as induced by the same process described in the previous chapter, i.e., the effect of shear on the northernmost onshore part of the rift, produced by the transtensional right-lateral Grimsey Lineament.
The graphs of Figures 7A,C,D demonstrate the predominance of opening directions rotated in a clockwise direction relative to fracture strike, which are expressed as the dominant right-lateral component. From a quantitative point of view, 52% of extension fractures show a right-lateral component, 31% pure extension, and 17% a left-lateral component. Figure 7A is particularly insightful, as it shows that the presence of the right-lateral component cannot be accounted for in terms of a rotation of the strike of fractures. The dominant right-lateral component, instead, is systematically observed whatever the orientation of the fractures. The red stripes indicating how the lateral component is expected to change in case of variations in the fracture strike with respect to the plate spreading direction, separates a large amount of values where the right-lateral component is present. The reason for this opening pattern, marked by a widespread right-lateral component, can be found in the model of Figure 1A, where oblique rifting is normally found at rifts with structures whose strike is not perpendicular to the direction of plate separation (Fournier and Petit, 2007). A dextral shear along a fossil dyke swarm was found at Álftafjörður, in eastern Iceland, by Eriksson et al. (2015), and was interpreted as due to a paleo-spreading direction of N120-125 • , which is rotated slightly more northerly than present-day. In regard to the ThFS, another possibility is that its kinematics is also guided by the presence of the nearby Krafla rift and the interposed tectonic block. The latter has a wedge shape with the apex toward the south, where the two rifts come closer. A further complication may arise from different opening directions and/or amount along the Krafla rift. In order to look into these possibilities, we set out to map, with a very high detail, the Krafla structures as well.

Along-Rift Variation of Dilation
Along the ThFS, the cumulated dilation tends to decrease northward and southward from the rift's central area, as can be observed both considering deformation only in Post-LGM units and considering also Pre-LGM units ( Figure 9C). The variation in the number of fractures is more complex, but a general decrease outward from the rift's center is recognizable, though there are other peaks at transects 10-31 ( Figure 9D). This combination of the outward-directed decrease in fracture number and overall opening indicates that the fracturing process is more intense at the rift center near the Theistareykir volcano. In map view, the exact location of the volcano (situated between transects 25 and 27) corresponds to a lower amount of fractures, but this is due to the local presence of historic lavas that hide most structures. In Post-LGM units, the width of the deformation zone is also maximum in correspondence of the central part of the rift, while a minimum is reached in correspondence of the intersection with the HFF; going southward, values decrease ( Figure 9F). Considering also Pre-LGM units, the width of the area affected by dilation is more stable, with a decrease in the southernmost sector, but an important decrease can still be noticed in correspondence of the intersection with the HFF. The outward decrease of all these features can be explained in terms of a process of faulting and fracturing whose intensity varies in space, under the influence of tectonic and magmatic stresses. In regard to tectonics, plate separation here can contribute to triggering faulting processes, especially in the case of the studied normal faults that display vertical offsets up to 200 m. For major normal faults within rift zones, Acocella and Trippanera (2016) and Trippanera et al. (2015a,b) suggested that several dyke events can induce a cumulated amount of vertical displacement at the surface, in the order of several tens of meters. However, according to the present level of knowledge of the ThFS rift, we cannot rule out the possibility that such major fault displacements could be compatible with the long-term effects of multiple rifting events induced in part by amagmatic tectonic plate spreading (fault slip not related to dyke intrusion) and in part by shallow dyke intrusion events.
Deformation may be concentrated in correspondence of the central zone of the rift because, here, the heat flux from the magma chamber of the Theistareykir volcano tends to focus the tectonic stresses, thus promoting fault slip near the volcano, as well as producing the above cited, outward-tapering pattern. A reorientation of rifting in correspondence of shallow magma chambers was documented at oceanic ridges, such as, for example at the East Pacific Rise where Rea (1978) and Martinez et al. (1997) documented a shift in the spreading center axis location. As observed also at the Afar depression (Lahitte et al., 2003), active magma chambers can lead to the concentration of the maximum tensile stress, due to high temperatures that increase ductility and reduce the yield strength and the thickness of the brittle lithosphere.
Regarding extension fractures, several may have been produced by multiple shallow dyking events. As dykes can propagate horizontally from the magma chamber, magma pressure decreases with distance, as magma is drawn from the chamber. This, in turn, brings about a decrease in the number of dykes that can reach far distances. As a consequence of this, the intensity of dyke-triggered deformation decreases outward from the central part of the rift (Buck et al., 2006).
The width variation of the deformation zone is different in a northward direction; in Post-LGM units, first there is a drop from 9 to 1 km width, and then a slight increase from 1 to 3 km toward the sea coast ( Figure 9F). The same pattern can be observed if also Pre-LGM units are considered, with a drop from 9 to 3 km width, and then an increase from 3 to 11,5 km going northward. This geometry can be explained in two ways: on one side, it can depend on the fact that with distance from the magma chamber, dykes tend to become deeper and the width of the related surface deformation zone widens, as observed at the Krafla rift in Iceland by Hollingsworth et al. (2013). According to a second hypothesis, the northward-directed widening of the deformation zone may be the consequence of other dyke intrusion episodes here; in fact, the presence of a submarine volcanic system (Mánáreyjar, Figure 2) north of the ThFS (Hjartardóttir et al., 2009;Magnúsdóttir et al., 2015; may account for the southward dyke propagation from Mánáreyjar along the ThFS during the Holocene. The calculated stretch ratio in the study area is as much as 1.0035 considering only Post-LGM units, and 1.107 if we consider also the deformation in Pre-LGM units. The first value is very similar to the stretch ratio of 1.009 estimated along the Krafla Fissure Swarm, resulting from the deformation accumulated during the last 10 ka (Dauteuil et al., 2001). Over this period, both tectonic stresses and stresses produced by dyking events have contributed to the total cumulated surface deformation. The importance of multiple dyking events in producing stretching at the surface was documented also by Paquet et al. (2007), who calculated a stretch ratio of 1.036-1.046 across the Tertiary Alftafjordur dyke swarm, produced solely by intrusion events during about 1 Ma.

CONCLUSION
Our work documents, for the first time and with unprecedented detail, the overall deformation pattern along one of the most active rifts of Iceland, the Theistareykir Fissure Swarm (ThFS). We mapped the whole population of 2622 faults and 4187 extension fractures by field surveys at 1:5000 scale and Unmanned Aerial Vehicle surveys with cm resolution. We identified 685, N-S to NE-SW striking structures affecting pre-LGM units across the western part of the study area. We focused most of our attention on the 6124, post-LGM structures, so as to compare structures of the same age.
Most post-LGM faults and extension fractures strike N-S to NNE-SSW, but in the central part of the ThFS and at its northernmost part, they rotate in anticlockwise (striking NNW-SSE) and clockwise (striking NE-SW) fashion, respectively. The first case is explained in terms of the presence, underneath Holocene lava units, of the Husavik-Flatey right-lateral strike-slip fault, which brings about fault rotation at the shallowest crustal level. The NE-SW strike of the structures in the northern ThFS is interpreted as due to the dragging effect along the block located south of the transtensional, right-lateral Grimsey Lineament, and to the structural inheritance of older faults in the substrate with the same NE-SW orientation that guided successive fracture propagation up to the surface, through post-LGM lava flows.
Along the ThFS, the cumulated dilation, expressed as the sum of the opening measured at all fractures, decreases southward and northward with respect to the central part of the rift, where the Theistareykir central volcano is located. Outward from this volcano, there is also an along-rift decrease in the number of structures, and the deformation zone becomes narrower. The above evidence can be explained by taking into consideration the outward propagation of dykes from the Theistareykir magmatic chamber; with distance, magma pressure decreases as magma is drawn away from the chamber, and this results in a decrease in the number of dykes that can reach far distances. As a consequence of this, the intensity of dyke-triggered deformation decreases outward from the central part of the rift.
The increase in the number of faults and in their dip-slip offsets at the central part of the rift, with special reference to the large faults located west of the Theistareykir volcano, is related to the heat flux connected with its magma chamber; the higher temperature increases ductility and reduces the yield strength and the thickness of the brittle lithosphere, focusing tectonic stresses and favoring fault slip. The proximity to the magma chamber is also associated with greater magma pressure, which may contribute to increasing dyke-induced fault slip especially along the rift axis, which is to say immediately north and south of the Theistareykir volcano.
Most extension fractures show components of right-lateral and left-lateral shear, with deviations as much as 40 • from pure extension. This kinematics may be explained with the earlier proposed models related to different processes: (i) repeated horizontal propagation of vertical dykes that produces strike-slip components of motions, accompanied by magma inflation and separation of the fracture walls; (ii) dragging along the Husavik-Flatey Fault and Grimsey Lineament; (iii) tectonic events related to regional stresses and plate motions. Cumulated episodes may lead to a resulting very complex kinematics in a volcanotectonic rift zone, where local processes interact with far-field tectonic stress.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
AT, FB, and FP wrote the first draft of the manuscript. NC, ER, PE, and ÁH wrote sections of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

FUNDING
This work was a contribution to the International Lithosphere Program -Task Force II. Funding for fieldwork were from ILP and from the MIUR project ACPR15T4_00098 -Argo3D (http: //argo3d.unimib.it/). Furthermore, this study has been carried out in the framework of the Project MIUR -Dipartimenti di Eccellenza 2018-2022, and European Space Agency project no. 38829. Finally, this effort was also an outcome of GeoVires, the Virtual Reality Lab for Earth Sciences hosted at the Department of Earth and Environmental Sciences, University of Milan Bicocca (Italy). Agisoft Metashape was also acknowledged for photogrammetric data processing.