## Abstract

Rifting occurs as episodes of active deformation in individual rift segments of the Northern Volcanic Zone (NVZ) in Iceland. Here, we simulate deformation around the Krafla central volcano and rift system in the NVZ using a 3-D numerical model in order to explain synthetic aperture radar data acquired by the ERS and Envisat satellite missions between 1993 and 2008. The deformation is non-linear in time over the observed interval. The observed deformation can be explained by a combination of three processes, including: (i) secular plate spreading between the North American and Eurasian plates at a rate of 18.2 mm yr−1, (ii) viscoelastic relaxation following the Krafla Fires rifting episode between 1975 and 1984 and (iii) inflation/deflation of shallow magma chambers beneath the Theistareykir and Krafla central volcanoes. We minimize the misfit between the observed and modelled values of the range change gradient, averaged over all samples, using a simulated annealing algorithm that uses a first-order Taylor series to approximate the fitting function. The calibration parameters include the locking depth of the plate boundary and the rheological properties of the lower crust and mantle. The 68-per cent confidence intervals for the parameters in the solution that best fits the data are: (i) a locking depth of 8.0 to 9.5 km, (ii) a viscosity of 19 to 49 EPa.s (1 EPa.s =1018 Pa.s) in the lower crust at depths between 8 and 24 km and (iii) a viscosity of 5 to 9 EPa.s in the upper mantle below 24 km.

## INTRODUCTION

The Northern Volcanic Zone (NVZ) in Iceland accommodates 18.2 ± 0.4 mm yr−1 of spreading at the divergent boundary separating the North American and Eurasian plates (DeMets et al.2010). It extends from the Tjornes Fracture Zone (TFZ) in the north to the Vatnajokull ice cap in central Iceland, and contains several north-northeast trending volcanic rift systems each composed of a central volcano and an associated fissure swarm, as shown in Fig. 1. Episodic rifting at these volcanoes, fed by deep sources, leads to the creation of new crust (Buck et al.2006). In the last 300 yr, three major rifting episodes have occurred within the NVZ, two at Krafla and one at Askja (Einarsson 2008). The most recent rifting episode, known as the Krafla Fires, occurred at Krafla between 1975 and 1984, after ∼250 yr of quiescence. During the rifting episode, there were about 20 distinct diking events, of which nine were eruptive (Sigmundsson 2006). The diking occurred along a 80-km-long rift segment, resulting in an average total opening of 5 m (Tryggvason 1984) that corresponds to ∼250 yr of interrifting deformation.

Figure 1.

Map showing the location of the Northern Volcanic Zone (NVZ) in Iceland. The ∼80-km long segment that accommodated rifting during the 1975–1984 rifting episode at Krafla is shown by the thick red line. Outlines of central volcanoes and their associated fissure swarm are shown by solid black lines (Johannesson 2009). Dashed line shows the study area. Solid red dot shows the location of the caldera at Krafla. Coordinates are in longitude and latitude. Inset shows boundary (red line) between the North American and European plates and their relative plate motion as calculated from the MORVEL plate model (DeMets et al.2010).

Figure 1.

Map showing the location of the Northern Volcanic Zone (NVZ) in Iceland. The ∼80-km long segment that accommodated rifting during the 1975–1984 rifting episode at Krafla is shown by the thick red line. Outlines of central volcanoes and their associated fissure swarm are shown by solid black lines (Johannesson 2009). Dashed line shows the study area. Solid red dot shows the location of the caldera at Krafla. Coordinates are in longitude and latitude. Inset shows boundary (red line) between the North American and European plates and their relative plate motion as calculated from the MORVEL plate model (DeMets et al.2010).

Like earthquakes, rifting events transfer stress to the surrounding region. Since the weak layers of the lithosphere and asthenosphere, such as the lower crust and upper mantle, cannot support high shear stresses, they relax over time, deforming the elastic upper crust. This transient deformation can last several decades, especially after large earthquakes (e.g. Hu et al.2004; Ali & Freed 2010), and has been used to study the rheology of lithosphere both at subduction zones and transform boundaries (Bürgmann & Dresen 2008). In contrast, spreading centres have been the focus of relatively few studies because most of the rifting occurs under water, where it is difficult to measure deformation. Modern satellite geodetic techniques have measured the deformation associated with large subaerial rifting events at only three locations: Krafla in Iceland (1975–1984), and Asal-Ghoubbet (1978) and Dabbahu (2005–2010) in Afar, as recently reviewed by Wright et al. (2012). Of these, Krafla has the longest time-series of observations recording post-rifting deformation, making it an excellent natural laboratory for a rheological study.

Using rift normal displacement observations from GPS surveys conducted in 1987 and 1990, following the Krafla Fires, Foulger et al. (1992) and Heki et al. (1993) estimated viscosity beneath the NVZ using 1-D models that assume a thin elastic layer over a viscous layer. Hofton & Foulger (1996) and Pollitz & Sacks (1996) used additional data from a GPS survey in 1992 to estimate viscosities using semi-analytical models. To date, the study of Pollitz & Sacks (1996) is the only one to account for steady plate spreading. In these studies, a number of candidate rheologies were tested, but no formal inverse modelling was performed to estimate the viscosity. More recently, Zeeuw-van Dalfsen et al. (2004) analysed interferometric synthetic aperture radar (InSAR) observations made between 1993 and 1998 over the NVZ. Using inverse modelling in combination with elastic half-space models, they suggested that accumulation of magma at a depth of 20 km, rather than viscoelastic relaxation, could explain the widespread uplift. As we show below, InSAR is well suited for studying viscoelastic relaxation at rift zones because it captures the dominant vertical component of the displacement field near the rift. In this paper, we combine InSAR measurements of deformation in the NVZ with viscoelastic models to study post-rifting relaxation following the Krafla Fires. Unlike Zeeuw-van Dalfsen et al. (2004), we test the hypothesis that viscoelastic relaxation following Krafla Fires can explain the widespread uplift observed using InSAR.

We describe the InSAR measurements made between 1993 and 2008 around Krafla using a model that accounts for all major tectonic and volcanic processes driving deformation during the post-rifting time interval, that is, post-rifting viscoelastic relaxation following the Krafla Fires, steady plate spreading and inflation/deflation of shallow magma chambers beneath the two central volcanoes. To do so, we use a 3-D numerical model along with a new strategy for inverse modelling that uses the gradient of range change as the observable quantity (Ali & Feigl 2012). In order to isolate the relative contribution of the underlying processes and to estimate geophysical parameters along with their formal uncertainties, that best explain the observations, we solve the non-linear inverse problem.

## DATA

We analyse SAR images acquired by the ERS (Duchossois & Guignard 1987) and Envisat (McLeod et al.1998) satellite missions of the European Space Agency on 14 distinct epochs between 1993 and 2008. Fig. 1 shows the geographic location of the subset of the ERS scene denoted by frame 2277 of Track 9. We combine the 14 distinct epochs to form nine ERS interferometric pairs (Carr 2008) and two Envisat pairs, as listed in Table 1. These pairs form three independent sets, as shown in Fig. 2. The 11 pairs constitute a minimal set spanning the observed time interval, as shown by the lack of closed loops in the incidence graph. Our data set includes most of the ‘good’ pairs analysed by DiCaprio (2010).

Figure 2.

Orbital separation versus time for the interferograms. Horizontal axis displays the acquisition date (epoch) of each image, labels next to circles denote orbit numbers, and the vertical axis shows the orbital separation at the acquisition epoch. Solid lines connect epochs used to form interferometric pairs in independent sets (also known as ‘species’) A (red), B (blue) and C (green).

Figure 2.

Orbital separation versus time for the interferograms. Horizontal axis displays the acquisition date (epoch) of each image, labels next to circles denote orbit numbers, and the vertical axis shows the orbital separation at the acquisition epoch. Solid lines connect epochs used to form interferometric pairs in independent sets (also known as ‘species’) A (red), B (blue) and C (green).

Table 1.

Interferometric pairs analysed in this study showing SAR images acquired by the ERS and Envisat (ENV) satellites. The term Ha denotes the altitude of ambiguity (Massonnet & Rabaute 1993).

Pair Sat. Ha Species First epoch (tiSecond epoch (tjSpan
(m)  Orbit Year Orbit Year (days)
ERS −2628.2 10 174 1993.4822 17398 1998.6274 1879
ERS −53.8 10 174 1993.4822 5875 1996.4235 1074
ERS −63.1 10 675 1993.5781 11386 1997.4767 1424
ERS 64.9 10 675 1993.5781 23410 1999.7781 2264
ERS −198.6 11 677 1993.7699 22408 1999.5863 2124
ERS 70.3 11 677 1993.7699 6877 1996.6147 1039
ERS 190.2 5875 1996.4235 23410 1999.7781 1225
ERS 392.6 6376 1996.5191 22408 1999.5863 1120
ERS 82.2 6877 1996.6147 17899 1998.7233 770
10 ENV 268.8 18 099 2005.6219 28620 2007.6356 735
11 ENV −370.6 22 608 2006.4849 28620 2007.6356 420
Pair Sat. Ha Species First epoch (tiSecond epoch (tjSpan
(m)  Orbit Year Orbit Year (days)
ERS −2628.2 10 174 1993.4822 17398 1998.6274 1879
ERS −53.8 10 174 1993.4822 5875 1996.4235 1074
ERS −63.1 10 675 1993.5781 11386 1997.4767 1424
ERS 64.9 10 675 1993.5781 23410 1999.7781 2264
ERS −198.6 11 677 1993.7699 22408 1999.5863 2124
ERS 70.3 11 677 1993.7699 6877 1996.6147 1039
ERS 190.2 5875 1996.4235 23410 1999.7781 1225
ERS 392.6 6376 1996.5191 22408 1999.5863 1120
ERS 82.2 6877 1996.6147 17899 1998.7233 770
10 ENV 268.8 18 099 2005.6219 28620 2007.6356 735
11 ENV −370.6 22 608 2006.4849 28620 2007.6356 420

We generate the interferograms using the DIAPASON software developed by the French Space Agency, CNES, as described previously (Massonnet & Rabaute 1993; Massonnet & Feigl 1998). The topographic contribution to the interferograms is removed using a digital elevation model that has been resampled to 100 m posting and 20 m accuracy (Arnason 2006). The coordinate system in easting X and northing Y is called ISN93, as defined by Rennen (2002). The observed wrapped phase change values for the ERS pair spanning time interval 1993.48–1998.62 is shown in Fig. 3(a). The observed wrapped phase change values for pairs 2–11 are shown in the first column of Figs 4 and 5. One fringe of phase change in these interferograms corresponds to 28 mm of range change along the radar line of sight between the satellite and the ground (shown by the black arrow). The interferograms exhibit multiple signatures with different characteristic scales in space and time. For example, a decrease in range in the zone roughly 50 km in diameter north of the central volcano at Krafla in Fig. 3(a) has been attributed to inflation at a depth of 20 km, whereas small concentric fringes at the central volcano result from deflation of the shallow magma chamber (Zeeuw-van Dalfsen et al.2004). Also, contributing to the range change field is the deformation due to steady plate spreading.

Figure 3.

(a) Interferogram showing observed values of wrapped phase change Δϕ for pair 1, spanning the time interval 1993.48–1998.63. One coloured fringe corresponds to one cycle of phase change, or 28 mm of range change. Black arrow shows the projection of the look vector (from sensor to target) onto the horizontal surface. (b) Observed values of the range gradient ψ in 28-mm cycles per 100-m pixel such that 0.05 cycles pixel−1 corresponds to 1.4×10−4 or 140 microstrain after resampling by the quad-tree algorithm. (c) Modelled values of the range change gradient ψ calculated from the final estimate of the parameters in the model for this individual pair. (d) Residual range gradient formed by subtracting the modelled values from the observed values. Gray areas have been excluded by the quadtree resampling procedure. All coordinates are in easting and northing in kilometers in the ISN93 Lambert projection (Rennen 2002).

Figure 3.

(a) Interferogram showing observed values of wrapped phase change Δϕ for pair 1, spanning the time interval 1993.48–1998.63. One coloured fringe corresponds to one cycle of phase change, or 28 mm of range change. Black arrow shows the projection of the look vector (from sensor to target) onto the horizontal surface. (b) Observed values of the range gradient ψ in 28-mm cycles per 100-m pixel such that 0.05 cycles pixel−1 corresponds to 1.4×10−4 or 140 microstrain after resampling by the quad-tree algorithm. (c) Modelled values of the range change gradient ψ calculated from the final estimate of the parameters in the model for this individual pair. (d) Residual range gradient formed by subtracting the modelled values from the observed values. Gray areas have been excluded by the quadtree resampling procedure. All coordinates are in easting and northing in kilometers in the ISN93 Lambert projection (Rennen 2002).

Figure 4.

Interferograms showing (from left to right in each row) observed phase as well as the observed, modelled and residual values of the gradient of range change for pairs 2–6. Plotting conventions as in Fig. 3.

Figure 4.

Interferograms showing (from left to right in each row) observed phase as well as the observed, modelled and residual values of the gradient of range change for pairs 2–6. Plotting conventions as in Fig. 3.

Figure 5.

Interferograms showing (from left to right in each row) observed phase as well as the observed, modelled and residual values of the gradient of range change for pairs 7–11. Plotting conventions as in Fig. 3.

Figure 5.

Interferograms showing (from left to right in each row) observed phase as well as the observed, modelled and residual values of the gradient of range change for pairs 7–11. Plotting conventions as in Fig. 3.

Following Ali & Feigl (2012), we calculate the dimensionless gradient of range change ψ for the kth pixel as:

(1)
$$\psi ^{(k)}=\frac{\Delta \rho ^{(k+1)}-\Delta \rho ^{(k-1)}}{X^{(k+1)}-X^{(k-1)}},$$
where Δρ is the range change, and X, the easting coordinate, both of which have dimensions of distance. We choose the east component of the range change gradient because it is sensitive to the details of the deformation. Accordingly, we neglect the north component. Just as the range change measured by InSAR is sensitive to displacement, the range change gradient is sensitive to strain. Mathematically, the range change gradient is equivalent to a particular component of the deformation gradient tensor (Malvern 1969). Since the range change gradient is a continuous and differentiable quantity (Sandwell & Price 1998), using it as the observable avoids the pitfalls associated with phase unwrapping techniques.

The range change gradient values are derived from wrapped phase data by a quad-tree resampling procedure, as summarized by Ali & Feigl (2012). It reduces the computational burden and mitigates correlations between neighbouring pixels. For example, the quad-tree resampling procedure reads the phase values shown in Fig. 3(a) to extract the range change gradient shown in Fig. 3(b). The original InSAR data set is reduced to 1500–3250 values of the range change gradient for each of the interferometric pairs.

## MODELLING STRATEGY

To describe the signatures observed in the interferograms, we account for three distinct processes, including: (i) post-rifting relaxation following the Krafla Fires, (ii) steady plate spreading and (iii) inflation/deflation of the magma chambers beneath the Krafla and Theistareykir central volcanoes, as sketched in Fig. 6. We model these three processes separately and sum the solutions to calculate the total deformation field.

Figure 6.

Cartoon showing active sources of deformation at the plate boundary during the time interval between 1993 and 2008. Post-rifting relaxation is modelled by imposing opening (during 1975–1981) on a 80-km rift segment in 3-D viscoelastic media. The viscosities of the lower crust, transitional layer and mantle are ηL, ηT and ηM, respectively. Interrifting deformation rate or steady plate spreading is modelled by imposing an annual slip of 9.1 mm on each of the two horizontal detachment surfaces at a locking depth Dl, a free parameter estimated in the inversion. Volcanic deformation is modelled using a spherical source in an elastic half-space at depth D, corresponding to DK or DT for Krafla or Theistareykir, respectively. The cumulative deformation field is the sum of deformation fields that result from the three sources.

Figure 6.

Cartoon showing active sources of deformation at the plate boundary during the time interval between 1993 and 2008. Post-rifting relaxation is modelled by imposing opening (during 1975–1981) on a 80-km rift segment in 3-D viscoelastic media. The viscosities of the lower crust, transitional layer and mantle are ηL, ηT and ηM, respectively. Interrifting deformation rate or steady plate spreading is modelled by imposing an annual slip of 9.1 mm on each of the two horizontal detachment surfaces at a locking depth Dl, a free parameter estimated in the inversion. Volcanic deformation is modelled using a spherical source in an elastic half-space at depth D, corresponding to DK or DT for Krafla or Theistareykir, respectively. The cumulative deformation field is the sum of deformation fields that result from the three sources.

### Post-rifting relaxation

We model post-rifting viscoelastic relaxation following the Krafla Fires using Defmod (Ali 2011a), an open-source, fully unstructured, implicit/explicit, 2/3-D, parallel finite element code. Defmod is especially well suited for modelling deformation due to earthquakes and volcanoes over timescales ranging from milliseconds to thousands of years. Fault slip and/or tensile opening is accommodated via additional linear constraint equations, implemented using Lagrange multipliers. The quasi-static viscoelastic problem is solved using sparse implicit solvers available in the portable, extensible toolkit for scientific computation (Balay et al.2012). The time-stepping algorithm for viscoelastic relaxation is similar to that used by Melosh & Raefsky (1980).

The hexahedral mesh used in our model spans 576 km in length, 576 km in width and 288 km in depth, containing ∼750 thousand degrees of freedom. The domain is large enough to avoid edge effects (Dubois et al.2008). The length of a single element's edge varies from 72 km at the bottom boundary to 0.1 km at the tip of the rift near the brittle-ductile transition (as shown in Supporting Information Fig. S1). The small elements are necessary for accurate modelling of viscoelastic relaxation. We assume a brittle elastic crust, 8 km in thickness, based on the depth of regional seismicity (Wright et al.2012), a 16-km-thick viscoelastic lower crust, a 8-km-thick viscoelastic layer that forms the transition between lower crust and mantle and finally a viscoelastic mantle below 32 km. The 24–32 km depth of the Mohorovičić discontinuity in our model is consistent with the seismic studies of Brandsdóttir et al. (1997) and Darbyshire et al. (2000). We assume a constant Poisson's ratio of 0.25 for all layers and a Young's modulus of 75 GPa for the elastic upper crust, 80 GPa for the lower crust, 85 GPa for the transition zone and 90 GPa for the mantle. These values are consistent with a shear modulus of 30–36 GPa. The viscosities of the three ductile layers (i.e. ηL, ηT and ηM) are estimated during the inversion. To allow large variations in viscosity, the parametrization uses a (base-10) logarithmic scale such that η* = Log10(η/1 Pa.s).

The strike of the vertical rift is fixed in the model at N6°E to follow the distribution of seismicity, observed during rifting (Wright et al.2012). This value is slightly more northerly than the value of N11°E estimated from surface observations (Tryggvason 1984). All boundary nodes, except those on the top free surface (i.e. maximum positive value of Z coordinate) are fixed. The rifting source model consists of a 80-km-long by 8-km-wide segment, discretized using 5115 coincident nodes that open steadily during a 6-yr interval between 1976 and 1981. The amount of opening is based on the observations of Tryggvason (1984) and tapers smoothly from 7.2 m near the northern tip of caldera, to 0 m at the three buried edges, that is the north, south and bottom boundaries of the rift segment. Although we could make the amount of tensile opening a free parameter, doing so causes a trade-off with the viscosity of the ductile layers.

### Interrifting deformation

To describe the interrifting deformation due to motion between the North American and European plates at a constant rate of 18.2 mm yr−1 (DeMets et al.2010), we use a kinematic model. As sketched along a 2-D cross-section in Fig. 6, the situation is analogous to two plates sliding away from each other on two horizontal detachment surfaces separated by the bottom edge of the rift zone. Assuming a half-space with uniform elastic properties, the interrifting model specifies horizontal displacements as two rectangular dislocations on a horizontal plane at depth Dl, as formulated by Okada (1985). The annual displacement on each dislocation plane is 9.1 mm. The azimuth of the dislocation vector is 276 degrees on the North American side of the plate boundary and 096 degrees on the European side. The red curve in Fig. 7 plots the resulting horizontal interrifting velocities along a profile crossing the rift. The sigmoidal shape is analogous to that produced by interseismic deformation across a strike-slip fault (Feigl & Thatcher 2006). Its half-width scales with the locking depth Dl and the height of the ‘step’ scales with the long-term plate spreading rate.

Figure 7.

Red curve shows horizontal velocities, in the rift-perpendicular direction, across the rift, due to interrifting deformation using the approach described in Fig. 6 with a locking depth Dl = 8.75 km. The blue curve shows the velocity that would result due to steady opening of a dike that extends from the surface to a depth of 8.75 km. The black curve shows velocity that would result by adding the red and the blue curves.

Figure 7.

Red curve shows horizontal velocities, in the rift-perpendicular direction, across the rift, due to interrifting deformation using the approach described in Fig. 6 with a locking depth Dl = 8.75 km. The blue curve shows the velocity that would result due to steady opening of a dike that extends from the surface to a depth of 8.75 km. The black curve shows velocity that would result by adding the red and the blue curves.

At the free surface, the displacement field calculated using this approach is equivalent to that calculated by a conventional model. By imposing the tensile opening u as a dislocation on a vertical rift zone extending from the free surface to depth Dl, the conventional approach can calculate the co-rifting displacement field from a single diking episode. One such episode occurs, on average, with a recurrence interval of T years. Then, the conventional approach calculates the interrifting velocity field as the difference between a rigid-body translational displacement field ±uT/2 and the co-rifting displacement field. After normalizing by the recurrence interval T, the conventional approach yields the interrifting velocity field. It exactly agrees with the velocity field calculated using our kinematic approach. We prefer our approach because it is less sensitive to inaccuracies in the trace of the rift separating the two plates.

Using a kinematic model to calculate interrifting deformation offers a considerable advantage in computational speed over a dynamic model driven by far field boundary conditions. When applied in subduction zones, dynamic models using the finite element method must evaluate a long time interval, including many large earthquakes at a fixed recurrence interval, before reaching a steady state. Known as ‘spinning up’ (Hetland & Hager 2006), this procedure is expensive in terms of computation time. Ali & Freed (2010) have shown that the kinematic and dynamic approaches yield equivalent results in a subduction zone, provided that the appropriate configuration is chosen. In particular, the dynamic formulation must use an appropriate, possibly non-linear, rheology (Ali & Freed 2010).

### Inflation and deflation of magma chamber beneath volcanoes

The InSAR phase data (Fig. 3a) show rapid deformation, concentrated within a circular area roughly ∼5 km wide, located inside the central caldera at Krafla (outlined in black in Fig. 3). Following previous studies (Zeeuw-van Dalfsen et al.2004; Ali & Feigl 2012), we interpret this signal as deflation of a shallow magma chamber beneath the caldera. Similarly, the intermittent deformation observed at Theistareykir volcano, northwest of Krafla, has been attributed to inflation of a shallow magma chamber (Metzger et al.2011). To model these short-wavelength signals, we assume that a spherical source changes volume at a rate of $$\Delta \dot{V}$$ cubic meters per year, a parameter that we estimate in the inverse modelling described below. The depth of the source is described by another free parameter D. For Krafla, the horizontal position of the source is fixed at ISN93 coordinates (X, Y) = (601.605, 580.917) km, as estimated separately from InSAR data (Ali & Feigl 2012), while the depth DK is estimated during the inversion. For the magma chamber beneath Theistareykir, the source coordinates, including depth DT, are estimated during the inversion.

To calculate the deformation caused by the inflation or deflation of each magma chamber, we apply the formulation attributed to Mogi (1958) as expressed by Segall (2010). It assumes uniform material properties everywhere in an elastic half space. Although this assumption can slightly bias estimates of the depth D (Masterlark 2007), the half-space approximation is justified because (i) the modelled depths of the two sources (DKDT ∼ 5 km) are much greater than the topographic relief of the volcanic edifice (Δh ∼ 300 m), reducing the sensitivity of the deformation to the structural details of the system (Cayol & Cornet 1998), and (ii) the ∼5-km length scales of the inflation/deflation signatures at Krafla and Theistareykir are at least an order of magnitude shorter than that of the post-rifting signature.

## OPTIMIZATION

To account for the cumulative deformation due to post-rifting viscoelastic relaxation, steady interrifting spreading and inflation/deflation of the magma chamber beneath the central volcanoes, we sum the displacement fields from the respective calculations. Thus configured, our model defines the fitting function that calculates the range change gradient from a set of parameters. Now, the goal is to estimate the set of parameters that best explain the InSAR measurements. We do so by using the strategy developed and validated by Ali & Feigl (2012), as summarized below.

To describe the observed range change gradient field ψ, we seek a modelled field $$\tilde{\psi }$$ defined in terms of a set of m parameters $$\vec{p}$$. We evaluate the fitting function $$\tilde{\psi }(\vec{p})$$ at each resampled pixel k over each time interval [ti, tj]. This procedure is also called solving the forward problem. To quantify the misfit between the observed ψ and modelled $$\tilde{\psi }$$ values of the range gradient, the objective function calculates the cost $$\bar{\omega }^{\prime }$$ in terms of their difference, averaged over all n points in the resampled data set, that is

(2)
$$\bar{\omega }^{\prime }=\frac{1}{n}\sum _{k=1}^{n} {arc}(\psi ^{(k)},\tilde{\psi }^{(k)}),$$
where the arc function, defined by Feigl & Thurber (2009, eqs (14)–(16) and references therein), returns the lesser (in absolute value) of the two angles separating the direction of the observed range gradient from that of the modelled range gradient. The objective function is minimized using a parallelized version (Ali 2011b) of the simulated annealing algorithm (Kirkpatrick et al.1983) implemented by Goffe (1996). This algorithm performs several thousand evaluations of the fitting function to find the optimal solution: the set of model parameters that produces the lowest value of cost. We constrain the viscosity to decrease with depth, such that ηL ≥ ηT ≥ ηM. To do so, we modify the objective function to include a penalty function such that:
(3)
$$\omega ^{\prime \prime }=\bar{\omega }^{\prime }+c\left[{\rm max}(\eta ^{*}_{M}-\eta ^{*}_{T},0)^{2}+{\rm max}(\eta ^{*}_{T}-\eta ^{*}_{L},0)^{2})\right],$$
where c = 100 is a constant.

For computational efficiency, the fitting function is approximated using a first-order Taylor series. To evaluate its first-order terms, we calculate the partial derivative of the range change gradient $$\delta \tilde{\psi }/\delta p$$ at each pixel k with respect to each parameter p using a forward-finite difference approximation. Although these calculations require (m + 1) evaluations of the exact (and slow) fitting function, they may be performed independently in parallel.

The simulated annealing algorithm then evaluates the approximate (and fast) version of the fitting function. As the algorithm performs these two steps several times, the parameter values converge on the final estimate, typically after a few iterations. In practice, the scheme is comparable to Newton's method for finding the (local) minimum of a function in the neighbourhood where the first-order Taylor approximation is valid. As discussed by Ali & Feigl (2012), the inversion scheme performs well when the number of parameters m is small. For large numbers of parameters in the viscoelastic finite-element model, adjoint state methods would be more suitable (e.g. Al-Attar & Tromp 2013).

In our case, the free parameters include: the locking depth Dl of the plate boundary, the viscosity of the lower crust ηC, the viscosity of the transition zone ηT, the viscosity of the mantle ηM, the depth of the deflating magma chamber at Krafla DK, its rate of volume change $$\Delta \dot{V}_{K}$$ and depth DT, as well as the latitude and longitude (X, Y)T of the magma chamber beneath Theistareykir, along with its rate of volume change $$\Delta \dot{V}_{T}$$.

To calculate the uncertainty of the parameters, we apply bootstrap resampling (Efron & Tibshirani 1986). We create 100 realizations of the data set by resampling, with replacement, the n measurements. For each realization, we minimize the objective function using the simulated annealing algorithm and the approximated fitting function from the final iteration. We then calculate the critical value of the cost as the sum of the optimum (minimum) cost calculated from the final iteration plus the sample standard deviation of the cost obtained from bootstrapping. The critical value of the cost sets a ‘height’ in the valley-like curve describing the cost as a function of a parameter (Fig. 8 ). For cost values below this critical ‘height’, the corresponding values of the parameter are statistically equivalent to the optimal value. For each parameter, the 68-per cent confidence interval is delimited by the minimum and maximum of the trial values for which the cost is less than the critical value (and necessarily greater than the minimum, optimal value) of the cost. The trial values of the parameters are those generated by the simulated annealing algorithm in the final iteration using the approximate fitting function and the original data-set (without bootstrap resampling).

Figure 8.

Cost as a function of parameter value: (a) lower crustal viscosity, (b) transition zone viscosity, (c) mantle viscosity and (d) locking depth. The critical value of cost, at 68-per cent confidence, is 0.03017 cycle pixel−1.

Figure 8.

Cost as a function of parameter value: (a) lower crustal viscosity, (b) transition zone viscosity, (c) mantle viscosity and (d) locking depth. The critical value of cost, at 68-per cent confidence, is 0.03017 cycle pixel−1.

## RESULTS

We now apply the strategy to estimate the parameters in the model, as described in Section 3. The initial estimate for each free parameter is listed in Table 2. The initial estimates of viscosity are of the same order of magnitude as those of Hofton & Foulger (1996), that is ∼1 EPa.s (or 1018 Pa.s) for the lower crust, transition zone, and mantle. The lower and upper bounds for the parameters and the step sizes used for calculating the partial derivatives are listed in Table 2. To account for errors in the orbital trajectories, we estimate an additional unknown, called the ‘nuisance’ parameter, for all but one of the epochs in a species. The result appears as a ramp in the phase field or equivalently, as an additive constant in the range change gradient. We perform the inversion using the resampled data for each pair individually. We also combine the resampled data from all 11 pairs into an ‘ensemble’ solution. To account for time-dependent deformation at the Krafla central volcano, in the ensemble solution, we use a piecewise linear approximation to the exponential temporal function given by Ali & Feigl (2012). For Theistareykir central volcano, we assume a linear time function (Metzger 2012).

Table 2.

Initial estimate, bounds and step size used for the optimization. $$\eta ^{*}_{L}$$, $$\eta ^{*}_{T}$$ and $$\eta ^{*}_{M}$$ (on a logarithmic scale) denote viscosities of lower crust (8–24 km), transitional zone (24–32 km) and mantle (below 32 km). Dl denotes the locking depth. DK and DT denote the depth of shallow magma chambers beneath the Krafla and Theistareykir central volcanoes, respectively. $$\Delta \dot{V}$$ denotes the annual rate of volume change.

Parameter name $$\eta ^{*}_{L}$$ $$\eta ^{*}_{T}$$ $$\eta ^{*}_{M}$$ Dl DK $$\Delta \dot{V}_{K}$$ DT $$\Delta \dot{V}_{T}$$
units Log10(η/1 Pa.s) (km) (km) ( × 105m3 yr−1(km) ( × 105 m3 yr−1
Lower bound 18.00 18.00 18.00 7.5 4.50 −50.0 1.00 0.0
Upper bound 21.00 21.00 21.00 8.5 5.50 0.0 7.50 15.0
Step size 0.25 0.25 0.25 0.5 0.25 1.0 0.25 1.0
Initial Est. 18.50 18.25 18.00 8.0 4.50 0.0 1.00 0.0
Parameter name $$\eta ^{*}_{L}$$ $$\eta ^{*}_{T}$$ $$\eta ^{*}_{M}$$ Dl DK $$\Delta \dot{V}_{K}$$ DT $$\Delta \dot{V}_{T}$$
units Log10(η/1 Pa.s) (km) (km) ( × 105m3 yr−1(km) ( × 105 m3 yr−1
Lower bound 18.00 18.00 18.00 7.5 4.50 −50.0 1.00 0.0
Upper bound 21.00 21.00 21.00 8.5 5.50 0.0 7.50 15.0
Step size 0.25 0.25 0.25 0.5 0.25 1.0 0.25 1.0
Initial Est. 18.50 18.25 18.00 8.0 4.50 0.0 1.00 0.0

The results for pair 1 are shown in Fig. 3 and the values of estimated parameters are listed in Table 3. Fig. 3(b) shows the observed values of the range change gradient ψ after quadtree resampling. The modelled values of the range change gradient $$\tilde{\psi }$$ calculated using the final estimate are shown in Fig. 3(c). The inversion takes five iterations before each parameter converges to a steady value, during which the cost decreases from ω = 0.0427 cycles per pixel to ω = 0.0273 cycles per pixel. The main features, such as the two large elliptical lobes on either side of the rift axis, are reproduced by the model. The modelled values of the range gradient from the three physical processes that contribute to the total deformation field over the period of observation are shown individually in Fig. 9. The deformation resulting from viscoelastic relaxation is greater than that caused by plate spreading during the 5-yr post-rifting interval spanned by the interferogram. The contribution from the shallow deflating magma chamber at Krafla is small and focused near the central volcano. Similarly, the model accounts for the deformation at Theistareykir.

Figure 9.

Maps for pair 1, spanning the time interval 1993.48–1998.63, showing (a) modelled values of the range change gradient calculated from the final estimate of the parameters in the model; (b) contribution of viscoelastic relaxation to the modelled values calculated using the finite element code Defmod (Ali 2011a); (c) contribution of plate spreading to the modelled values calculated using an elastic dislocation model; (d) contribution of subsidence and uplift at the Krafla and Theistareykir central volcanoes to the modelled values calculated using an elastic model. All coordinates are in easting and northing in kilometers in the ISN93 Lambert projection (Rennen 2002).

Figure 9.

Maps for pair 1, spanning the time interval 1993.48–1998.63, showing (a) modelled values of the range change gradient calculated from the final estimate of the parameters in the model; (b) contribution of viscoelastic relaxation to the modelled values calculated using the finite element code Defmod (Ali 2011a); (c) contribution of plate spreading to the modelled values calculated using an elastic dislocation model; (d) contribution of subsidence and uplift at the Krafla and Theistareykir central volcanoes to the modelled values calculated using an elastic model. All coordinates are in easting and northing in kilometers in the ISN93 Lambert projection (Rennen 2002).

Table 3.

Model parameters estimated using pair 1 and pairs 1-11 (individually and jointly) along with their bootstrap uncertainties. [1] Over the period 1993–1999.

Parameter name $$\eta ^{*}_{L}$$ $$\eta ^{*}_{T}$$ $$\eta ^{*}_{M}$$ Dl DK $$\Delta \dot{V}_{K}$$ DT $$\Delta \dot{V}_{T}$$
units Log10(η/1 Pa.s) (km) (km) ( × 105m3 yr−1(km) ( × 105 m3 yr−1
Pair 1 only 19.46 18.91 18.91 9.5 4.5 −7.6 7.2 15.0
−1σ Uncertainty 0.15 0.08 0.10 1.2 0.0 5.4 2.1 7.0
+1σ Uncertainty 0.12 0.12 0.08 0.0 1.0 3.8 0.3 0.0
Pairs 1–11 (mean) 19.81 18.96 18.95 9.1 4.8 −7.9 5.1 8.0
−1σ Uncertainty 0.76 0.22 0.33 2.2 0.3 8.0 3.3 10.5
+1σ Uncertainty 0.33 0.56 0.16 0.4 0.7 5.6 1.6 2.5
Pairs 1–11 (ensemble) 19.47 18.81 18.81 9.5 4.5 −6.11 7.5 15.0
−1σ Uncertainty 0.19 0.11 0.13 1.5 0.0 5.2 1.9 10.7
+1σ Uncertainty 0.22 0.28 0.14 0.0 1.0 3.6 0.0 0.0
Parameter name $$\eta ^{*}_{L}$$ $$\eta ^{*}_{T}$$ $$\eta ^{*}_{M}$$ Dl DK $$\Delta \dot{V}_{K}$$ DT $$\Delta \dot{V}_{T}$$
units Log10(η/1 Pa.s) (km) (km) ( × 105m3 yr−1(km) ( × 105 m3 yr−1
Pair 1 only 19.46 18.91 18.91 9.5 4.5 −7.6 7.2 15.0
−1σ Uncertainty 0.15 0.08 0.10 1.2 0.0 5.4 2.1 7.0
+1σ Uncertainty 0.12 0.12 0.08 0.0 1.0 3.8 0.3 0.0
Pairs 1–11 (mean) 19.81 18.96 18.95 9.1 4.8 −7.9 5.1 8.0
−1σ Uncertainty 0.76 0.22 0.33 2.2 0.3 8.0 3.3 10.5
+1σ Uncertainty 0.33 0.56 0.16 0.4 0.7 5.6 1.6 2.5
Pairs 1–11 (ensemble) 19.47 18.81 18.81 9.5 4.5 −6.11 7.5 15.0
−1σ Uncertainty 0.19 0.11 0.13 1.5 0.0 5.2 1.9 10.7
+1σ Uncertainty 0.22 0.28 0.14 0.0 1.0 3.6 0.0 0.0

Fig. 3(d) shows the residual difference $$(\psi -\tilde{\psi })$$ between the observed and modelled values of the range change gradient. Indeed, the modelled values matches the observed values quite well, to within ω = 0.0273 cycles per 100-m pixel (equivalent to 8.3 microstrain) in range gradient. The residual range gradient field in Fig. 3(d) shows some small areas of misfit with spatial scales of 2–5 km. Some of the misfit may be due to tropospheric effects. Some of the mismatch near the southeastern corner of the study area could also be related to localized subsidence due to extension associated with the fissure swarms. Fig. 10 shows (i) the observed wrapped phase change values, (ii) the modelled wrapped phase change values calculated from the final estimate of the parameters and (iii) the residual wrapped phase change values calculated by subtracting the modelled values from the observed values.

Figure 10.

Interferograms showing (a) observed, (b) modelled and (c) residual values of wrapped phase change for pair 1, spanning the time interval 1993.48 − 1998.63. Modelled values are calculated from the final estimate of the parameters. Residual values are calculated by subtracting the modelled values from the observed values. Plotting conventions as in Fig. 3(a).

Figure 10.

Interferograms showing (a) observed, (b) modelled and (c) residual values of wrapped phase change for pair 1, spanning the time interval 1993.48 − 1998.63. Modelled values are calculated from the final estimate of the parameters. Residual values are calculated by subtracting the modelled values from the observed values. Plotting conventions as in Fig. 3(a).

We repeat the procedure for the other pairs (2–11) using the same initial estimate, as listed in Table 2. In all cases, we find that each parameter converges to a steady value after ∼5 iterations. For brevity, we show only the range gradient (Figs 4 and 5, second column) calculated from the observed phase change (Figs 4 and 5, first column) after quadtree resampling, the range gradient calculated from the final estimate (Figs 4 and 5, third column) and the corresponding residuals (Figs 4 and 5, fourth column) calculated by subtracting the modelled values of the range gradient from the observed values. The values of the observed, modelled and residual wrapped phase change are shown in Supporting Information Figs S2 and S3. Fig. 11 shows how the estimated parameters vary across different pairs. The weighted mean value for each of the parameters estimated from pairs 1–11 individually is listed in Table 3. Pairs with large decorrelated regions result in a poor fit, increasing the uncertainty of the estimated parameters. Combining all 11 interferometric pairs in an ensemble solution spanning the post-rifting interval from 1993 to 2008, we estimate a lower crustal viscosity of $$\eta ^{*}_{L}=19.47\pm 0.2$$ between depths of 8 and 24 km, a transition zone viscosity of $$\eta ^{*}_{T}=18.81^{-0.1}_{+0.3}$$ between depths of 24 and 32 km, and a mantle viscosity of $$\eta ^{*}_{M}=18.81\pm 0.14$$ below 32 km. The estimated locking depth Dl at the plate boundary along the NVZ ranges between 8.0 and 9.5 km. This 68-per cent confidence interval is delimited by the upper bound imposed during optimization. The individual parameter estimates agree within their uncertainties (Table 3) to those estimated from other pairs, thus demonstrating that the inversion procedure is robust. The variations in parameters as a function of cost are shown in Fig. 8. When we perform inversion without the penalty function, that is removing the constraint that viscosity of the viscoelastic layers must decrease with depth, we obtain a lower crustal viscosity of $$\eta ^{*}_{L}=19.35{\rm -}19.44$$, a transition zone viscosity of $$\eta ^{*}_{T}=18.37{\rm -}20.67$$, and a mantle viscosity of $$\eta ^{*}_{M}=18.46{\rm -}18.71$$, for the 11-pair ensemble data set. This suggests that the deformation field is not as sensitive to the viscosity of the transition zone as it is to the viscosity of the lower crust or the mantle.

Figure 11.

Variation of parameters with pair number. Solid red line represents the value estimated from the 11 pair ensemble data set and the shaded region represents the corresponding 68 per cent confidence interval.

Figure 11.

Variation of parameters with pair number. Solid red line represents the value estimated from the 11 pair ensemble data set and the shaded region represents the corresponding 68 per cent confidence interval.

To evaluate the validity of alternative models, we compare their costs. The final cost of the 11-pair ensemble solution just presented is ω = 0.0300 cycles pixel−1. A model with no plate spreading results in a final cost of ω = 0.0306 cycles pixel−1. Similarly, a model with no post-rifting viscoelastic relaxation results in a final cost of ω = 0.0351 cycles pixels−1. Each of these increments is significant (with at least 68-per cent confidence) because the cost of each of the poorly fitting models is greater than the critical value of ω = 0.0302 cycles/pixel. We therefore conclude that a complete description of the deformation field must account for both post-rifting viscoelastic relaxation and steady plate spreading.

To compare the viscosity estimates from the 11 pairs taken individually, we plot the vertical uplift due to rifting and subsequent viscoelastic relaxation since 1976 at a distance of 1 km away from the centre of rift axis, in Fig. 12. All pairs show uplift rates that decay with time. Pair 11 shows uplift rates that are higher than those calculated from pairs 1–10. This discrepancy could be due to the relatively short time span (420 days) and the low signal-to-noise ratio for this pair. As much as 15 per cent of the cumulative uplift observed over the 32-yr interval from 1976 to 2008 results from post-rifting viscoelastic relaxation. In general, the later pairs (10–11) result in a slightly higher estimate of lower crustal viscosity (albeit with larger uncertainties) than the earlier pairs (1–9).

Figure 12.

Vertical uplift of a point 1 km away from rift axis due to rifting and subsequent viscoelastic relaxation, calculated from best-fitting values of parameters estimated from pairs 1 − 11, individually and as an ensemble.

Figure 12.

Vertical uplift of a point 1 km away from rift axis due to rifting and subsequent viscoelastic relaxation, calculated from best-fitting values of parameters estimated from pairs 1 − 11, individually and as an ensemble.

## DISCUSSION

The 3-D numerical model that accounts for post-rifting viscoelastic relaxation following the Krafla Fires, steady plate spreading and inflation/deflation of the shallow magma chambers explains most of the InSAR observations between 1993 and 2008. A viscosity of ηL = 19 − 49 EPa.s ($$\eta ^{*}_{L}=19.47^{+0.22}_{-0.19}$$) in the lower crust between 8 and 24 km, ηT = 5–12 EPa.s ($$\eta ^{*}_{T}=18.81^{+0.28}_{-0.11}$$) between 24 and 32 km and ηM = 5–9 EPa.s ($$\eta ^{*}_{M}=18.81^{+0.14}_{-0.13}$$) below 32 km provides the best fit to the observed range gradient values resampled from the 11 interferometric pairs. These viscosities are of the same order as those estimated by Pollitz & Sacks (1996) using GPS observations made between 1987 and 1992, although they are much higher than the ∼1 EPa.s estimated for the viscoelastic half-space by Hofton & Foulger (1996). We also compare our viscosity estimates to those determined from post-glacial rebound, south of NVZ near the Vatnajokull ice cap. Using GPS observations made between 1996 and 2004, following glacial isostatic rebound that started in 1890, in combination with an elastic-over-viscoelastic half-space model, Pagli et al. (2007) estimated a half-space viscosity of 4–10 EPa.s. A similar result was obtained by Auriac et al. (2013) using InSAR measurements made between 1995 and 2009. This difference in the viscosity of ductile layers beneath the NVZ and Vatnajokull could be due to lateral variations in the crustal thickness (Allen et al.2002) and/or a depth-dependent viscosity structure (Yamasaki & Houseman 2012). Our estimates, however, are closer to those of Árnadóttir et al. (2009), who used GPS observations made between 1993 and 2004. Their preferred model, which explains most of the uplift signal in central and southeast Iceland due to glacial isostatic rebound, has a viscosity of 100 EPa.s between depths of 10 and 40 km, and 10 EPa.s below 40 km.

The estimated locking depth of 8.0–9.5 km is greater than the value of 5.0 ± 2.0 km estimated for the NVZ by Árnadóttir et al. (2009). It is slightly closer to the depth of $$6.3^{+1.7}_{-1.2}$$ km estimated for the TFZ in the north from GPS observations made between 2006 and 2010 (Metzger et al.2011). Since the study of Árnadóttir et al. (2009) is based on a few GPS stations, the deformation near the rift is not well resolved. However, both of these GPS-based studies: (i) neglect the contribution from post-rifting viscous relaxation following the Krafla Fires and (ii) assume a uniform elastic half-space, unlike our model where the elastic properties vary with depth. A greater locking depth, corresponding to a broader interrifting ‘deformation-accommodating zone’, is also supported by the fact that interrifting deformation at NVZ is not accommodated by Krafla alone but also by the Theistareykir and Fremrinamar volcanic rift systems subparallel to it on either side.

In order to verify our model, we compare its predictions to independent measurements from precise levelling, tiltmeters, and GPS surveys. Fig. 13 shows the cumulative uplift recorded at the Krafla power plant by levelling and tiltmeters during the rifting episode between 1975 and 1981 (Tryggvason 1984; Wright et al.2012). The modelled uplift, calculated at the power plant at geographic coordinates N65.703°, W16.775° and shown by the dashed black line, is in good agreement with the observed uplift (solid red line). The modelled uplift occurs steadily because the rifting rate is assumed to be constant between 1975 and 1981, unlike the observed uplift, which occurs during multiple episodes of inflation and deflation (Ewart et al.1990). We also assume that cumulative observed uplift is only due to the opening of the rift, and any uplift or subsidence due to charging and subsequent discharging of magma chamber during and between diking averages to zero. Metzger et al. (2012) have used GPS to measure vertical deformation rates between 1997–2010, north of the Krafla caldera. Their results, which have been corrected for seasonal oscillations and transient uplift at Theistareykir, are shown in Fig. 14 (black arrows) along with results calculated using our model (green arrows, red to blue shading) that includes deformation due to plate spreading and post-rifting viscoelastic relaxation. Our modelled vertical velocities, averaged over 1997–2010, are in good agreement with those observed by GPS. The pattern of the modelled deformation (e.g. uplift on either side near the rift and subsidence about 60 km to the west, and to the north of rift) as well as its magnitude, are consistent with the observations. A model with plate spreading alone, that does not account for viscoelastic relaxation following the Krafla Fires, produces subsidence and cannot explain the observed uplift near the rift.

Figure 13.

Comparison of observed (solid red line) and modelled (dashed black line) uplift at the Krafla power plant over 1975–1981.

Figure 13.

Comparison of observed (solid red line) and modelled (dashed black line) uplift at the Krafla power plant over 1975–1981.

Figure 14.

Vertical velocities in mm yr−1 as observed by GPS (black arrows, Metzger et al.2012) and modelled using the best fitting estimates of the parameters (green arrows, red to blue shading), averaged over 1997–2010.

Figure 14.

Vertical velocities in mm yr−1 as observed by GPS (black arrows, Metzger et al.2012) and modelled using the best fitting estimates of the parameters (green arrows, red to blue shading), averaged over 1997–2010.

The good fit of our model to the InSAR observations indicates a relatively stronger lower crust at the NVZ, assuming that viscosity decreases with depth and neglecting any lateral variations in viscosity which might result due to cooling of the lithosphere with time, away from the rift or due to variation in the geometry of the plate boundary, along strike (Pedersen et al.2009). Ideally, a model with more free parameters, where each finite element in the discretized mesh has its own viscosity estimate, would yield a better fit to the data. Since calculating the partial derivatives requires evaluating the exact fitting function at each iteration, the inversion procedure would become prohibitively expensive for more than a few free parameters. Even if we were to limit the number of free parameters by imposing a constraint on the smoothness of the viscosity distribution, such a solution is beyond the scope of this study. The viscosity estimates are also sensitive to the distribution of tensile opening in the co-rifting model which we have fixed a priori in order to avoid a trade-off between parameters. We also neglect non-Newtonian rheologies in order to avoid the computationally expensive step of ‘spinning up’ the 3-D model. The assumption of a linear rheology is justified because the observed deformation occurs at least a decade after the end of rifting (Freed & Bürgmann 2004).

As recently reviewed by Wright et al. (2012), the Krafla Fires rifting episode bears some similarities to the rifting episodes in Afar. For example, the deformation following rifting at the Asal-Ghoubbet rift in Djibouti can be described by a model that primarily includes steady dike opening, along with some viscoelastic relaxation (Cattin et al.2005), or by viscoelastic relaxation alone (Vigny et al.2007). Similarly, the rifting episode that began in 2005 at the Dabbahu-Manda Hararo rift segment shows complex behaviour (Nooner et al.2009; Grandin et al.2010; Ali et al.2012).

## CONCLUSIONS

Using a strategy that utilizes the gradient of range change along with a 3-D numerical model and a non-linear inversion scheme, we estimate geophysical parameters for the NVZ in Iceland from SAR data acquired between 1993 and 2008. The results indicate that the observed deformation can be explained by a combination of three processes, including: (i) secular plate spreading, (ii) viscoelastic relaxation following the Krafla Fires rifting episode (1975–1984) and (iii) deflation of a shallow magma chamber beneath the central volcano. The best-fitting model favours a lower crustal viscosity of 19 to 49 EPa.s, a mantle viscosity of 5 to 9 EPa.s and a locking depth of 8.0 to 9.5 km. These results demonstrate that post-rifting viscoelastic relaxation at spreading centres can last for decades, dominating the deformation field near the rift.

We thank Cliff Thurber, Chuck DeMets, Hélène Le Mével, Eric Calais, Tim Wright and Thora Arnadottir for constructive discussions. We also thank two anonymous reviewers for helpful comments that improved the clarity of the manuscript. Data acquired by ERS and Envisat satellites were provided by the European Space Agency under the terms of Category-I projects. Freysteinn Sigmundsson acknowledges support from the University of Iceland Research Fund and the Icelandic Research Fund. Computing resources for benchmarking and testing Defmod were provided under XSEDE award TG-EAR110020 to S. Tabrez Ali. This research was supported in part by the U.S. National Science Foundation award EAR-0810134 and a grant from the Graduate School at the University of Wisconsin-Madison.

*
Now at: School of Earth and Space Exploration, Arizona State University, Tempe, AZ 85281, USA.
Now at: Department of Geology and Geological Engineering, South Dakota School of Mines & Technology, Rapid City, SD 57701, USA.

## REFERENCES

Al-Attar
D.
Tromp
J.
Geophys. J. Int.
,
2013

doi:10.1093/gji/ggt395
Ali
S.T.
Defmod – Finite element code for modeling crustal deformation
2011a

Available at: http://defmod.googlecode.com, last accessed 1 May 2012.
Ali
S.T.
Par-SA – Parallel simulated annealing using MPI
2011b

Ali
S.T.
Feigl
K.L.
A new strategy for estimating geophysical parameters from InSAR data: application to the Krafla central volcano in Iceland
Geochem. Geophys. Geosyst.
,
2012
, vol.
13

doi:10.1029/2012GC004112
Ali
S.T.
Freed
A.M.
Contemporary deformation and stressing rates in Southern Alaska
Geophys. J. Int.
,
2010
, vol.
183

2
(pg.
557
-
571
)
Ali
S.T.
Hamling
I.
Feigl
K.L.
Calais
E.
Wright
T.
Geodetic measurements and numerical models of the Afar rifting sequence 2005-2010. Abstract G42A-04 presented at 2012 Fall Meeting, American Geophysical Union (AGU), San Francisco, California, USA
2012
Allen
R.M.
, et al.  .
Plume-driven plumbing and crustal formation in Iceland
J. geophys. Res.
,
2002
, vol.
107

doi:10.1029/2001JB000584
T.
Lund
B.
Jiang
W.
H.
Björnsson
H.
P.
Sigurdsson
T.
Glacial rebound and plate spreading: results from the first countrywide GPS observations in Iceland
Geophys. J. Int.
,
2009
, vol.
177

2
(pg.
691
-
716
)
Arnason
K.
A Digital Elevation Model for Iceland with 25-m Grid Posting and 10-m Accuracy
,
2006
Reykjavik, Iceland
Icelandic Land Survey
Auriac
A.
Spaans
K.H.
Sigmundsson
F.
Hooper
A.
Schmidt
P.
Lund
B.
Iceland rising: solid earth response to ice retreat inferred from satellite radar interferometry and visocelastic modeling
J. geophys. Res.
,
2013
, vol.
118

doi: 10.1002/jgrb.50082
Balay
S.
, et al.  .
PETSc (portable, extensible toolkit for scientific computation)
2012

Available at: http://www.mcs.anl.gov/petsc, last accessed 1 May 2012.
Brandsdóttir
B.
Menke
W.
P.
White
R.S.
Staples
R.K.
Färoe-Iceland ridge experiment 2. Crustal structure of the Krafla central volcano
J. geophys. Res.
,
1997
, vol.
102

B4
(pg.
7867
-
7886
)
Buck
W.R.
P.
Brandsdottir
B.
Tectonic stress and magma chamber size as controls on dike propagation: constraints from the 1975-1984 Krafla rifting episode
J. geophys. Res.
,
2006
, vol.
111

doi:10.1029/2005JB003879
Bürgmann
R.
Dresen
G.
Rheology of the lower crust and upper mantle: evidence from rock mechanics, geodesy, and field observations
Annu. Rev. Earth Planet. Sci.
,
2008
, vol.
36
(pg.
531
-
567
)
Carr
B.B.
Geodetic measurements and numerical models of rifting in Northern Iceland for 1993-1999
MS thesis
,
2008
Cattin
R.
Doubre
C.
de Chabalier
J.B.
King
G.
Vigny
C.
Avouac
J.P.
Ruegg
J.C.
Numerical modelling of quaternary deformation and post-rifting displacement in the Asal-Ghoubbet rift (Djibouti, Africa)
Earth planet. Sci. Lett.
,
2005
, vol.
239

3
(pg.
352
-
367
)
Cayol
V.
Cornet
F.H.
Effects of topography on the interpretation of the deformation field of prominent volcanoes: application to Etna
Geophys. Res. Lett.
,
1998
, vol.
25

11
(pg.
1979
-
1982
)
Darbyshire
F.A.
Priestley
K.F.
White
R.S.
Stefánsson
R.
Gudmundsson
G.B.
Jakobsdóttir
S.S.
Crustal structure of central and northern Iceland from analysis of teleseismic receiver functions
Geophys. J. Int.
,
2000
, vol.
143

1
(pg.
163
-
184
)
DeMets
C.
Gordon
R.G.
Argus
D.F.
Geologically current plate motions
Geophys. J. Int.
,
2010
, vol.
181

1
(pg.
1
-
80
)
DiCaprio
C.J.
Measuring and modeling viscoelastic relaxation of the lithosphere with application to the northern volcanic zone, Iceland
PhD thesis
,
2010
California Institute of Technology
Dubois
L.
Feigl
K.L.
Komatitsch
D.
T.
Sigmundsson
F.
Three-dimensional mechanical models for the June 2000 earthquake sequence in the south Iceland seismic zone
Tectonophysics
,
2008
, vol.
457

1
(pg.
12
-
29
)
Duchossois
G.
Guignard
J.-P.
Proposed uses of ERS-1
,
1987
, vol.
7

11
(pg.
293
-
298
)
Efron
B.
Tibshirani
R.
Bootstrap methods for standard errors, confidence intervals, and other measures of statistical accuracy
Stat. Sci.
,
1986
, vol.
1

1
(pg.
54
-
75
)
P.
Plate boundaries, rifts and transforms in Iceland
Jökull, Icelandic J. Earth Sci.
,
2008
, vol.
58
(pg.
35
-
58
)
Ewart
J.A.
Voight
B.
Bjornsson
A.
Ryan
M.
Dynamics of Krafla caldera, north Iceland: 1975–1985
Magma Transport and Storage
,
1990
John Wiley and Sons
(pg.
225
-
276
)
Feigl
K.L.
Thatcher
W.
Geodetic observations of post-seismic transients in the context of the earthquake deformation cycle
Comptes Rendus - Geosci.
,
2006
, vol.
338
(pg.
1012
-
1028
)
Feigl
K.L.
Thurber
C.H.
A method for modelling radar interferograms without phase unwrapping: application to the M 5 Fawnskin, California earthquake of 1992 December 4
Geophys. J. Int.
,
2009
, vol.
176

2
(pg.
491
-
504
)
Foulger
G.R.
Jahn
C.H.
Seeber
G.
P.
Julian
B.R.
Heki
K.
Post-rifting stress relaxation at the divergent plate boundary in Northeast Iceland
Nature
,
1992
, vol.
358

6386
(pg.
488
-
490
)
Freed
A.M.
Bürgmann
R.
Evidence of power-law flow in the Mojave desert mantle
Nature
,
2004
, vol.
430

6999
(pg.
548
-
551
)
Goffe
W.L.
Simann: a global optimization algorithm using simulated annealing
Stud. Nonlinear Dyn. Economet..
,
1996
, vol.
1

3
(pg.
169
-
176
)
Grandin
R.
Socquet
A.
Doin
M.-P.
Jacques
E.
de Chabalier
J.-B.
King
G.
Transient rift opening in response to multiple dike injections in the Manda Hararo rift (Afar, Ethiopia) imaged by time-dependent elastic inversion of interferometric synthetic aperture radar data
J. geophys. Res. (1978–2012)
,
2010
, vol.
115

B9
Heki
K.
Foulger
G.R.
Julian
B.R.
Jahn
C.H.
Plate dynamics near divergent boundaries: geophysical implications of postrifting crustal deformation in NE Iceland
J. geophys. Res.
,
1993
, vol.
98
(pg.
14279
-
14297
)
Hetland
E.A.
Hager
B.H.
Interseismic strain accumulation: spin-up, cycle invariance, and irregular rupture sequences
Geochem. Geophys. Geosyst.
,
2006
, vol.
7

doi:10.1029/2005GC001087
Hofton
M.A.
Foulger
G.R.
Postrifting anelastic deformation around the spreading plate boundary, north Iceland, 1. Modeling of the 1987-1992 deformation field using a viscoelastic Earth structure
J. geophys. Res.
,
1996
, vol.
101
(pg.
25
-
25
)
Hu
Y.
Wang
K.
He
J.
Klotz
J.
G.
Three-dimensional viscoelastic finite element model for postseismic deformation of the great 1960 Chile earthquake
J. geophys. Res.
,
2004
, vol.
109

doi:10.1029/2004JB003163
Johannesson
H.
Tectonic Map of Iceland (1:600,000 scale)
,
2009
Mal Og Menning
Kirkpatrick
S.
Gelatt
C.D.
Jr.
Vecchi
M.P.
, et al.  .
Optimization by simmulated annealing
Science
,
1983
, vol.
220

4598
(pg.
671
-
680
)
Malvern
L.E.
Introduction to the Mechanics of a Continuum Medium
,
1969
Prentice-Hall
Massonnet
D.
Feigl
K.L.
Radar interferometry and its application to changes in the earth's surface
Rev. Geophys.
,
1998
, vol.
36

4
(pg.
441
-
500
)
Massonnet
D.
Rabaute
T.
IEEE Trans. Geosci. Remote Sens.
,
1993
, vol.
31

2
(pg.
455
-
464
)
Masterlark
T.
Magma intrusion and deformation predictions: sensitivities to the Mogi assumptions
J. geophys. Res.
,
2007
, vol.
112

B6
pg.
6419

McLeod
I.H.
Cumming
I.G.
Seymour
M.S.
Envisat ASAR data reduction: impact on SAR interferometry
IEEE Trans. Geosci. Remote Sens.
,
1998
, vol.
36

2
(pg.
589
-
602
)
Melosh
H.J.
Raefsky
A.
The dynamical origin of subduction zone topography
Geophys. J. R. astr. Soc.
,
1980
, vol.
60

3
(pg.
333
-
354
)
Metzger
S.
Seismic potential of the Húsavík-Flatey fault and kinematics of the Tjörnes Fracture Zone, North Iceland, studied using InSAR and GPS
PhD thesis
,
2012
Eidgenössische Technische Hochschule Zürich
Metzger
S.
Jónsson
S.
Danielsen
G.
Hreinsdóttir
S.
Jouanne
F.
Giardini
D.
Villemin
T.
Present kinematics of the Tjörnes Fracture Zone, North Iceland, from campaign and continuous GPS measurements
Geophys. J. Int.
,
2012
, vol.
192
(pg.
441
-
455
)
Metzger
S.
Jónsson
S.
H.
Locking depth and slip-rate of the Húsavík Flatey fault, North Iceland, derived from continuous GPS data 2006–2010
Geophys. J. Int.
,
2011
, vol.
187
(pg.
564
-
576
)
Mogi
K.
Relations between the eruptions of various volcanoes and the deformations of the ground surfaces around them
Bull. Earthquake Res. Inst.
,
1958
, vol.
36
(pg.
99
-
134
)
Nooner
S.L.
Bennati
L.
Calais
E.
Buck
W.R.
Hamling
I.J.
Wright
T.J.
Lewi
E.
Post-rifting relaxation in the Afar region, Ethiopia
Geophys. Res. Lett.
,
2009
, vol.
36
pg.
L21308

Y.
Surface deformation due to shear and tensile faults in a half-space
Bull. seism. Soc. Am.
,
1985
, vol.
75

4
(pg.
1135
-
1154
)
Pagli
C.
Sigmundsson
F.
Lund
B.
Sturkell
E.
H.
P.
T.
Hreinsdóttir
S.
Glacio-isostatic deformation around the Vatnajökull ice cap, Iceland, induced by recent climate warming: GPS observations and finite element modeling
J. Geophys. Res.
,
2007
, vol.
112

doi:10.1029/2006JB004421
Pedersen
R.
Sigmundsson
F.
Masterlark
T.
Rheologic controls on inter-rifting deformation of the Northern Volcanic Zone, Iceland
Earth planet Sci. Lett.
,
2009
, vol.
281

1
(pg.
14
-
26
)
Pollitz
F.F.
Sacks
I.S.
Viscosity structure beneath northeast Iceland
J. geophys. Res.
,
1996
, vol.
101

B8
(pg.
17771
-
17793
)
Rennen
M.
Basics on coordinates and their reference
Landmælingar Íslands
,
2002
Sandwell
D.T.
Price
E.J.
Phase gradient approach to stacking interferograms
J. geophys. Res.
,
1998
, vol.
103

B12
(pg.
30183
-
30204
)
Segall
P.
Earthquake and Volcano Deformation
,
2010
Princeton Univ. Press
Sigmundsson
F.
Iceland Geodynamics: Crustal Deformation and Divergent Plate Tectonics
,
2006
Springer-Verlag
Tryggvason
E.
Widening of the Krafla fissure swarm during the 1975–1981 volcano-tectonic episode
Bull. Volcanol.
,
1984
, vol.
47

1
(pg.
47
-
69
)
Vigny
C.
de Chabalier
J.-B.
Ruegg
J.-C.
Huchon
P.
Feigl
K.L.
Cattin
R.
Asfaw
L.
Kanbari
K.
Twenty-five years of geodetic measurements along the Tadjoura-Asal rift system, Djibouti, East Africa
J. geophys. Res. (1978–2012)
,
2007
, vol.
112

B6
Wright
T.J.
, et al.  .
Geophysical constraints on the dynamics of spreading centres from rifting episodes on land
Nat. Geosci.
,
2012
, vol.
5

4
(pg.
242
-
250
)
Yamasaki
T.
Houseman
G.A.
The signature of depth-dependent viscosity structure in post-seismic deformation
Geophys. J. Int.
,
2012
, vol.
190
(pg.
769
-
784
)
Zeeuw-van Dalfsen
E.
Pedersen
R.
Sigmundsson
F.
Pagli
C.
Satellite radar interferometry 1993–1999 suggests deep accumulation of magma near the crust–mantle boundary at the Krafla volcanic system, Iceland
Geophys. Res. Lett.
,
2004
, vol.
112

doi:10.1029/2004GL020059