Interseismic deformation of the Húsavík–Flatey fault (North Iceland) from two decades of GNSS data

SUMMARY

The Tjörnes Fracture Zone in North Iceland is one of two transform zones in Iceland capable of generating earthquakes of magnitude |$\sim$|7. More than 150 years have passed since the last two major earthquakes occurred on the Húsavík–Flatey fault, one of the two main transform structures within the Tjörnes Fracture Zone. Given the seismic hazard posed to Húsavík and adjacent coastal communities, accurately determining the slip rate and locking depth of this fault is crucial for a robust earthquake hazard and risk assessment. In this study, we significantly expand the existing global navigation satellite system data set for the Tjörnes Fracture Zone by incorporating more than a decade of additional data and doubling the number of stations. This expansion not only improves the spatial coverage of the network, but also refines the station velocities. We present an updated interseismic velocity field for North Iceland and implement a backslip model with nine dislocation segments to describe the plate boundary deformation. Additionally, we include a point pressure source for the ongoing broad uplift signal in the study area. Our findings indicate a locking depth of |$7.3\substack{+0.9 \\ -0.7}$| km and an average slip rate of |$6.9\pm 0.2$| mm yr⁻¹ for the Húsavík–Flatey fault. With our updated approach, we can narrow down model parameter constraints from previous studies and thereby provide an enhanced understanding of the earthquake potential of this region.

Satellite geodesy, Numerical approximations and analysis, Time-series analysis, Kinematics of crustal and mantle deformation, Oceanic transform and fracture zone processes

1 INTRODUCTION

Located on the northern Mid-Atlantic Ridge and atop a mantle plume, Iceland provides a unique setting where plate extension between Eurasia and North America can be observed above sea level (Sæmundsson 1974; Einarsson 2008; Karson 2017). The westwards migration of the Mid-Atlantic Ridge relative to the plume has caused eastwards jumps of the subaerially exposed plate boundary, with new rifts developing in zones of weakness created by the hotspot (Einarsson 2008). Consequently, the active rifts maintain their position approximately over the plume (Burke et al. 1973; Einarsson 2008; Karson 2017). The offsets generated by these eastward jumps are accommodated in transform zones in northern and southern Iceland (see inset in Fig. 1).

Figure 1.

Tectonic setting, seismicity and GNSS stations in North Iceland. The inset shows how the mid-Atlantic plate boundary is offset to the east by the South Iceland Seismic Zone (SISZ) and the Tjörnes Fracture Zone (TFZ). The bathymetry is from the GEBCO_2023 Grid (GEBCO Bathymetric Compilation Group 2023). Locations of M > 6 earthquakes (yellow stars, Stefánsson et al. 2008) in the past 300 years and instrumental seismicity 2009–2022 (magenta dots) are shown, as well as campaign (purple triangles) and continuous (pink triangles) GNSS stations. Fissure swarm fractures and faults (brown lines), central volcanoes (orange outlines) and onshore fissure swarms (yellow areas) are also shown. KR: Kolbeinsey Ridge, ER: Eyjarfjarðaráll Rift, GOR: Grímsey Oblique Rift, HFF: Húsavík–Flatey Fault, NVZ: Northern Volcanic Zone, EVZ: Eastern Volcanic Zone, RR: Reykjanes Ridge, Th: Þeistareykir, Kr: Krafla, Fr: Fremrinámar, DL: Dalvík Lineament, NA: North-American plate, TB: Tjörnes block, EU: Eurasian plate.

To the north, the Tjörnes Fracture Zone (TFZ) is an |$\sim \!$|120 km oceanic transform offset that connects the offshore Kolbeinsey Ridge to the onshore Northern Volcanic Zone (NVZ; Einarsson 1991; Stefánsson et al. 2008; Homberg et al. 2010). The relative plate motion at this latitude is 18 mm yr⁻¹ with an azimuth of N104°E, according to the MORVEL plate-motion model (DeMets et al. 2010), and is accommodated in the TFZ by two subparallel structures: The Húsavík–Flatey Fault (HFF) and the Grímsey Oblique Rift (GOR) (see Fig. 1).

The GOR, located entirely offshore, is characterized by oblique tectonics (Einarsson 1976, 1991, 2008) and accommodates both normal and strike-slip displacements on steeply dipping faults arranged en échelon and striking roughly N–S. Graben-like troughs show the same trend in bathymetric, magnetic and seismic data (McMaster et al. 1977). The GOR exhibits Holocene volcanism distributed across four volcanic systems (Brandsdóttir et al. 2005; Magnúsdóttir et al. 2015), with at least two eruptions in historical times, northwest of Grímsey in 1372 and within the Mánáreyjar volcanic system in 1867–1868 (Thorarinsson 1965; Magnúsdóttir & Brandsdóttir 2011; Magnúsdóttir et al. 2015). Present-day seismicity is mostly associated with N–S oriented, left-lateral strike-slip faults (Rögnvaldsson et al. 1998). Large earthquakes in the GOR include a magnitude M6.3 event in 1885 towards its southeastern part and an M7 earthquake in 1910 in its central region (Stefánsson et al. 2008). In 1976, a M6.2 event in Öxafjörður was triggered during the initial stage of the Krafla rifting episode that started in 1975 (Einarsson 1979; Tryggvason 1980, 1984; Björnsson 1985; Passarelli et al. 2013).

The HFF is a right-lateral strike-slip fault system, roughly extending 100 km from the southern continuation of the Kolbeinsey Ridge, the Eyjafjarðaráll Rift, to the NVZ (Fig. 1). Unlike the GOR, the HFF has no associated Holocene volcanism (Einarsson 2008). The HFF can be well traced in bathymetric data (Brandsdóttir et al. 2005), and while mostly offshore, the HFF enters land just north of the town of Húsavík, extending for about 25 km into the Þeistareykir volcanic system (Sæmundsson 1974; Gudmundsson et al. 1993; Garcia & Dhont 2005; Pasquarè Mariotto et al. 2015). The orientation of the fault changes slightly between its offshore and onshore segments, becoming more oblique-transtensional onshore. The presence of four pull-apart basins along the length of the fault (Brandsdóttir et al. 2005) further indicates that the HFF accommodates transtensional deformation (Garcia & Dhont 2005).

The HFF system is well defined by seismicity in the west (Einarsson 1991, 2008). Passarelli et al. (2018) reported six earthquake swarms between 1997 to 2015, mostly confined towards the western portion of the HFF and the Eyjafjarðaráll Rift, which were possibly associated with transient aseismic slip. However, no evidence for aseismic strain release could be deduced from global navigation satellite system (GNSS) observations. Another energetic earthquake swarm took place in 2020 with moment magnitude (Mw) Mw5.4 and Mw5.7 earthquakes near the western end of the HFF, followed by a Mw6.0 normal-faulting earthquake on the eastern side of the Eyjafjarðaráll Rift (Viltres et al. 2022). Historical earthquakes on the HFF include one M7 event in Skjálfandi in 1755, one M6.5 event in 1838 towards its western end, and two M6.5 events in 1872 close to Húsavík and Flatey (Stefánsson et al. 2008). After the Krafla rifting episode, the seismicity on the HFF has been mostly limited to its northwestern half (Tryggvason 1980; Einarsson 1991; Maccaferri et al. 2013). Maccaferri et al. (2013) concluded that the decrease in seismicity of the eastern half of the fault was due to a stress shadow cast by the Krafla rifting episode that took place between 1975–1984.

The existence of a third seismic zone within the TFZ, the Dalvík Lineament (DL), see Fig. 1, or Dalvík Zone (between DL and HFF), has been a matter of debate. Although, there is little structural evidence of a strike-slip fault parallel to the HFF in this area (Einarsson 1991; Rögnvaldsson et al. 1998; Garcia & Dhont 2005), the zone is marked by diffuse seismicity 30 km south of the HFF (Einarsson 1991, 2008; Fig. 1). Two major earthquakes are known in this zone, one of M6.2 in 1934 close to Dalvík, and the M7 Skagafjörður earthquake in 1963 (Einarsson 1976, 1991, 2008; Stefánsson et al. 2008, Fig. 1). Rögnvaldsson et al. (1998) proposed that faulting in the Dalvík Zone might bear some similarities to that of the GOR, with en échelon faults with normal and left-lateral components striking N–S. This is supported by N–S lineations of earthquakes (Stefánsson et al. 2008). However, geodetic studies have not yet reported strain accumulation across the Dalvík zone (Metzger et al. 2011, 2013; Metzger & Jónsson 2014).

The NVZ is a 200-km-long and 50-km-wide extensional plate boundary zone, with seven volcanic systems, each consisting of a central volcano and an associated fissure swarm (Hjartardóttir et al. 2016) (Fig. 1). The volcanic systems are arranged en échelon due to the obliquity of the plate spreading with respect to their trend (Einarsson & Brandsdóttir 2021) and extend to the north coast, linking to the TFZ. The four northernmost volcanic systems are Þeistareykir, Krafla, Fremrinámar and Askja. The main rift axis appears to go through the central volcanoes of Krafla, Fremrinámar and Askja (Drouin et al. 2017a). The last three major rifting episodes occurred in Askja (1874–1875), Krafla (1975–1984) (Tryggvason 1980, 1984; Björnsson 1985), and farther south in Bárðarbunga-Holuhraun (2014) (Sigmundsson et al. 2015; Ruch et al. 2016). The Krafla episode corresponded to 275 yr of plate spreading and was characterized by more than 20 intrusive events and crustal widening of over 5 m (Tryggvason 1994). It also caused a transient deformation pulse that propagated and decayed in amplitude away from the rift axis (Foulger et al. 1992; Heki et al. 1993; Hofton & Foulger 1996). By 2000, or at the latest by 2005, the deformation rates returned to pre-rifting rates (Völksen 2000; Árnadóttir et al. 2009).

The 2014 Bárðarbunga-Holuhraun dyke intrusion (Sigmundsson et al. 2015; Ruch et al. 2016) caused detectable deformation up to 80 km away from the dyke (Ofeigsson et al. 2015). Jónsson et al. (2019) pointed out displacements even farther away, to distances of over 200 km, with transient motions of about 6 mm to the northwest, indicating that the deformation affected all of North Iceland. Additional deformation signals within the TFZ include: (1) a broad inflation signal north of the Krafla caldera since at least 1993 (de Zeeuw-van Dalfsen et al. 2004; Metzger & Jónsson 2014), (2) localized deflation signals of likely geothermal origin southwest of the Krafla caldera (Drouin et al. 2017a), (3) inflation at Þeistareykir between 2007 and 2008, and (4) deflation at Þeistareykir following the previous inflation period. The two latter deformation events were modelled as point pressure sources with an estimated depth of |$\sim$|8.5 km (Metzger et al. 2011, 2013; Metzger & Jónsson 2014; Drouin et al. 2017a).

In addition to tectonic and volcanic deformation, glacial isostatic adjustment (GIA) plays a dominant role in the large-scale deformation over Iceland, particularly in the vertical component, and also generates horizontal displacements away from the diminishing load (Árnadóttir et al. 2009; Drouin et al. 2017a; Cao et al. 2023). The ongoing retreat of the Icelandic ice caps since at least 1890 (Adalgeirsdóttir et al. 2020) generates uplift rates of up to 30 mm yr⁻¹ and uplift accelerations 1–2 mm yr⁻² in the Central Highlands, as evidenced by GNSS measurements between 1995 and 2014 (Compton et al. 2015). Uplift rates in the southern part of the NVZ reach 20 mm yr⁻¹, and the effects of the GIA extend to the northern coast (Drouin et al. 2017a). Drouin et al. (2017a) and Cao et al. (2023) modelled the GIA in GNSS and Interferometric Synthetic Aperture Radar (InSAR) data by scaling the GIA model (Auriac 2014) to the geodetic observations. Their results require different scaling factors for the horizontal and vertical components.

Previous studies of northern Iceland provided valuable insights into the slip rate and locking depth of the HFF. Based on seismicity data collected between 1994 and 1997, Rögnvaldsson et al. (1998) inferred that the base of the brittle crust, which can be considered a proxy for the locking depth, is at 10 km depth. They also observed an increase in the depth of hypocentres from 5 to 10 km with increasing distance from the Eyjafjarðaráll Rift in the western part of the HFF. Earthquake depths of a seismic swarm that took place in 2013 on the western HFF reached 8.0–11.5 km depth (Jónsdóttir et al. 2016). Additionally, local earthquake tomography results in the TFZ revealed low velocities associated with the HFF and GOR to approximately 10 km depth (Abril et al. 2021). The same study revealed velocity anomalies at 5–10 km deep that were interpreted as supercritical fluids in the fractures (Abril et al. 2021) beneath most of the HFF.

The HFF has probably existed since the NVZ became active 8.0–8.5 Ma (Garcia et al. 2003). This would indicate 150 km transfer motion if the HFF had been the only active transform fault in North Iceland. However, there were two parallel rifts in North Iceland until |$\sim$|3 Ma (Garcia et al. 2003) and since then, the GOR has accommodated part of the transfer motion. Thus, the total offset of the HFF is smaller and has been estimated to be up to |$\sim$|60 km (Sæmundsson 1974). Deducing slip rates from this offset is difficult given the complex rifting history. Based on a single offset measurement of |$\sim$|12 ka-old lava in the vicinity of the Þeistareykir fissure swarm, Rust & Whitworth (2019) proposed a Holocene slip rate of 3.8 mm yr⁻¹ for the HFF. More recently, (Matrau 2023) estimated a slip rate of 6.0 to 6.5 mm yr⁻¹ based on the Pleistocene and postglacial fault offsets.

Árnadóttir et al. (2009) used nationwide episodic GNSS (eGNSS) data from 1993 and 2004 to estimate a locking depth of |$\sim$|5 km for the HFF and a slip rate of less than 5 mm yr⁻¹. Geirsson et al. (2006) used the first three continuous GNSS (cGNSS) sites in North Iceland and found that the HFF accommodated about 40 per cent of the motion across the TFZ. Jouanne et al. (1999, 2006) analysed velocities from eGNSS sites between 1995 and 2002 and constrained the slip rate to 8 mm yr⁻¹. Their profile was not long enough to estimate the locking depth of the fault; however, they suggested it was slightly larger than the 10–12 km that could be inferred from microseismicity (Jouanne et al. 2006). More recent geodetic studies (Metzger et al. 2011, 2013; Metzger & Jónsson 2014) combined data for 14 cGNSS stations from 2006–2010 and 44 eGNSS sites between 1997 and 2011 with InSAR data from 1992–2009, and deduced a locking depth of 6–10 km and a slip rate of 6–9 mm yr⁻¹ for the HFF.

Previous estimates of the slip rate and locking depth of the HFF were derived from geodetic data with a scattered distribution and limited temporal coverage. In this study, we utilize two decades of GNSS measurements that offer enhanced spatial coverage to model plate boundary deformation in the TFZ. Our goal is to improve the constraints on the slip rate and locking depth of the HFF; parameters that are crucial to understanding the seismic hazard in the region.

2 GNSS DATA PROCESSING

A continuously operating GNSS station network was first established in Iceland in 1995, with station REYK in Reykjavík, followed by station HOFN on the southeastern coast of the country in 1997 (Geirsson et al. 2006). The first cGNSS stations in North Iceland (AKUR, RHOF and ARHO; Fig. 1) were installed between 2001 and 2002. In North Iceland, eGNSS measurements began in 1986, focusing on deformation following the Krafla rifting episode (Foulger et al. 1992; Heki et al. 1993; Hofton & Foulger 1996). An eGNSS network to study the interseismic deformation in the TFZ was installed in 1995 (Jouanne et al. 1999) and, since then, the network has been steadily expanded, reaching 24 continuous and 92 episodic sites by 2024 (Barreto et al. 2022). In addition to a higher station density, the current network has a better spatial coverage, particularly towards the west covering the Dalvík Zone.

We processed data from a total of 23 cGNSS stations in North Iceland, spanning from 2001 to 2022, and 108 eGNSS sites, including data from 16 sites of the Icelandic nationwide GNSS network (ISNET) acquired in 1993, 2004 and 2016. All of the 92 eGNSS sites in North Iceland have been measured at least twice, but most have been measured at least five times. Along with the data collected in 1997, 1999, 2002, 2007 (Jouanne et al. 2006) and 2009, 2010, 2011 (Metzger et al. 2013), we processed new episodic data from 2013, 2016, 2019 and 2022. Both the continuous and episodic data were sampled every 15 s. We processed the Receiver Independent Exchange files using the GAMIT/GLOBK 10.71 software suite (Herring et al. 2018) following the steps described in Floyd et al. (2010), and Kogan et al. (2012). The continuous and episodic data were processed separately, using a shared set of 19 International GNSS Service (IGS) stations (ALGO, CHUR, GODE, GRAS, MAR6, NRC1, NYA1, ONSA, REYK, HOFN, KIRU, YEBE, PDEL, SCH2, STJO, THU2, TRO1, WTZR, ZIMM), and the data were combined later.

First, we produced loosely constrained daily solutions of the station coordinates and their covariance matrices from double-difference GNSS phase observations. Second, we used the GLOBK Kalman filter (Herring et al. 2018) to align the daily solutions into the International Terrestrial Reference Frame, ITRF14 (Altamimi et al. 2017), and generated daily position time-series. For the continuous data, outliers were removed iteratively and offsets were accounted for. Then, we calculated a six-parameter Helmert transformation that minimized the difference between the ITRF14 (Altamimi et al. 2017) and the 19 IGS stations, defining a stable reference frame to produce more accurate daily position time-series.

Episodic GNSS time-series were also inspected for outliers, and downweighted appropriately to reduce unwanted effects in the velocity estimation. We used GLOBK to produce a single position solution per episodic survey and to combine episodic and continuous data. We used the software HECTOR 1.9 (Bos et al. 2013) to account for the temporally correlated noise, and estimated the coloured-noise component of the cGNSS time-series as a combination of Flicker and white noise. Then, we added the median coloured noise of the cGNSS stations in the region to the episodic stations as equivalent random-walk magnitude during the Kalman-filter stage in GLOBK (Wedmore et al. 2021). This allowed us to obtain more realistic uncertainties of the estimated station velocities.

3 TRAJECTORY MODELLING AND VELOCITY ESTIMATION

GNSS signals are usually dominated by a secular trend of tectonic origin and seasonal oscillations due to hydrological loading. In tectonically active regions, factors such as volcanic unrest, earthquake coseismic displacements, fault-creep events and post-seismic relaxation can also impact the time-series (Freymueller 2021). The trajectory model for the eGNSS time-series only contained a secular trend (see Figs S1 and S2 in the Supporting Information). Seasonal effects were suppressed by carrying out the surveys during the same time of the year (July to early September).

For the cGNSS, we performed the time-series decomposition assuming the observations are a sum of a deterministic model plus stochastic noise (Flicker + white noise). Following the work of Bevis & Brown (2014), our cGNSS trajectory model included a secular trend, annual and semi-annual oscillations, and instantaneous jumps caused by earthquakes or antenna changes. The final trajectory model can be expressed as:

$$\begin{eqnarray} x(t) & =& x_0 + v\cdot (t - t_0) \\ && + \sum _{i=1}^{n} \left[ A_i \sin (\omega _i t) + B_i \cos (\omega _i t) \right] \\ && + \sum _{j=1}^{m} D_j H(t - t_j) , \end{eqnarray}$$

(1)

where |$x(t)$| is the position at time t, |$x_0$| is the initial position and v is the velocity. The seasonal terms are modelled by a sum of sines and cosines, with n corresponding to the number of frequencies used to model the periodic terms (i.e. |$n = 2$| for annual plus semi-annual oscillations). |$A_i$| and |$B_i$| correspond to the amplitudes of the oscillations and |$\omega _i = 2\pi /\tau _i$|⁠. The period |$\tau _i$| is defined for the annual and semi-annual terms. m is the number of offsets at |$t_j$| of magnitude |$D_j$| modelled by a Heaviside function |$H(t - t_j)$|⁠. The equation is applied to each velocity component (east, north, up) of each time-series independently. We used the software package HECTOR 1.9 (Bos et al. 2013) to estimate the time-series components and the temporally correlated noise.

3.1 Modelled time-series

The trajectory modelling captures the dominant linear trend and seasonal oscillations in the observed time-series (Fig. 2 and Fig. S3 in the Supporting Information). The results also show that some residual seasonal signals remain, for example, at stations AKUR and KOSK. This suggests year-to-year or multiyear differences in the hydrological cycle might be better modelled by incorporating a time-varying seasonal term or adding more sine and cosine terms to the original trajectory model. However, the residuals of most stations do not exhibit seasonality (Fig. S4 in the Supporting Information). Additional signals, such as transients of volcanic origin [e.g. a brief phase of volcanic unrest at Þeistareykir between 2007 and 2008 affecting station GRAN (Metzger et al. 2011)] are considered to be too short—relative to the total length of the time-series—to affect the velocity estimates. However, more prominent transients of likely volcanic origin affect station MYVA (Figs S3 and S4 in the Supporting Information).

Figure 2.

Time-series and residuals for six continuous GNSS stations in North Iceland: AKUR, GRAN, KVIS, ARHO, KOSK and RHOF (see station locations in Fig. 1) The east, north and up components are presented relative to stable North America as defined by (Altamimi et al. 2017). Daily positions with their respective uncertainties (light grey) and residuals (dark grey) are displayed in mm. The best-fitting trajectory model prediction (magenta) accounts for a linear trend, seasonal terms and offsets (eq. 1). Position outliers were removed during pre-processing.

The residual time-series for the east and north components plot tightly around zero, while the up component shows significant deviations. The vertical residuals display a clear concave-down shape, particularly evident for stations AKUR, GRAN, ARHO and RHOF (Figs S3 and S4 in the Supporting Information), which can be attributed to the effects of accelerating GIA (Compton et al. 2015; Cao et al. 2023). Accelerating GIA might also explain the slight concave residuals observed in the horizontal components of some stations (AKUR, ARHO, RHOF, SJUK; Fig. S4 in the Supporting Information).

To account for the concave residuals, we tested a trajectory model including a second-degree polynomial term in the trajectory model: |$\frac{1}{2} \cdot a(t-t_0)^2$|⁠, where the quadratic term corresponds to half the acceleration |$(a)$| in the time-series (Bos et al. 2013; Bevis & Brown 2014). While acceleration parameters in the east and north components are almost negligible, most stations show negative accelerations between 0.1 and 0.4 mm yr⁻² in the vertical component (Figs S5–S7 in the Supporting Information). Outlier accelerations are observed at some stations with time-series that are too short to get reliable estimates (i.e. less than five years). For instance, stations such as BRIK, BRTT and HRAC exhibit notably large negative accelerations. Conversely, station MANA, which also has a relatively short time span, shows a positive acceleration. Short time-series are also associated with outlier accelerations in the other components (Fig. S8 in the Supporting Information). Given their consistency, we suggest the observed accelerations might be due to a common process that affects all stations in North Iceland, possibly related to the accelerating GIA displacements reported by Compton et al. (2015). However, the acceleration values need to be addressed carefully, as they contain large uncertainties, and may not be entirely related to a single physical process.

Including a second-degree polynomial in the trajectory model improves the fit to the observations and removes the concavity from the residuals (Fig. S6 in the Supporting Information), but does not significantly alter the estimated secular trends of the time-series. The differences in velocity estimates between the models with and without acceleration are less than 1 mm yr⁻¹ for all components of all stations. The mean absolute difference between the both sets of velocities is 0.11 mm yr⁻¹ for the east component, 0.08 mm yr⁻¹ for the north component, and 0.32 mm yr⁻¹ for the up component, which is well below the average of the velocity uncertainties. Since we are primarily interested in the secular velocities, we decided to use the original trajectory model presented in eq. (1).

3.2 Updated velocity field

A gradual spatial change in velocities can be observed in Fig. 2 and Fig. S3 in the Supporting Information, in which the cGNSS stations are arranged from top to bottom from the North-American to the Eurasian plate. Using a stable North-American reference frame (Altamimi et al. 2017), AKUR and GRAN on the North-American plate still show slight movement to the west and to the north of approximately 2 mm yr⁻¹ to the west and 3-4 mm yr⁻¹ to the north (Fig. 2). This systematic northwestward plate velocity has been reported in previous studies (Árnadóttir et al. 2009; Geirsson et al. 2010; Metzger et al. 2011, 2013; Michalczewska et al. 2013) and is attributed to either a small reference frame problem in North Iceland or the study region not being completely stable. Therefore, all velocities are corrected for the small reference frame shift (see Section 4).

The addition of over ten years of data since the last results were published (Metzger & Jónsson 2014) in the area yields lower uncertainties in our estimated velocities, particularly for the vertical component. For the common stations in the previous and current studies, the mean uncertainties decreased from 0.3 to 0.2 mm yr⁻¹ for the east and north components and from 1.6 to 0.6 mm yr⁻¹ for the up component. The mean uncertainties for our updated velocity field and increased station network are 0.30, 0.41, and 0.94 mm yr⁻¹ for the east, north and up components.

The overall velocity field exhibits consistently increasing rates in the ESE direction, from velocities near zero on the North-American plate to rates approaching the full expected plate motion towards the ESE on the Eurasian plate (Fig. 3a), such as station RHOF which has a velocity of 18.6 mm yr⁻¹ in an ESE direction. Velocities at the northern tip of the Tjörnes peninsula display roughly half the total rate observed for the sites farthest away from the plate boundary on the Eurasian plate, for instance station ARHO located on the Tjörnes peninsula has an ESE velocity of 8.4 mm yr⁻¹.

Figure 3.

(a) Horizontal GNSS velocities and 95 per cent confidence ellipses relative to stable North America (Altamimi et al. 2017), corrected for a small reference frame shift that affects all stations in northwestern Iceland (Section 4). (b) Vertical GNSS velocities. The larger circles correspond to cGNSS and the small circles correspond to eGNSS. The dashed brown lines indicate the numbered rectangular dislocation segments representing the plate boundary (Section 4.1). Segment 1 represents the Kolbeinsey Ridge (KR), which extends north of the map. Stations within the dashed purple rectangle in (a) are used to plot the fault-parallel velocity profile of Fig. 7. The green star marks the location of the point pressure source used in the kinematic model. Other features are the same as in Fig. 1. The inset marked by the dashed green square shows more details of the stations located at the Tjörnes peninsula.

Several stations near the Krafla volcanic system (GRAE, MYVA, HVIT, KRAF) exhibit deviations from the general trend of plate divergence, reflecting local deformation (Fig. 3a). Other sites with anomalous velocities include RAND within the Þeistareykir volcanic system, and ISHO, the southernmost station in our study area. Station RAND shows a larger southward component than its neighbours, and station ISHO presents an increased NNE velocity. Additionally, two sites on the Tjörnes block close to the HFF (REYB, NYKU) show increased southward velocities, though these values are likely biased as these sites have a limited number of observations.

Most vertical rates range between −5 and 6 mm yr⁻¹, with uncertainties between 0.30 to 1.78 mm yr⁻¹ for cGNSS stations and from 0.46 to 4.85 mm yr⁻¹ for eGNSS sites. The vertical velocity field shows a transition from uplift at inland stations to subsidence near the coast, especially towards the northwest of the study area (Fig. 3b). Although the vertical velocity field appears consistent and aligns with InSAR results from Cao et al. (2023), vertical rates need to be analysed carefully, particularly those of campaign stations.

The southernmost stations, ADAL and ISHO, exhibit significant uplift, with ISHO displaying the largest vertical velocity in the study area of around 12 mm yr⁻¹. Uplift rates higher than 5 mm yr⁻¹ can also be found close to the NVZ towards the south at stations GRAE and KRME. However, stations within the Krafla volcanic system, such as MYVA, HVIT and KRAF, show vertical velocities close to zero, with HVIT experiencing subsidence of about 2 mm yr⁻¹. Two clear outliers can be identified in Fig. 3(b): stations TH03 and TR10, which show anomalously high (one order of magnitude) subsidence and uplift, respectively. However, only two measurements (2019 and 2022) are available for these stations. Stations ADAL, ISHO, GRAE and KRME only have three measurements available (1993, 2004 and 2016).

We did not correct our velocities for GIA. The scaling method of the GIA model from Auriac (2014) does not predict the subsidence observed in the northwest of our study area, as evidenced by the negative residuals in North Iceland in the work of Drouin et al. (2017a) and Cao et al. (2023). Moreover, Drouin et al. (2017a) and Drouin et al. (2017b) showed that the modelled GIA contributions in the Krafla area are only a few mm yr⁻¹ for the vertical component and less than one mm yr⁻¹ for the horizontal components. These GIA corrections are even smaller farther north. Therefore, as our primary interest are the tectonic movements causing horizontal deformation, and GIA corrections mainly affect the vertical component, we proceed to model the plate boundary deformation without including a GIA correction. However, we still retain the vertical component in our modelling.

4 BACKSLIP MODELLING

4.1 Model description

Kinematic modelling of the observed velocity field allows for the estimation of key fault parameters, such as the locking depths and slip rates of the main structures accommodating plate divergence. The interseismic deformation produced by the relative motion of two blocks separated by a vertical strike-slip fault, which is locked to a certain depth but moves freely below (Savage & Prescott 1978), was modelled by superimposing steady backslip creep on rigid block motion. The backslip is of equal magnitude but in the opposite direction to that of the rigid block motion, and is confined to the locked portion of the fault. This results in velocities close to zero near the locked interface, while the far-field shows full plate velocities. Additionally, we allowed the model to include an opening component to simulate the transtensional character of the TFZ (Metzger et al. 2011).

The TFZ and the rifts segments are described by nine rectangular dislocation segments in an elastic half-space (Okada 1985). Each segment is fully characterized by ten parameters: seven describing the geometry and three describing the displacement on the segment. Of the seven geometrical parameters (dislocation length, width, depth, strike, dip angle, east and north location), we fix six and only estimate the locking depth. The dislocation segments are illustrated in Fig. 3(b), and the fixed parameters are described in Table 1 following Metzger & Jónsson (2014). Of the three kinematic parameters (dip-slip, strike-slip and opening), we do not allow for dip-slip and derive the other two parameters from the relative plate motion and the strike of each segment. In the TFZ, a portion of the total plate motion between the Eurasian and North-American plates is accommodated by motion of the Tjörnes block (TB in Fig. 3). Consequently, the combined slip rates of the segments east (2–4) and west (5–8) of the Tjörnes block sum up to the full plate motion (Fig. 3b). Segments 1 and 9 remain unaffected by the motion of the Tjörnes block and accommodate the full plate extension.

Table 1.

Location, length and orientation of the vertical model dislocation segments, as well as the location and depth of a point pressure source following Metzger & Jónsson (2014). The segment coordinates are given in Universal Transverse Mercator, zone 28W.

Segment	Segment	Easting	Northing	Length	Strike
#	Type	[km]	[km]	[km]	[N|$^\circ$| E]
1. KR	|$ridge$|	420.79	7774.04	672.2	13.6
2. GOR	|$ridge$|	352.25	7421.66	55.7	157.4
3. GOR	|$transform$|	395.95	7371.85	81.7	126.1
4. NVZ	|$ridge$|	425.19	7327.91	40.4	10.8
5. KR	|$ridge$|	334.41	7419.87	56.8	14.6
6. ER	|$ridge$|	328.44	7375.54	33.8	175.9
7. HFF	|$transform$|	361.62	7342.93	71.3	116.3
8. HFF	|$transform$|	407.50	7317.60	33.7	124.4
9. NVZ	|$ridge$|	347.03	6984.44	664.1	12.9
	Mogi	Easting	Northing	Depth
	Source	[km]	[km]	[km]
	Krafla North	420.963	7302.235	21.5

With this model setup, only five parameters need to be estimated to describe the plate boundary kinematics: the azimuth and magnitude of the total plate motion, the proportion of the plate motion accommodated by the HFF (with the rest by the GOR), and two locking depths for the transform-type and ridge-type segments in Table  1. To account for the reference frame shift that affects the velocities at all stations (Metzger et al. 2011, 2013; Michalczewska et al. 2013; Metzger & Jónsson 2014) two correction parameters (shifts in the east and north directions) were also estimated. For the broad uplift signal that influences the stations north of the Krafla central volcano, we estimated the volume change of a point pressure source (Mogi 1958) at a fixed depth of 21.5 km (Metzger & Jónsson 2014). Thus, the final model (henceforth called Model 1) contains eight free model parameters.

We excluded the following stations with anomalous velocities from our kinematic modelling: GRAE, MYVA, HVIT and KRAF near the Krafla volcano due to significant local deformation. We also excluded stations REYB and NYKU on the Tjörnes block, as they exhibit increased southward velocities compared to their neighbours, and stations TH03, TR10 and ISHO, which show anomalous vertical rates (Fig. 5a). The vertical component of station ADAL was downweighted, but retained in our analysis. The excluded stations displayed clear inconsistencies with neighbouring stations and large residuals at the early stages of the modelling in at least two of their three components, degrading the fit to the model prediction considerably. The decision to retain ADAL was based on it showing acceptable horizontal residuals.

Figure 4.

(a) Histograms and scatter plots representing 1-D and 2-D marginal probability density functions and the estimated model parameters. The scatterplots include the best-fitting model realizations (orange square) from the 1000 randomize-then-optimize model runs (grey points). ldRID: locking depth of the ridge segments, ldHFF: locking depth of the HFF and all other transform segments, vTot: total plate motion; prcnt: percentage of the total plate motion taken up by the HFF, vAzi: azimuth of the total plate motion, mogV: volume change rate of the volcanic source. (b) Best-fitting solutions and confidence levels for the slip rate of the onshore HFF segment (N124|$^\circ$| E) versus the locking depth of the HFF. The best-fit model parameters are highlighted in orange. The solid and dashed magenta lines delineate the 68 and 95 per cent confidence levels, respectively.

Figure 5.

(a) Comparison between the best-fitting model (red) and GNSS velocities (black), relative to stable North America. (b) Residuals between the horizontal GNSS velocities and the model predictions. Velocities and residuals are depicted with 95 per cent uncertainty ellipses. Dashed lines represent the model dislocation segments along the plate boundary and the green star indicates the location of the Mogi source. Other features are the same as in Fig. 1. The inset marked by the dashed green square shows more details for the stations located at the Tjörnes peninsula.

4.2 Optimization routine and uncertainty estimation

To estimate the eight model parameters, we used a two-step optimization approach (Metzger et al. 2011, 2013; Metzger & Jónsson 2014). In the first step, we employed a Monte Carlo-type simulated annealing method (Cervelli et al. 2001) where the entire model space was thoroughly scanned to find the global minimum. In the second step, we used the solution obtained from simulated annealing as a prior for a derivative-based optimization routine to constrain the set of model parameters that best fit the GNSS velocity field. We used the velocity uncertainties as data weights in the optimization and ran the entire optimization five times to ensure reproducible results.

We then estimated the uncertainties of the model parameters using the randomize-then-optimize approach described by Metzger et al. (2011). In this approach, Gaussian random errors, scaled by the velocity uncertainties, are added to the observed velocities to slightly modify them and a new best-fitting solution is found. Repeating this error propagation process multiple times (1000 iterations) provides a statistical assessment of the model parameter uncertainties due to errors in the data; these uncertainties are comparable to those obtained using Bayesian inference (Metzger et al. 2013). However, neither this assessment nor Bayesian inference account for errors associated with the model assumptions themselves (e.g. simplified earth structure, straight faults, etc.).

4.3 Best-fit model parameters

We tested an additional model in which we allowed the locking depth of the onshore HFF (segment 8) to be independent of locking depth of the other transform segments. This refined model, henceforth called Model 2, has nine free parameters. The results of Model 1 and Model 2 are summarized in Table 2, along with three reference models from the literature Metzger et al. (2011) and Metzger & Jónsson (2014). We evaluate our models by calculating the root mean square error (RMSE). Model 1 and Model 2 have almost identical RMSEs: the east, north and up components have RMSEs of 0.57, 0.88, and 1.98 mm yr⁻¹, respectively, for both models. The differences in the RMSEs are smaller than 0.001 mm yr⁻¹, considerably smaller than the mean uncertainty of our velocity field, and are thus not statistically significant.

Table 2.

Best-fitting model parameters and uncertainty estimates for several plate-boundary models constrained by GNSS data. The difference between Models 1 and 2 is that the HFF is modelled as two segments with different locking depths in Model 2. Reference models R0, A1 and C1 are from Metzger & Jónsson (2014). Values in brackets are fixed. For model C1, the total plate motion and azimuth were set according to the MORVEL model values (DeMets et al. 2010). All models include a point source to account for broad uplift north of the Krafla volcano, with source depth from Metzger & Jónsson (2014).

Parameter	Model 1	Model 2	Reference R0	Reference A1	Reference C1
Locking depth ridge [km]	|$5.8 \pm 0.3$|	|$5.7\substack{+0.4 \\-0.3}$|	|$3.2\pm 0.2$|	|$4.5 \pm 0.3$|	|$3.1 \pm 0.2$|
Locking depth HFF 1 [km]	|$7.3\substack{+0.9 \\-0.7}$|	|$6.4\substack{+1.0 \\-0.6}$|	|$6.2\substack{+0.8 \\-0.7}$|	|$8.5\substack{+1.8 \\-1.4}$|	|$8.7\substack{+1.2 \\-0.9}$|
Locking depth HFF 2 [km]	–	|$9.4\pm 1.4$|	–	–	–
Total plate motion [mm/yr]	|$20.0\pm 0.2$|	|$19.9 \substack{+0.3 \\-0.2}$|	|$20.3\substack{+0.4 \\-0.3}$|	|$19.9\substack{+0.6 \\-0.5}$|	[18.4]
Azimuth [N|$^\circ$| E]	|$112.5 \pm 0.6$|	|$112.4\pm 0.6$|	|$109.4 \pm 0.7$|	|$111.9\substack{+1.0 \\-1.1}$|	[104.0]
Partial motion on HFF [%]	|$34.8\substack{+0.9 \\-1.0}$|	|$35.9\substack{+1.0 \\-1.2}$|	|$33.4 \pm 1.4$|	|$35.7\substack{+1.6 \\-1.7}$|	|$37.1\substack{+1.7 \\-2.1}$|
Krafla vol change [×106 m3]	|$125.8\substack{+7.9 \\-7.1}$|	|$125.2\substack{+8.1 \\-6.7}$|	–	|$203\substack{+42 \\-43}$|	|$114\substack{+53 \\-43}$|
Krafla source depth [km]	[21.5]	[21.5]	–	|$26.3\pm 0.25$|	|$25.6\substack{+4.4 \\-4.7}$|
Þeistareykir vol change [×106 m3]	–	–	[8.5]	|$20.5\substack{+7.7 \\-5.6}$|	|$48.0\substack{+13.7 \\-10.9}$|
Þeistareykir source depth [km]	–	–	[25]	|$8.6\substack{+1.4 \\-1.5}$|	|$11.9\substack{+1.5 \\-1.6}$|

The differences in the best-fit parameters between Model 1 and Model 2 are also marginal. In Model 2, the ridge locking depth is 0.1 km shallower. Additionally, the total plate motion rate is 0.1 mm yr⁻¹ smaller with an azimuth that is 0.1|$^\circ$| smaller, while the partial plate motion in the HFF is 1.1 per cent higher (Table 2). The Krafla volume change rate is lower by 0.6 × 10⁶ m³ yr⁻¹. All of these differences fall below the uncertainties of the estimated parameters. The best-fit estimate for the proportion of motion accommodated by the HFF increased by |$\sim$|1 per cent from 34.8 per cent in Model 1 to 35.9 per cent in Model 2.

Regarding the locking depths of the HFF, Model 1 predicts |$7.3\substack{+0.9 \\-0.7}$| km for the entire length of the HFF while Model 2 estimates |$6.4\substack{+1.0 \\-0.6}$| km for the offshore segment and |$9.4\pm 1.4$| km for the onshore segment. Although the estimated locking depth is deeper for the onshore HFF than for the offshore segment, it is difficult to draw any meaningful conclusions regarding this difference. The presence of a Mogi source particularly influences the locking depths, and in Model 2 this influence is restricted mainly to the onshore segment of the HFF.

F-tests have previously been used to assess whether the inclusion of additional plates significantly improves model fit in plate motion inversions (Stein & Gordon 1984). Similarly, we employed an F-test to assess whether there is a statistically significant difference between Model 1 and Model 2, which are nested models that differ by only one degree of freedom. We calculated the F-statistic using eq. (2):

$$\begin{eqnarray} F = \frac{(\mathrm{ RSS}_1 - \mathrm{ RSS}_2)/(p_2 - p_1)}{\mathrm{ RSS}_2/(n - p_2)} , \end{eqnarray}$$

(2)

where |$\mathrm{ RSS}_1$| and |$\mathrm{ RSS}_2$| are the residual sum of squares for Model 1 and Model 2, |$p_1 = 8$| and |$p_2 = 9$| are the number of free parameters (degrees of freedom) in Model 1 and Model 2, and |$n = 121$| is the total number of observations. The F-statistic assesses whether the increment in explained variance from Model 1 to Model 2 justifies the added complexity. The F-test indicates that Model 2 does not provide a statistically significant improvement in fit over the simpler Model 1 at the 95 per cent or even 90 per cent confidence levels. Additionally, Model 2 introduced greater uncertainty in estimating the HFF locking depth compared to other models. Consequently, we will focus on the results of Model 1.

We compared the best-fitting model parameters of our models with those from three reference models which use as input only GNSS data (Metzger et al. 2013; Metzger & Jónsson 2014), Table 2. Reference model R0 corresponds to the reference model in Metzger & Jónsson (2014), which includes one Mogi source to account for inflation at Þeistareykir at a fixed depth and volume change rate. Reference model A1 from Metzger & Jónsson (2014) includes two Mogi sources, one to represent the broad uplift observed north of Krafla, and one for the inflation at Þeistareykir during their study period. Reference model C1, also by Metzger & Jónsson (2014), has a fixed plate motion of 18.4 mm yr⁻¹ at N104|$^\circ$| E, based on (MORVEL DeMets et al. 2010). Further details of the construction of these models can be found in Metzger & Jónsson (2014).

The uncertainties of the estimates for total plate motion, azimuth and partial motion on the HFF are 15 per cent to 60 per cent lower in Model 1 compared to the uncertainties from the reference models R0, A1 and C1. These reductions reflect the decrease in the uncertainty of our velocity field, which can be attributed to our extended time-series. However, it is important to note that these uncertainties are model-specific, and alterations in the model setup or input data, such as excluding the vertical component or changing the number and location of the Mogi sources can influence the parameters beyond what is indicated by the estimated uncertainties (see Fig. S9 in the Supporting Information). This impact is particularly evident in the ridge locking depth parameter, where Model 1 indicates a depth of |$5.8 \pm 0.3$| km, compared to |$3.2 \pm 0.2$| km in model R0, |$4.5 \pm 0.3$| km in model A1 and |$3.1 \pm 0.2$| km in model C1. These estimates are not within the uncertainties for the ridge locking depth of the other models, and the uncertainty for the ridge locking depth in Model 1 is 50 per cent higher compared to models R0 and C1. The HFF locking depth estimate in Model 1 falls within uncertainties of the reference models, being deeper than in model R0 but shallower than in models A1 and C1. The uncertainty for the HFF locking depth was reduced by 24 per cent to 50 per cent in Model 1 with respect to models A1 and C1, but showed a slight increase of |$\sim$|7 per cent compared to model R0. This evidences the sensitivity of the locking depth parameter to the model setup.

Most model parameters exhibit minimal correlation with each other, as evidenced by correlation coefficients below 0.3 (Fig. 4a). The highest correlation coefficients observed include 0.50 between the total plate velocity and the ridge locking depth, 0.53 between the total plate velocity and the HFF locking depth, and 0.74 between the ridge locking depth and the volume change rate of the Krafla Mogi source. The histograms displayed in Fig. 4(a) generally show a normal distribution, with the exception of the histogram for the HFF locking depth, which is slightly skewed to the right, favouring locking depths shallower than 9 km.

Model 1 yields a total plate motion of |$20.0\pm 0.2$| mm yr⁻¹ with an azimuth of |$112.5^\circ \pm 0.6$|⁠, of which |$34.8\substack{+0.9 \\-1.0}$| per cent is accommodated by the HFF. This corresponds to a slip rate of |$7.0\pm 0.2$| mm yr⁻¹ on the offshore segment of the HFF (116|$^\circ$| N), and |$6.8\pm 0.2$| mm yr⁻¹ on the onshore segment (N124|$^\circ$| E). The estimated locking depth for the HFF in Model 1 is |$7.3\substack{+0.9 \\-0.7}$| km. A scatterplot illustrating the relationship between the HFF slip rate and the HFF locking depth indicates insignificant trade-off between these parameters (Fig. 4b).

5 INTERSEISMIC DEFORMATION

The agreement between the observed and predicted GNSS velocities is generally good (Fig. 5 and Tables S1 and S2 in the Supporting Information). The most notable discrepancies between observed and modelled velocities occur towards the NVZ and GOR. On the Eurasian plate, the model overpredicts the southward motion of stations south of 66|$^\circ$|N; for instance, station AUSB has a residual of −2.5 mm yr⁻¹ on the north component. This residual pattern could be attributed to unmodelled GIA, which is stronger in the south, closer to the Vatnajökull glacier (Cao et al. 2023).

On the Tjörnes block, the observed velocities at the offshore stations MANA and GMEY, near the GOR segment, are directed more to the SSE than the predicted velocities, resulting in south and southwestwards residuals (Fig. 5). Several factors may contribute to the misfit at these stations, including earthquake swarms, volcanic deformation within the GOR, structural complexities of the GOR itself that are not accounted for in our model, and the limited ability to resolve kinematics and fault parameters in this area due to a lack of offshore geodetic data. The remaining residuals on the Tjörnes peninsula are smaller. Previous studies hinted at a clockwise rotation of the residuals on the peninsula (Metzger et al. 2011). In our results, the residuals by the western coast of the Tjörnes peninsula, from station HOTJ all the way to station KVIS and HEDI and stations to the southeast suggest a slight clockwise rotation (Fig. 5b). Within the Þeistareykir volcanic system, the northernmost residuals (ARHO, BREV, INGL) point towards the SSE, while residuals farther south (BANG, AUDB, HEHO, GRFJ) point towards the SSW and are larger in magnitude (Fig. 5b). This pattern could be attributed to internal deformation within the block that we cannot resolve with our current model. The southernmost and easternmost residuals (BLAS, GAKE, HOLL) do not follow this rotation and are likely influenced by deformation related to the HFF and NVZ segments (Fig. 5b).

On the North-American plate, the southernmost residuals around Akureyri (AKUR, SPOR, GODA, BADA) exhibit a northward trend, while the northernmost residuals, east of Skagafjörður, (SELV, HRAC), tend to point south. Northward residuals at stations such as BADA and GODA are in the order of 2.5 mm yr⁻¹, while southward residuals at stations on the North-American plate such as SELV and HRAC are approximately 1.8 mm yr⁻¹. These observations may also correspond to GIA, similar to the residuals in the southeastern part of the network. Southeast of the HFF, there is generally good velocity agreement in terms of azimuth, although a few stations are slightly under or overestimated (Fig. 5b).

Our backslip model based on block rotations does not produce strong vertical displacements. Vertical motion is primarily induced by the introduction of the point-pressure source and the opening component along the segments. Since our observations were not corrected for GIA, our vertical residuals likely reflect the effects of GIA. The residuals for the vertical component are consistently small towards the northeast of the study area in the Eurasian plate and along the onshore HFF, particularly within the Þeistaykir volcanic system. This consistency may indicate that the combined effects of the opening component in the HFF segment and the point-pressure source are accurately captured in this area (Fig. 6). Stations closer to the Mogi source within the Krafla volcanic system still show uplift in their residuals with rates of around 3 mm yr⁻¹ (Fig. 6). Uplift is also prominent in the south and southwest of the studyarea is attributed to GIA. The sub region, with the highest rates close to station ADAL. Towards the northwest, the vertical residuals indicate subsidence in the order of 2 mm yr⁻¹ at station HRAC, 3.6 mm yr⁻¹ at station GJOG, and 4 mm yr⁻¹ at station GMEY (Fig. 6). The unpredicted uplift in the south of the study area is attributed to GIA. The subsidence in the northwest is also likely to be associated with GIA, although the current GIA models do not predict subsidence along the coast. Some of the subsidence could also be related to crustal thinning at the divergent plate boundary.

Figure 6.

Residuals between observed and modelled velocities for the up component. The data points were interpolated by kriging to produce a map of residual uplift (red) and subsidence (blue). The dashed brown lines represent the model dislocation segments, and the green star indicates the location of the Mogi source used in the modelling. Onshore fissure swarms are highlighted in yellow, and central volcanoes are outlined in orange.

To analyse the interseismic rate gradient across the TFZ, we plotted the fault-parallel (N116|$^\circ$| E) velocities, and the velocity estimates from Model 1, for 11 cGNSS stations and 26 eGNSS stations from the southwest towards the northeast across the TFZ (Fig. 7). Approaching the HFF from the southwest, the fault-parallel velocity starts to increase from values close to zero (stable North America) to values just below 5 mm yr⁻¹. In the middle of the Tjörnes block, 20 km northeast of the HFF, velocities reach |$\sim$|9 mm yr⁻¹. The velocities continue to increase more steeply as the GOR is approached, but once crossed, the velocities tend to stabilize. Sites on the Eurasian plate show velocities just below 18 mm yr⁻¹. Finally, although affected by the overlapping deformation from the GOR, the S-shaped profile across the HFF clearly shows that the fault is locked.

Figure 7.

Fault-parallel (N116|$^\circ$| E) GNSS velocities (depicted in dark blue with grey uncertainty bars) on a profile across the TFZ, perpendicular to the HFF (purple rectangle in Fig. 3). The shaded grey area indicates the upper and lower boundaries of the best-fitting model prediction, plotted in cyan. Vertical dashed red lines mark the locations of the HFF and GOR. Sites southwest of the HFF are situated on the North-American plate, sites between the HFF and the GOR are on the Tjörnes Block, and sites northeast of the GOR are on the Eurasian plate. NA: North America, TB: Tjörnes Block, EU: Eurasia.

The data fit of our preferred velocity model is best within the Tjörnes block and close to the model segments, while the rates far into the Eurasian and North-American plates are slightly over- and underestimated, respectively. The decrease in the amount of fault-parallel velocity around 40 km southwest of the HFF, particularly, evident for station ENBL, was previously attributed to uplift at Þeistareykir (Metzger et al. 2011). Given our longer time-series and results, we believe this decrease in velocity might be more strongly related to GIA. Station LEIT, to the northeast of the HFF, shows a mean fault-parallel velocity sligthly lower than nearby stations, but it still within the boundaries of the model prediction considering its uncertainty. The values for station LEIT are likely biased, as they were estimated from only two measurements.

Earthquake locations south of the HFF appear to form several N–S or NNE–SSW lines, which have been interpreted as likely source zones of past large earthquakes (Stefánsson et al. 2008). However, profiles of the fault-parallel velocities in Tröllaskagi and Flateyjarskagi do not show fault-parallel deformation across the DL (Fig. S9 in the Supporting Information), although the DL appears to mark the southern limit of the interseismic deformation.

5.1 Volcanic signals

Volcanic unrest in the northern volcanic zone affects geodetic observations of secular plate spreading. The broad uplift north of the Krafla central volcano, observed since 1993, can be attributed to magma accumulation at a depth of approximately 21 km (Sturkell & Sigmundsson 2000; de Zeeuw-van Dalfsen et al. 2004), or to viscoelastic relaxation following the Krafla rifting episode (Ali et al. 2014). This broad signal was modelled as a point-pressure source at a fixed depth, following Metzger & Jónsson (2014) and Drouin et al. (2017a). This inflation produces north and northwestwards displacements in the Tjörnes block, and somewhat decreases the ESE rates representative of the plate motion. In a North-American-fixed reference frame, the inflation increases the eastwards rate on stations in the Eurasian plate, while the positions of stations on the North-American plate would be pushed towards the west and WNW.

Other volcanic signals, such as the deflation inside the Krafla central volcano, observed since 1989 and attributed to cooling and contraction (de Zeeuw-van Dalfsen et al. 2004), were not included in our model as there is insufficient data to properly resolve the overlapping deformation signals. The source representing inflation at Þeistareykir (Metzger et al. 2011, 2013; Metzger & Jónsson 2014) was also omitted because Drouin et al. (2017a) reported subsidence at Þeistareykir following the 2007–2009 inflation period, and Drouin et al. (2020) and Dawson (2022) indicated that the transients at Þeistareykir were short-lived and likely do not contribute significantly to present-day stress changes.

6 DISCUSSION

The extension of the time-series by more than a decade compared to previous studies (Metzger et al. 2011, 2013; Metzger & Jónsson 2014) yields a clear reduction in the uncertainties of our GNSS velocities, even when accounting for temporally correlated noise. Our trajectory modelling reveals that short transient signals of volcanic origin, which were significant in the period studied by Metzger et al. (2011), are of secondary importance to our 20+ year-long time-series. However, some stations, for example, MYVA (Fig. S3 in the Supporting Information), continued to exhibit transient deformation throughout our study period and were consequently excluded from the modelling.

Our simple trajectory modelling also reveals residuals of concave shape in the time-series, indicative of accelerated motion, particularly in the vertical component. In our study area, the accelerations are mostly negative, indicating accelerating subsidence for stations moving downwards. The consistency in the acceleration terms may suggest a connection to accelerating GIA, as reported in previous studies (Compton et al. 2015; Cao et al. 2023). We found that incorporating a second-degree polynomial into the trajectory model significantly improved the fit of our time-series. However, the second-degree polynomial term was small and did not significantly affect the estimated linear trends, which remained virtually unchanged.

Although we demonstrated that the contribution of accelerating GIA to the GNSS velocities in our study area is negligible, a linear contribution from GIA might still significantly affect our observations, as suggested by the pattern of subsidence towards the northwest and uplift towards the southeast in our horizontal residuals Fig. 5. We did not correct for GIA, considering the current corrections are below 1 mm in the horizontal component and a few millimeters in the vertical component around Krafla (Drouin et al. 2017a), and would only decrease farther north and northeast. Our residuals are a little higher (⁠|$\sim$|2.5 mm yr⁻¹) than what we might expect for a horizontal GIA correction term. Additionally, the correction methodology using the current GIA model (Auriac 2014) does not reproduce either the subsidence or the southwards motion in the northwest of our study area, which we also associate with GIA. It is important, however, to be cautious when interpreting these results as they may not be solely related to GIA.

We tested estimating an independent locking depth for the onshore segment of the HFF (Model 2), as seismotectonic studies have suggested it might be deeper towards the west (Rögnvaldsson et al. 1998; Jónsdóttir et al. 2016). However, the difference between Model 1 and Model 2 was not statistically significant. At present, our ability to further resolve model parameters is limited due to the spatial distribution of our data, particularly the lack of constraints offshore. Therefore, we conclude that the HFF has a uniform locking depth of |$7.3\substack{+0.9 \\-0.7}$| km and accommodates |$34.8\substack{+0.9 \\-1.0}$| per cent of the total plate motion, which is |$20.0\pm 0.2$| mm yr⁻¹ at an azimuth of |$112.5^\circ \pm 0.6$|⁠. Our model parameters are within the uncertainty levels of previous studies, and—with the exception of the ridge locking depth—our values exhibit lower uncertainties compared to earlier work (Metzger et al. 2013; Metzger & Jónsson 2014), as summarized in Table 2. The locking depth parameters are particularly influenced by the model setup, notably by the Mogi sources and the vertical component of the velocity field (see Fig. S10 in the Supporting Information). The model might be performing poorly in the case of the ridge locking depth constrained by stations in the NVZ, as the NVZ is composed of several volcanic systems rather than our simple vertical line. The assumptions and simplifications of our model preclude us from further resolving the locking depths and slip rates of the offshore segments, where lack of data becomes the main limiting factor. Despite this, our model still provides a useful approximation in the vicinity of the HFF.

Assuming complete stress relaxation and steady stress accumulation at the locked interface of the HFF since the 1872 earthquakes, the accumulated slip deficit on the fault would be |$1.04\pm 0.02$| m. An earthquake releasing the full slip deficit and rupturing the entire approximately 110-km-long fault down to a constant (locking) depth of |$7.3\substack{+0.9 \\-0.7}$| km would correspond to M_w  |$\sim \!6.9$| on the moment magnitude scale (Hanks & Kanamori 1979). Improving the estimates of key fault parameters is crucial for seismic hazard assessment. The slip rate of the HFF is now better constrained than in previous studies; however, its locking depth remains poorly resolved. While our results are in line with recent tomography studies (Abril et al. 2021), further constraining the HFF fault locking depth remains of critical importance for future studies.

7 CONCLUSIONS

We generated time-series for 22 cGNSS stations from 2001 to late 2022 and for 99 eGNSS stations measured at least twice but typically at least five times between 1993 and 2022. The network has significantly improved spatial coverage compared to previous studies (Jouanne et al. 2006; Metzger et al. 2011, 2013; Metzger & Jónsson 2014), particularly towards the west and around the onshore portion of the HFF. From these time-series, we produced the most up-to-date GNSS velocity field for North Iceland (Fig. 3), which was then used to constrain a plate boundary deformation model for the study area and to estimate key fault parameters, such as the locking depth and slip rate of the HFF.

The interseismic velocity field indicates that the HFF is locked. To model the observed velocity field, we implemented a kinematic backslip model consisting of nine dislocation segments. The segments represent plate divergence in a double transform zone that encloses the Tjörnes block, combined with a point pressure source that mimics the broad volcanic uplift north of Krafla volcano. Our best-fitting model indicates that |$\sim \!$|35 per cent of the transform motion is focused on the HFF, with the rest accommodated by the GOR. The estimated relative plate motion of |$\sim \!$|20 mm yr⁻¹ with an azimuth of |$\sim \!$|N112|$^\circ$| E between the North-American and Eurasian plates is slightly higher and more SSE-directed than predicted by the MORVEL plate motion model (18 mm yr⁻¹, N104|$^\circ$| E). However, our results are in line with previous geodetic studies in the region (Metzger et al. 2011, 2013; Metzger & Jónsson 2014; Drouin et al. 2017a). From the model parameters, we find an average slip rate of |$6.9\pm 0.2$| mm yr⁻¹ on the HFF, with a locking depth of |$7.3\substack{+0.9 \\-0.7}$| km.

Except for the parameter describing ridge locking, all model parameter uncertainties are smaller than in earlier studies, due to the longer GNSS time-series, and improved network coverage. The HFF slip rate and locking depth remained stable and well-constrained within the uncertainties of the different models evaluated. However, further refinement of the estimated parameters remains somewhat elusive, mainly due to the lack of offshore data. However, our results indicate that the strain accumulated at the HFF since the last major earthquakes of 1872 are equivalent to a magnitude |$\sim \!6.9$| earthquake.

FUNDING

This research was supported by the King Abdullah University of Science and Technology (KAUST), under Award Number BAS/1/1353-01-01

Acknowledgements

We thank Ulas Avsar, Teng Wang, Rishabh Dutta, Wenbin Xu, Jon Harrington, Ayrat Abdullin and Max Mai for taking part in the 2013 measurement campaign, Joel Ruch, Hannes Vasyura-Bathke, Tim Sonnemann and James Scott Berdahl for taking part in the 2016 campaign, Laura Parisi, Daniele Trippanera, Xing Li, Sturla Thengilsson, Laila Mai, Kári Daníel Alexandersson and Margrét Hrafnsdóttir for taking part in the 2019 campaign and Jochen Woessner, Matthieu Ribot, Adriano Nobile, Adrien Moulin, Sturla Thengilsson, Paul Mai, Margrét María Hallgrímsdóttir and Margrét Hrafnsdóttir for taking part in the 2022 campaign. We thank Guðmundur Valsson (National Land Survey of Iceland) for providing the ISNET GNSS data.

DATA AVAILABILITY

The bathymetry (https://doi.org/10.5285/f98b053b-0cbc-6c23-e053-6c86abc0af7b) is from the GEBCO 2023 Grid courtesy of GEBCO Compilation Group (2023). The digital elevation model, the fissure swarm fractures, central volcanoes and onshore fissure swarms are from the National Land Survey of Iceland and can be obtained from their website (https://www.lmi.is, last accessed October 2023). The offshore fractures and other features were mapped from multibeam bathymetry data, courtesy of Bryndís Brandsdóttir from the University of Iceland. The IMO seismic catalogue is from the Icelandic Meteorological Office and can be obtained from their website (https://www.vedur.is, last accessed October 2023). We provide a supplementary material file containing the observed and predicted velocities of the GNSS stations in North Iceland.