## 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 =10^{18} 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.

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).

Pair | Sat. | Ha | Species | First epoch (t_{i}) | Second epoch (t_{j}) | Span | ||
---|---|---|---|---|---|---|---|---|

(m) | Orbit | Year | Orbit | Year | (days) | |||

1 | ERS | −2628.2 | A | 10 174 | 1993.4822 | 17398 | 1998.6274 | 1879 |

2 | ERS | −53.8 | A | 10 174 | 1993.4822 | 5875 | 1996.4235 | 1074 |

3 | ERS | −63.1 | A | 10 675 | 1993.5781 | 11386 | 1997.4767 | 1424 |

4 | ERS | 64.9 | A | 10 675 | 1993.5781 | 23410 | 1999.7781 | 2264 |

5 | ERS | −198.6 | B | 11 677 | 1993.7699 | 22408 | 1999.5863 | 2124 |

6 | ERS | 70.3 | B | 11 677 | 1993.7699 | 6877 | 1996.6147 | 1039 |

7 | ERS | 190.2 | A | 5875 | 1996.4235 | 23410 | 1999.7781 | 1225 |

8 | ERS | 392.6 | B | 6376 | 1996.5191 | 22408 | 1999.5863 | 1120 |

9 | ERS | 82.2 | B | 6877 | 1996.6147 | 17899 | 1998.7233 | 770 |

10 | ENV | 268.8 | C | 18 099 | 2005.6219 | 28620 | 2007.6356 | 735 |

11 | ENV | −370.6 | C | 22 608 | 2006.4849 | 28620 | 2007.6356 | 420 |

Pair | Sat. | Ha | Species | First epoch (t_{i}) | Second epoch (t_{j}) | Span | ||
---|---|---|---|---|---|---|---|---|

(m) | Orbit | Year | Orbit | Year | (days) | |||

1 | ERS | −2628.2 | A | 10 174 | 1993.4822 | 17398 | 1998.6274 | 1879 |

2 | ERS | −53.8 | A | 10 174 | 1993.4822 | 5875 | 1996.4235 | 1074 |

3 | ERS | −63.1 | A | 10 675 | 1993.5781 | 11386 | 1997.4767 | 1424 |

4 | ERS | 64.9 | A | 10 675 | 1993.5781 | 23410 | 1999.7781 | 2264 |

5 | ERS | −198.6 | B | 11 677 | 1993.7699 | 22408 | 1999.5863 | 2124 |

6 | ERS | 70.3 | B | 11 677 | 1993.7699 | 6877 | 1996.6147 | 1039 |

7 | ERS | 190.2 | A | 5875 | 1996.4235 | 23410 | 1999.7781 | 1225 |

8 | ERS | 392.6 | B | 6376 | 1996.5191 | 22408 | 1999.5863 | 1120 |

9 | ERS | 82.2 | B | 6877 | 1996.6147 | 17899 | 1998.7233 | 770 |

10 | ENV | 268.8 | C | 18 099 | 2005.6219 | 28620 | 2007.6356 | 735 |

11 | ENV | −370.6 | C | 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.

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

*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.

### 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 η* = Log_{10}(η/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 *D*_{l}, 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 *D*_{l} and the height of the ‘step’ scales with the long-term plate spreading rate.

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 *D*_{l}, 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 *D*_{K} is estimated during the inversion. For the magma chamber beneath Theistareykir, the source coordinates, including depth *D*_{T}, 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 (*D*_{K} ∼ *D*_{T} ∼ 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 [*t*_{i}, *t*_{j}]. 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

*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:

*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 *D*_{l} 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 *D*_{K}, its rate of volume change $$\Delta \dot{V}_{K}$$ and depth *D*_{T}, 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).

## 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 10^{18} 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).

Parameter name | $$\eta ^{*}_{L}$$ | $$\eta ^{*}_{T}$$ | $$\eta ^{*}_{M}$$ | D_{l} | D_{K} | $$\Delta \dot{V}_{K}$$ | D_{T} | $$\Delta \dot{V}_{T}$$ |
---|---|---|---|---|---|---|---|---|

units | Log_{10}(η/1 Pa.s) | (km) | (km) | ( × 10^{5}m^{3} yr^{−1}) | (km) | ( × 10^{5} m^{3} 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}$$ | D_{l} | D_{K} | $$\Delta \dot{V}_{K}$$ | D_{T} | $$\Delta \dot{V}_{T}$$ |
---|---|---|---|---|---|---|---|---|

units | Log_{10}(η/1 Pa.s) | (km) | (km) | ( × 10^{5}m^{3} yr^{−1}) | (km) | ( × 10^{5} m^{3} 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.

Parameter name | $$\eta ^{*}_{L}$$ | $$\eta ^{*}_{T}$$ | $$\eta ^{*}_{M}$$ | D_{l} | D_{K} | $$\Delta \dot{V}_{K}$$ | D_{T} | $$\Delta \dot{V}_{T}$$ |
---|---|---|---|---|---|---|---|---|

units | Log_{10}(η/1 Pa.s) | (km) | (km) | ( × 10^{5}m^{3} yr^{−1}) | (km) | ( × 10^{5} m^{3} 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.1^{1} | 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}$$ | D_{l} | D_{K} | $$\Delta \dot{V}_{K}$$ | D_{T} | $$\Delta \dot{V}_{T}$$ |
---|---|---|---|---|---|---|---|---|

units | Log_{10}(η/1 Pa.s) | (km) | (km) | ( × 10^{5}m^{3} yr^{−1}) | (km) | ( × 10^{5} m^{3} 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.1^{1} | 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.

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 *D*_{l} 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.

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).

## 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.

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.

## REFERENCES

## SUPPORTING INFORMATION

Additional Supporting Information may be found in the online version of this article:

**Figure S1.** Part of the finite element mesh used for calculating deformation following Krafla Fires. Colours represent the east component of the displacement field at the end of rifting in 1981.

**Figure S2.** Interferograms showing (a) observed, (b) modelled and (c) residual values of wrapped phase change for pairs 2–6. Modelled values are calculated from the final estimate of the parameters. Residual values are calculated by subtracting the modelled values from the observed values. 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.

**Figure S3.** Interferograms showing (a) observed, (b) modelled and (c) residual values of wrapped phase change for pairs 7–11. 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. S2 (http://gji.oxfordjournals.org/lookup/suppl/doi:10.1093/gji/ggt462/-/DC1).

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.