Present kinematics of the Tj ¨ornes Fracture Zone, North Iceland, from campaign and continuous GPS measurements

campaign GPS (EGPS) data points, which have been regularly observed since 1997. We extract the steady-state interseismic velocities within the TFZ by correcting the GPS data for volcanic inﬂation of Theistareykir— the westernmost volcano of the NVZ—using a model with a magma volume increase of 25 × 10 6 m 3 , constrained by InSAR time-series analysis results. The improved velocity ﬁeld based on 58 GPS stations conﬁrms the robustness of our previous model and allows to better constrain the free model parameters. For the HFF we ﬁnd a slightly shallower locking depth of ∼ 6.2 km and a slightly higher slip-rate of ∼ 6.8 mm yr − 1 that again result in the same seismic potential equivalent to a M w 6.8 earthquake. The much larger number of GPS velocities improves the statistically estimated model parameter uncertainties by a factor of two, when compared to our previous study, a result that we validate using Bayesian estimation.


I N T RO D U C T I O N
The Tjörnes Fracture Zone (TFZ) in North Iceland, is one of two transform zones in Iceland that can produce magnitude ∼7 strikeslip earthquakes with a recurrence interval of decades to centuries (Einarsson 1991). Due to its mostly offshore location and its complex geometry, consisting of at least two parallel transform structures, the detailed tectonics and kinematics of this plate boundary zone are still not well understood. One of the main lineaments in the fracture zone, the Húsavík-Flatey fault (HFF), has ruptured in major earthquakes in 1755 and 1872, which indicates that the interseismic period of the earthquake cycle might end soon. It is therefore important to estimate the seismic potential of this strike-slip fault and its implications for the inhabitants of Húsavík, a town located right on top of the fault. One way to assess the seismic potential of a locked fault is to estimate the dimensions of the locked fault plane (length and locking depth) and the accumulated slip deficit since the last major earthquake, derived from interseismic deformation rates (Wesnousky 1986).
Several studies have provided information about the slip rate and the locking depth of the HFF. An estimation of max. 60 km accumulated slip during the last 7-9 Myr from Saemundsson (1974) and Gudmundsson et al. (1993) results in an average slip rate of a max. 7-8 mm yr −1 . Rögnvaldsson et al. (1998) analysed relative locations of 1400 earthquakes in 60 swarms within the TFZ 1994-1997 with a vertical uncertainty <2 km and found that less than 10 per cent of all C The Authors 2012. Published by Oxford University Press on behalf of The Royal Astronomical Society. events occur deeper than 10 km, deducing a base of the brittle crust at 10 km. Velocities from the first three continuous GPS stations in North Iceland that were installed in [2001][2002] suggest that the HFF accommodates 40 per cent of the total plate motion . Using a plate spreading velocity of 18 mm yr −1 (MORVEL, DeMets et al. 2010) this corresponds to a slip rate of 7 mm yr −1 . Other slip rates derived from campaign GPS data are in the range of ∼5 mm yr −1 (Árnadóttir et al. 2009) to 8 mm yr −1 (Jouanne et al. 2006). In addition, Árnadóttir et al. (2009) used modelling to constrain the locking depth of the HFF at shallow 5 km, but the data used were sparse near the HFF and thus the locking depth was poorly constrained.
In our previous study, we presented data from 14 continuous GPS (CGPS) stations in North Iceland, most of which we installed in 2006 to monitor the surface deformation within the TFZ (Metzger et al. 2011). We modelled the kinematics of the plate boundary in North Iceland with an interseismic back-slip model based on nine dislocations in an elastic half-space plus a Mogi source representing inflation of the Theistareykir volcano (Fig. 1). Our aim was to estimate both the locking depth and the slip rate of the HFF. We found a slip rate of 6.6 ± 0.6 mm yr −1 and a surprisingly shallow locking depth of 6.3 +1.7 −1.2 km. Furthermore, the modelling revealed that the HFF accommodates only 34 ± 3 per cent of the total plate motion, while the majority of the transform motion must be taken up by the offshore Grímsey Oblique Rift (GOR), the other main lineament of the TFZ. Assuming complete stress unloading in the last major earthquakes in 1872 and steady stress accumulation since then, the seismic potential of a complete rupture of the HFF corresponds to a M w 6.8 ± 0.1 event (Metzger et al. 2011).
In this paper, we improve the kinematic model of the Tjörnes Fracture Zone by adding 44 campaign GPS (EGPS) observations to the 14 CGPS data to further constrain the model parameters of the TFZ. In addition, we analyse the time evolution of the transient volcanic uplift at Theistareykir using InSAR time-series analysis and remove this signal from the GPS data to obtain corrected interseismic deformation rates. Due to the enhanced network density near the onshore part of HFF and across the rift zone along the Krafla system we are able to better constrain all model parameters.
After an overview section of the tectonic setting in North Iceland we present a review of historical earthquakes in the region. We then document the acquisition and processing of the campaign GPS data, explain how we resolve the transient deformation signal using InSAR time-series analysis, and how we estimate the steady interseismic velocities from both the EGPS and CGPS time-series. The last part of this paper explains how we model the kinematics of the TFZ and derive best-fitting model parameters and their uncertainties using two different approaches of uncertainty estimation. Finally, the obtained results are discussed and compared with earlier published kinematic parameters.

T E C T O N I C S E T T I N G O F N O RT H I C E L A N D
Iceland is located on the mid-Atlantic Ridge where the Eurasian and North American Plate are separating at a rate of 18 mm yr −1 in a N104.5 • E direction, according to the MORVEL plate-motion model (DeMets et al. 2010). The plate boundary zone in Iceland is a few tens of kilometre wide, crossing the island from the southwest to the northeast, and is mostly defined by volcanic zones [i.e. the western, eastern and Northern Volcanic Zones (NVZ)] that accommodate the spreading motion within several parallel and overlapping volcanic systems ( Fig. 1). However the onshore part of the plate boundary is offset ∼100 km towards east, resulting in two transform zones at the southern and the northern shore of the island. The southern transform zone is called the South Iceland Seismic Zone and connects Reykjanes peninsula, a continuation of the Reykjanes Ridge, to the Eastern Volcanic Zone. It is completely onshore and consists of a series of parallel N-S trending strike-slip faults, sometimes referred to as bookshelf faulting (Einarsson 1991). The Northern transform zone is usually called the Tjörnes Fracture Zone (TFZ) and is of a completely different geometry. It consists of the 100-km-long HFF and a second lineament located northeast of the HFF, the Grímsey Oblique Rift (GOR), that both connect the NVZ to the Kolbeinsey Ridge (KR).
The plate boundary structures in the TFZ are primarily marked by fault surface traces and microseismicity. The HFF is a mostly off-shore strike-slip fault with a minor opening component. The WNW-oriented 90-km-long main fault strand extends from the Theistareykir fissure swarm to the Eyjafjar aráll Rift. This rift links the HFF and the Kolbeinsey Ridge, which is the northward continuation of the mid-Atlantic Ridge (Fig. 1). The HFF has been active for as long as 7-9 Myr with an estimated cumulative displacement of max. 60 km and is older than the GOR (Saemundsson 1974;Gudmundsson et al. 1993;Homberg et al. 2010). Fault surface traces between the coastal town Húsavík and the Theistareykir fissure swarm show several parallel branches and a slight bend at Húsavík leading to a NW-strike of the onshore part of the fault system. This change in strike results in an increasing opening component of the fault system, which is expressed in pull-apart structures (sag ponds) southeast of Húsavík. The southeastern part of the fault lacks microseismicity, particularly since the Krafla rifting-episode 1975-1984(Tryggvason 1980, 1984Björnsson 1985). One possible explanation for the low microseismicity is that the rifting episode released accumulated stress on the southeastern part of the HFF .
The second main lineament in the TFZ is the Grímsey Oblique Rift. It is similar to the Reykjanes Peninsula in southwest Iceland, consisting of a set of parallel, N-S oriented faults exhibiting both normal and strike-slip faulting and it is volcanically active (Rögnvaldsson et al. 1998). In our previous work we estimated that 66 ± 3 per cent of the total TFZ transform motion is accommodated by the GOR (Metzger et al. 2011) where as the rest is taken up by the HFF.
The NVZ is a 40 km wide rift zone that extends from Vatnajökull glacier in Southeast Iceland to the north coast and the TFZ, consisting of five volcanic systems (from NNW to SSE): Theistareykir, Krafla, Fremrinámar, Askja and Kverkfjöll (Fig. 1). It became active 7-9 Myr after an eastward jump of the plate boundary (Saemundsson 1974). The most recent rifting episode within the NVZ took place in the Krafla volcanic system in 1975-1984(Tryggvason 1980, 1984Björnsson 1985) with over 20 intrusive events and crustal widening amounting to an average of 5 m across the Krafla fissure swarm, corresponding to 275 yr of plate spreading (Tryggvason 1994). After 9 yr the rifting episode stopped, but the extensional pulse continued to propagate away from the rift axis and over time slowly decayed in amplitude (Foulger et al. 1992;Heki et al. 1993;Hofton & Foulger 1996). Fifteen years later the influence of the rifting episode had diminished (Völksen 2000). Other recent rifting episodes in the NVZ occurred in the Krafla volcanic system 1724-1729 and Askja volcanic system 1874-1875 (Sigurdsson & Sparks 1978). Two minor rifting events took place in 1618 and 1885 in Theistareykir, the westernmost volcanic system in the NVZ (Magnúsdóttir & Brandsdóttir 2011). Theistareykir has not erupted since 2500 BP (Karl Grönvold 2011, personal communication) but between 2007 and 2009 it had a period of inflation with a maximum uplift rate of ∼3 cm yr −1 in late 2008 (see Section 5.1).

H I S T O R I C A L E A RT H Q UA K E S I N T H E T JÖ R N E S F R A C T U R E Z O N E
Knowledge of large historical earthquakes provides vital information on the average length of the earthquake cycle. The first instrumentally derived magnitude and location estimations of Icelandic earthquakes were available in the beginning of the last century. For older events we have to rely on historical accounts, which are far from being complete, and treat the available information with caution. Information on historical earthquakes in North Iceland dates back to 1260, when a large earthquake took place close to Flatey island ( Fig. 1, Thoroddsen 1880Thoroddsen , 1925, but accounts from the following 3-4 centuries are limited (Thorgeirsson 2012). Below we list the largest known earthquakes in the TFZ during the past 300 yr with magnitude estimates as given in (Stefánsson et al. 2008).
In 1755 a M7 earthquake took place in Skjálfandi Bay, which was felt widely on Tjörnes Peninsula as well as on the two peninsulas west of Skjálfandi Bay and caused damage to many buildings, open cracks and rock fall along the fjords (Thoroddsen 1925). Another earthquake happened near the northwestern end of the HFF in 1838 (M6.5). The next two large earthquakes on the HFF occurred 6 hr apart in 1872 and the size of each has been estimated M6.5. The former event took place near Húsavík, causing fissures with up to 1 m of opening (Thoroddsen 1925, p. 400), but the latter occurred close to Flatey island. Other major events were the M6.2 earthquake close to the town of Dalvík in 1934 and a M7 earthquake in 1963 that was located ∼60 km northwest of Dalvík. Due to these two large events and some microseismicity Einarsson (1976) have suggested a third 'weak' Dalvík zone southwest of the HFF with several N-S trending faults, similar to the Grímsey Oblique Rift.
Known large earthquakes along the Grímsey Oblique Rift are a M6.3 earthquake in 1885 induced by the Theistareykir rifting event at that time, a M7 earthquake in 1910 midway between the island Grímsey andÖxarfjördur and a M6.2 earthquake in 1976 at the northern end of the Krafla fissure swarm. The last one mentioned was part of the initial phase of the Krafla rifting episode 1975-1984(Tryggvason 1980Björnsson 1985;Tryggvason 1984;Passarelli et al. 2012).
In order to estimate the seismic potential of the HFF we assume that during the last 250 years only three major earthquakes in 1755 and 1872 ruptured substantial parts of HFF. This indicates that the next big event might be due soon. We furthermore assume that the stress regime on the HFF was completely relaxed in the last two earthquakes of 1872 and that the slip deficit of the fault has steadily increased since then.

GPS DATA
The first GPS reference markers were installed in North Iceland in 1987 to study the post-rifting relaxation after the Krafla rifting episode 1975-1984(Foulger et al. 1992Heki et al. 1993). This network was remeasured several times in the following years (1990,1992,1995) to analyse the influence of the rifting episode and the return to steady-state deformation (Hofton & Foulger 1996;Völksen 2000. Another, partially overlapping network was installed in 1995 to study the TFZ and it has been remeasured completely five times (1997/1999/2002/2007/2010) and partially on several occasions (2000/2001/2005/2009/2011), that is, mostly the GPS markers close to Húsavík (Jouanne et al. 1999(Jouanne et al. , 2006. 15 new campaign GPS points were added in 2009-2011 to improve the GPS station density near the HFF. Today the TFZ GPS network consists of 14 continuous and 61 campaign GPS stations (Fig. 1, Table 3).

GPS campaigns 1997-2011
During GPS campaigns, each point was typically measured for two consecutive nights, except in 1997, when the observation time was somewhat shorter (∼24 hr). The GPS antennas were aligned to magnetic north in the earlier campaigns, whereas in the 2009 campaign and after the antennas were aligned with true north using a magnetic declination of −14 degrees. The instrument types and the number of points measured differed between campaigns, see summary information in Table 1. The publication by Jouanne et al. (2006) provides more information about the acquisition strategy of the 1997-2007 GPS data. In 2009 we installed four new markers on the profile across the HFF and in 2010-2011 we added 11 more data points that enhance the network density north and south of the  Table 2. Details of the GPS processing with the GAMIT-GLOBK v10.4 (Herring et al. 2010a,b,c) and Bernese V5.0 (Dach et al. 2007) software. on-land part of the HFF (Fig. 1). The short time-series of these 15 points are not included in the work we present here.

Data processing with GAMIT-GLOBK and BERNESE
The entire campaign GPS data set 1997-2010 was processed independently with GAMIT-GLOBK V10.4 (Herring et al. 2010a,b,c) and Bernese V5.0 (Dach et al. 2007) in order to the test robustness of the resulting GPS velocities. The details of the processing methods are given in Table 2. Both software packages constrain daily network solutions using a well-chosen set of reference (IGS) stations. In both cases we included IGS stations as well as CGPS stations in North Iceland (when available) in the data analysis and we corrected for offsets in the time-series of the continuous GPS stations due to earthquakes or antenna changes. We used the same ocean-loading model (FES2004), absolute antenna phase centre models, and the Saastamoinen model to correct for the atmospheric delay. The dry/wet delay was corrected with the GMF (in GAMIT-GLOBK) and NIELL (in BERNESE) mapping function. In the GAMIT processing we used a slightly higher cut-off-angle (10 • ) than in the BERNESE processing (5 • ). In GPS campaigns before 2009 the antennas were aligned with magnetic north instead of true north. Since the antenna phase centre usually does not coincide with the geometrical antenna centre and can deviate quite substantially, a correction to that shift should improve the results. We corrected for this shift in the GAMIT processing but not in the BERNESE processing. In addition, during GAMIT processing we already removed data outliers in the EGPS data, whereas with BERNESE we processed all the EGPS data without restrictions. Therefore, the BERNESE results needed some post-processing visual quality control and filtering (see last part of Section 5.3).

Continuous GPS data 2006-2011
In addition to the episodic GPS (EGPS) data 1997-2010 we use velocities of 14 continuous GPS (CGPS) stations in North Iceland plus one additional station in east Iceland (HEID). The CGPS network was densified significantly from 4 to 14 stations in 2006 as we described in our previous study (Metzger et al. 2011). In that study we estimated the velocities with a superposition of a linear trend and a sinusoidal oscillation term. Here we have a better constrained estimation of the interseismic deformation rates due to an additional year of data and an improved correction of the transient volcanic signal of Theistareykir volcano. We explain in the following section how we derive the interseismic deformation rates for all the GPS sites.

E S T I M AT I O N O F I N T E R S E I S M I C G P S V E L O C I T I E S
Tectonic plates move steadily over millions of years but at plate boundaries the deformation is highly episodic. Most plate boundaries are locked for long periods of time and then suddenly become active in discrete events such as earthquakes or rifting episodes. The earthquake deformation cycle consists of slow but steady interseismic, sudden co-seismic, and transient post-seismic deformation. In rift zones as the NVZ, rifting episodes cover the average spreading of hundreds of years within relatively short rifting episodes that typically last a few years, followed by decaying post-riftingrelaxation. In addition, volcanic centres within rift zones can go through periods of inflation or deflation, which may cause transient signals in GPS time-series. In this paper, we want to extract the steady-state deformation rates within and around the TFZ, in particular the interseismic deformation rates near the HFF to constrain model parameters for the fault, so we can estimate the slip deficit of the locked HFF accumulated since the last major earthquake. This means that we have to identify and eliminate possible transient signals within the GPS time-series to isolate the linear interseismic deformation rates of the TFZ.
The Krafla rifting episode ended 25 yr ago and its effect on GPS velocities appears to have diminished before 1997 (Völksen 2000). We therefore assume that the plate-boundary deformation rates during the observation time can be approximated as linear. On the other hand, inflation at Theistareykir central volcano 2007-2009 caused a transient deformation signal at neighbouring GPS sites that we have to take into account (Metzger et al. 2011). In addition to tectonic signals, annual oscillation signals are visible within the continuous GPS time-series. Fig. 2 schematically shows the different components apparent in the acquired GPS time-series of the TFZ. We model the CGPS time-series using the following function: where A + Bt represents the linear displacement rate, C the amplitude and φ the phase shift of annual oscillations and D the amplitude, t c the central time and t d the curvature factor of the transient volcanic inflation at Theistareykir central volcano. To extract the steady interseismic velocities, we have to first correct all the GPS data for the volcanic transient signal (Section 5.1) and-in case of CGPS data-annual oscillations (Section 5.2).

Uplift at Theistareykir volcano
In our previous work we detected uplift in ENVISAT interferograms at Theistareykir central volcano with a maximum uplift rate of 3 cm yr −1 between 2007 and 2008 (Metzger et al. 2011). Using both an ascending and a descending interferogram we were able to fit the spatial pattern of the inflation well with a Mogi source (Mogi 1958) and thus constrain the source location. However, the limited ENVISAT catalogue does not contain enough data to analyse well the temporal evolution of the volcanic inflation. The CGPS timeseries do not provide useful information about the transient either, because they started in the midst of the inflation period and the CGPS stations are located too far away from the inflation centre. We We produced a number of interferograms out of 38 descending orbit scenes of the ERS-1 and 2 satellites spanning the timeperiod from 1992 to 2010 using the GAMMA software (GAMMA v1.0 2006). We then performed time-series analysis of the unwrapped interferograms using the π -RATE software package (Poly-Interferogram Rate And Time-series Estimator, 2009) that was developed by Biggs et al. (2007), Elliot et al. (2008) and Wang et al. (2009). The expected steady-state deformation, that is, due to plate spreading across the plate boundary and interseismic strain accumulation near the HFF, covers the whole scene and would be partially removed by the software during the correction of orbital errors. We preserve the cross-boundary signal by removing predicted lineof-sight (LOS) displacements from the interferograms before the orbital and topographical error correction and then add these predicted displacements back in again. The predicted displacements were derived from the deformation model of Metzger et al. (2011). The InSAR time-series processing steps are explained in more detail in Metzger & Jónsson (in preparation).
We extract an 800 × 800 m area from the InSAR time-series results, centred at Theistareykir volcano, and calculate the mean LOS displacement for each image acquisition date between 1992 and 2010 (Fig. 3). We find a linear rate between 1997 and 2007 that is caused by steady-state plate-boundary processes. Sudden uplift of ∼78 mm in LOS occurred in 2007, that seems to have ended in late 2008. Before 1997 we can see another transient which was most probably caused by deep magma accumulation at the Krafla volcanic system as reported by de Zeeuw-van Dalfsen et al. (2004). We fit an inverse tangens functional to the temporal evolution of LOS displacements after 1997 (Fig. 3)

Annual oscillations within CGPS time-series
The CGPS time-series show substantial annual oscillations, particularly in the vertical component (Grapenthin et al. 2006). After removing the volcanic uplift signal from the continuous time-series (Section 5.1) we then fit the following function, with a linear velocity term A + Bt and an oscillation term C to each single component and station. Before the final non-linear optimization run, outliers were removed in two separate stages: First, all data points with a standard error three times larger than the mean error were dismissed. This affected only a couple of data points. After a first optimization run we excluded also all data points with a misfit three times larger than the mean misfit. The amount of excluded data for each station-component was never larger than 4 per cent. With a second optimization run we extracted the CGPS velocities of the cleaned data set. The resulting phase shift and amplitude parameters are plotted in Fig. 4 and show clearly, how the vertical component contains the strongest oscillation signal with a maximum value in winter time (φ = 0). For the east and north component the amplitude is much smaller and thus not as nicely aligned with its mean value. However it is interesting to see that the mean phase shift for the east and north components differ from the vertical component with a maximum in October and September, respectively. The average rms

Estimation of the EGPS velocities
The campaign data have always been acquired during summers so we ignore possible influence from annual oscillations. After correcting the EGPS time-series for the transient volcanic signal at Theistareykir we estimate the velocity for each station-component. For the EGPS data points we get a mean standard deviation from the assumed linear model of 3.4 mm for the horizontal and 22.3 mm for the vertical components (Fig. 5). The velocity variance was scaled by 1/T 2 , where T is the total duration of each time-series Metzger et al. 2011). The EGPS data set was analysed independently with both the BERNESE V5.0 and the GAMIT-GLOBK software packages as discussed in Section 4.2. To estimate the interseismic velocities of the EGPS stations we used SNX-files obtained from the BERNESE processing and the time-series resulting from the GLOBK analysis. Within GAMIT processing the data were already cleaned from potential outliers. This was not the case for the BERNESE data set where all data were processed and we therefore applied a visual quality control in the velocity estimation process to eliminate outliers in each time-series.
The resulting BERNESE and GAMIT velocity fields match well within uncertainties. We find that the velocities derived from GAMIT point slightly more towards the south (mean shift: 0.6 ± 0.6 mm yr −1 ) and the east (0.3 ± 0.6 mm yr −1 ) than the BERNESE solutions. The vertical components of GAMIT also show less upward motion (2.1 ± 1.5 mm yr −1 ). Although this indicates a slight shift between the two network solutions we consider it as neglectable. The standard deviations of the shifts are larger than the shift itself, and what is more important, we account for a possible reference frame offset in our modelling in any case. The EGPS velocities we used in the modelling (Section 6) were obtained using the GAMIT-GLOBK software (Table 3).

Interseismic deformation in North Iceland
The estimated velocities and uncertainties for both the CGPS and EGPS stations are plotted in Fig. 6 in a reference frame originally based on stable North American Plate (MORVEL, DeMets et al. 2010), but modified with a small offset estimated in the modelling (Section 6). This correction was necessary for the modelling as GPS points located on the North American Plate, that is, southwest of  The best-fitting model parameters in comparison to the results of previous studies. The fault-parallel slip rate along the two HFF segments (Fig. 1) is derived from the partial transform motion parameter and the parameter describing amplitude and azimuth of the relative plate motion. The studies denoted with an asterisk ( * ) share parts of the data set.  1997-2011 2006-2010 2001-2004 1993-2004 1997-2002 1993-1999 1994-1998 7- Rögnvaldsson et al. (1998); h Saemundsson (1974).
Húsavík, are not entirely stable in the North American reference frame and show a 2-4 mm yr −1 motion to the northwest. A similar motion pattern has been seen in previous studies (e.g. Árnadóttir et al. 2009) and it seems to apply to all GPS stations in Northwest Iceland, indicating a local reference frame problem. We find that the TFZ covers the full plate motion of 18 mm yr −1 and that the deformation gradient across the NVZ is particularly strong (Fig. 6). All data vectors are more or less perpendicular to the orientation of the NVZ, except for stations close to Krafla central volcano that are influenced by local subsidence. Station velocities south of Skjálfandi Bay towards the northern tip of Tjörnes Peninsula gradually increase, indicating a locked HFF (see Fig. 1 for geographic locations).
Due to potential inconsistencies of the antenna height measurements the vertical deformation rates of EGPS data have to be interpreted with caution. North of Krafla central volcano a broad uplift of up to 6 mm yr −1 is apparent (Fig. 6). This signal coincides with uplift seen in InSAR data during -1999(de Zeeuw-van Dalfsen et al. 2004) that was thought to be caused by deep magma accumulation north of Krafla central volcano. At Krafla (relative) subsidence is visible, which is still ongoing but has been slowly decaying since the end of the Krafla rifting episode (Sturkell et al. 2008).

M O D E L L I N G
We describe the surface deformation of the TFZ as it has been observed by episodic and continuous GPS measurements during the last 14 yr with an interseismic back-slip model. This model consists of a set of nine planar plate-boundary segments with a fixed (Cartesian) geometry in an elastic half-space, as described by Metzger et al. (2011). The location of the plate boundary segments is shown in Fig. 1. The boundary follows roughly the Krafla fissure swarm in the south, then separates into two subparallel discontinuities along the HFF1 and 2 and the Grímsey Oblique Rift (GOR) that again reunite north of the TFZ at the Kolbeinsey Ridge. The plate-boundary deformation is described by superimposing reverse slip ('back-slip') on the locked part of the plate boundary onto rigid plate motion. The model does not allow for any rotation, because the slip on each model segment is uniform and is defined by the overall relative plate motion. Unlike in our previous study we do not include volcanic inflation at Theistareykir in the modelling, as we already eliminated this transient signal in the time-series analysis (Section 5). Otherwise, the model parameter optimization and error estimation procedures follow Metzger et al. (2011). The model parameter uncertainties are estimated by stochastically propagating the data errors through the modelling and do therefore not include the impact of model assumptions and simplifications. We here refer to this method as 'error propagation' and compare it to a second, independent uncertainty estimation, based on Bayesian estimation (Section 6.1).
The model segments are described by 10 parameters each, but many of them are constrained due to (1) the fixed location and vertical dip, (2) only two types of locking depths for the strike-slip segments (HFF1/2 and GOR in Fig. 1) and 'rift-type' segments (all other segments) and (3) no dip-slip motion. The key parameter we solve for is the locking depth of the HFF (Table 4). This value, along with the fault's length, determines the size of the locked HFF plane. Additional parameters are the locking depth of the ridge segments, the partitioning of transfer motion among the two lineaments HFF and GOR and the magnitude and azimuth of the overall plate spreading motion (Table 4). Two auxiliary parameters allow for a small shift of the North American reference frame into a model frame, which assumes a stable reference point southwest of the HFF. Together, the spreading vector and the partitioning of motion between the HFF and GOR define the slip-rate on the HFF. The slip-rate provides information about the stressing rate on a fault plane since the last large earthquake and gives, together with the locking-depth, an estimation of the seismic potential of this fault plane.
We include horizontal velocity components of all the CGPS and EGPS stations in the model parameter optimization (Table 3, except for one CGPS and two EGPS stations located near Krafla central volcano (MYVA/KRAF/HVIT). These stations are influenced by local subsidence at Krafla, which we do not account for in our model, because it is far from the HFF and does not influence nearfault velocities.
We find a locking depth for the HFF (and the GOR) of 6.2 and 3.2 km for the ridge segments. The total relative plate motion of 20.3 mm yr −1 with an azimuth of 109.4 • E is separated between the HFF and the GOR in a ratio of 33/67 per cent. The resulting data fit obtained by the best-fitting model parameters (Table 4) is in general very good (Fig. 8, Table 3), except for the area north of Krafla central volcano where the steady-state interseismic deformation is slightly modified by uplift and extension. This is presumably due to deep magma accumulation during -1999(de Zeeuw-van Dalfsen et al. 2004).

Validation of the uncertainty estimation using a Bayesian estimation
We estimated the uncertainties of the best-fitting model parameters (Table 4) by analysing outcomes of 10 000 optimization runs with slightly modified ('noisy') input data. This method of estimating parameter uncertainties only considers the errors in the input data, but does not include errors from the modelling procedure or due to model simplifications. The resulting model parameter uncertainties should therefore be regarded as minimum uncertainties. Now, we validate this method of error propagation by applying Bayesian estimation, which provides a posterior probability distribution over the model parameters given the recorded data and serves as an independent assessment of the model.
We assume a linear M-dimensional model space M, and a linear D-dimensional data space D. The forward operator between these two spaces g, is assumed to be only mildly non-linear. The recorded measurements d and model parameters m are assumed to be realizations of the random variables D and M such that where is a realization of stochastic noise. The posterior density in the model space is calculated according to Bayes' formula, here given in the form of Tarantola (2005) gives a measure of how good a model m is for explaining the data. ρ D is the prior information on the data, and θ (d|m) represents the correlation between the data and model parameters. Assuming the theoretical relationship between model parameters and data to be exact, θ (d|m) = δ [d − g(m)], allows to solve eq. (5). From the independence of and m in eq. (4), and from the assumption of Gaussian errors it follows that (Tarantola 2005) This density σ M must then be evaluated by a Markov Chain Monte Carlo (MCMC) method due to the high dimensionality of model space M. An MCMC algorithm is an algorithm for constructing a Markov chain with an equilibrium distribution equal to a given probability density function. The Metropolis-Hastings (M-H) algorithm (Hastings 1970) is an MCMC algorithm that picks the following state from a proposal distribution that is simpler than the sampled distribution, but uses a condition for rejecting unlikely states with greater probability to more likely ones. To find suitable parameters for this proposal distribution we first used a parameterfree MCMC algorithm called the Gibbs sampler described in Geman & Geman (1984). The Gibbs sampler iteratively samples the conditional distributions of each variable, making it considerably slower for the problem at hand that lacks simple conditional distributions.
The M-H algorithm was started as multiple chains from arbitrary non-zero starting points in the model space. The first third of every chain is discarded to ensure convergence. The remaining samples are treated for autocorrelation by thinning, where only every τ th sample was picked, so that samples τ steps apart are uncorrelated. τ was determined by the Geyer IMSE heuristic described in Geyer (1992). The remaining samples were considered representative of the posterior density (eq. 7).
The resulting marginal distributions of the model parameters from the Bayesian estimation are in a perfect agreement with the uncertainty estimation obtained by the error propagation (Fig. 7). Since one Bayesian estimation needs much less calculation time than obtaining a significant statistic from the error propagation method, we compare the 10 000 optimization runs of the error propagation method with one million of samples from the Bayesian estimation. This explains the difference in smoothness of the marginal distributions in Fig. 7. The 68 and 95 per cent-confidence levels do not deviate from each other more than 2 per cent for any parameters. This result is reassuring and confirms the validity of the error propagation method for determining the part of the model parameter uncertainties that is caused by errors in the input data.

D I S C U S S I O N
The kinematic model presented in this paper is based on GPS timeseries from 14 continuous GPS stations running since 2006 and 44 GPS markers that have been remeasured at least five times since 1997. We can compare the obtained result directly to our previous results (Metzger et al. 2011) where we used only the 14 CGPS data points (Table 4). The CGPS stations lie primarily on a profile across the two lineaments HFF and GOR and are sparse near the volcanic rift zone, which previously resulted in a poorly constrained ridge locking depth. The campaign GPS observations complement the CGPS network nicely and fill the gaps in the NVZ and on the Flateyarskagi peninsula southwest of the HFF (Fig. 1). This is the reason why the uncertainties of all our model parameters are at least 50 per cent smaller (Table 4) and why we get significant changes for the locking depth of the ridge segments (from 4.8 to 3.2 km) and for the azimuth of plate motion (from 115 • E to 109.4 • E), which is now closer to the MORVEL plate motion azimuth of 105 • E. The full plate motion increases again slightly from our previous bestfitting estimate of 19.6-20.3 mm yr −1 , which might be due to the broad uplift signal north of Krafla central volcano (Fig. 8) that also influences the horizontal displacement rates. All other parameters express only a slight change within or close to the uncertainty limits of our last study.
However, thanks to the improved data set, the estimated model parameter uncertainties have become so small that model uncertainties, which are not assessed in this study, would probably outweigh the propagated data uncertainties. In other words, if we use a different geometry for the plate boundary segments or another Earth model than an elastic half-space, the estimated model parameters would likely change beyond the current model parameter confidence bounds. Therefore, the estimated model parameter uncertainties should be regarded as minimum uncertainties as they only include the effect of the data errors. Further important model assumptions are: A complete stress release in 1872 when the last large earthquakes hit the HFF, a steady stress accumulation since then and a constant locking depth along the HFF segments. This last simplifying assumption is needed as the locking depth is mostly constrained by data points close to the southeastern end of the HFF, leaving no control over the northwestern part of the fault. Finally, we do not know exactly the effect of the Krafla rifting episode on the stress-regime of the HFF. Coulomb failure stress calculations derived from an opening dyke model at Krafla suggest a stress drop on the eastern part of the HFF, which might have relieved some of the accumulated stress on the fault .
In comparison to our previous study, the implications on the seismic potential of the HFF fault do not change much and are in the same range of the 1755 and 1872 earthquakes. The slightly faster slip rate is neutralized by a slightly lower locking depth. To put the obtained modelling results in a general context, Table 4 compares the results of previous studies that provide information about the locking depth and slip rate of the HFF. Some of these studies are based in part on the same CGPS and EGPS data we use in this paper, while others are based on geological information (Saemundsson 1974), seismicity analysis (Rögnvaldsson et al. 1998) or ISNET campaign GPS data (Árnadóttir et al. 2009). All these studies generally agree on the slip rate and locking depth of the HFF. A notable exception is the locking depth inferred from the results of Rögnvaldsson et al. (1998) that is almost twice as large as the estimates based on geodetic data (Metzger et al. 2011, and this study). The estimate from Rögnvaldsson et al. (1998) is deduced from relocated earthquakes with most of the epicentres close to the northwestern end of the HFF while the geodetically derived depths were constrained primarily by data points located near the southeastern end of the fault. The temperature gradient in Iceland can reach 150-200 • C km −1 at the flanks of the rift zone and decreases to 40-50 • C km −1 in the oldest crust in East and West Iceland (Palmason & Saemundsson 1979). This means that the eastern end of the HFF that links to the NVZ likely has a higher temperature gradient than the Western, offshore end of the fault. Furthermore the seismic activity in Icelandic crust ceases at temperatures above 600-800 • C and the crust becomes partially molten at 1200 • C (Flóvenz & Saemundsson 1993;Björnsson 2006). A map showing the estimated depth to the 1200 • C-isotherm in Iceland published in Flóvenz & Saemundsson (1993) indicates increasing isotherm depth along the HFF from less than 15 km close to the triple junction at Theistareykir central volcano to more than 20 km at the northwestern end of the fault. The same applies for the heat flow that decreases from 140 to 175 mW m −2 , measured in boreholes close to the eastern end of the HFF in the Theistareykir area, to 80 mW m −2 on Flatey island. Assuming linear temperature gradient with depth, this would result in a variable thickness of the seismogenic zone and locking depth along the fault.
To calculate the seismic potential of the HFF we combine the fault slip-rate, length and locking depth of the HFF, which results in a tightly constrained moment that corresponds to a M w 6.81 ± 0.04 earthquake. The model segment of the HFF is connected to the Krafla rift segment but in reality the fault is ∼20 km shorter and ends in the Theistareykir fissure swarm (Fig. 1). As the seismic potential scales with the size of the fault plane and the accumulated slip (Aki 1966;Hanks & Kanamori 1979), the shorter fault length reduces the seismic potential by 20 per cent to M w 6.75. On the other hand, if the locking depth is not uniform along the fault as our model assumes, but increases from 6.2 km in the east to 12 km in the West Rögnvaldsson et al. (1998), the fault loading increases by 15 per cent and the seismic potential becomes M w 6.85. Given this variability, we therefore adjust the estimated seismic potential and its uncertainty to M w 6.8 ± 0.1.

C O N C L U S I O N S
The key objective of the work presented in this study was to derive a kinematic model of the Tjörnes Fracture Zone in North Iceland in order to estimate the locking depth and slip deficit of the HFF, and thus its seismic potential that has accumulated since the last large earthquakes in 1872. In our former paper we used only 14 CGPS stations to constrain parameters of an interseismic back-slip model of the plate boundary in North Iceland (Metzger et al. 2011). In this paper we almost quadruple the number of input GPS velocities by adding EGPS data dating back to 1997 and find similar optimal model parameters while all model parameter uncertainties are reduced by more than 50 per cent. We confirm our uncertainty estimations based on data error propagation with a Bayesian estimation. However, with the increased input data density, the derived model parameter uncertainties have become so small that the ambiguity caused by the choice of the model would probably outweigh the propagated data uncertainties and would likely significantly increase the overall model parameter uncertainties. Compared to our earlier study we find a slightly larger slip-rate and a slightly shallower locking depth of the HFF, resulting in an unchanged accumulated seismic moment, corresponding to a M w 6.8 ± 0.1 earthquake.
All model parameters, except for the ridge locking depth, changed within the estimated uncertainties. Relocated off-shore earthquakes along the western half of the HFF (Rögnvaldsson et al. 1998) suggest deeper locking (10-12 km) than what we obtain from on-land geodetic data near the eastern end of the fault (6.2 km). This possible along-strike variation in locking depth is supported by significant change in thermal gradient along the fault from a high gradient of 80 • C km −1 near the NVZ in the east to probably ∼50 • C km −1 in the west. If the locking depth increases gradually from 6.2 km in the east to 12 km in the west, it would mean that the accumulated stress is ∼15 per cent higher than estimated above, or equivalent to a M w 6.85 earthquake.

A C K N O W L E D G M E N T S
The collection of EGPS data up to 2007 was conducted by Thierry Villemin and François Jouanne, with financial support of the French