Improved source modelling through combined use of InSAR and GPS under consideration of correlated data errors: application to the June 2000 Kleifarvatn earthquake, Iceland

Simultaneous use of multiple independent data sets can improve constraints on earthquake source-model parameters. However, the ways in which data sets have been combined in the past are manifold and usually qualitative. In this paper we present a method to combine geodetic data in source model estimations, which includes characterizing the data errors and estimating realistic model-parameter uncertainties caused by these errors. We demonstrate this method in a case study of the June 2000 Kleifarvatn earthquake, which occurred on Reykjanes Peninsula in Iceland. We begin by showing to what extent additional data can positively inﬂuence the source modelling results, by combining both GPS and descending-orbit InSAR data, which were used in two earlier studies of that event, with InSAR data from an ascending orbit. We estimate the data error covariances of the InSAR observations and base the data weights in our model-parameter optimization on the corresponding data variance–covariance matrix. We also derive multiple sets of synthetic data errors from the estimated data covariances that we use to modify the original data to generate numerous data realizations. From these data realizations we estimate the model-parameter uncertainties. We ﬁrst model the Kleifarvatn earthquake as a simple uniform-slip fault and subsequently as a fault with variable slip and rake. Our fault model matches well with the ﬁeld observations of coseismic surface ruptures and its near-vertical dip (83 ◦ ) agrees with the regional faulting style as well as with aftershock locations. The two published source models of the event, on the other hand, both differ from our model as well as differing for one another. These studies, which were based on the descending InSAR data alone (the ﬁrst study) and on that same data and GPS data (the second study), both neglect correlations in the InSAR data and do not report model-parameter uncertainties. Therefore, to compare these results with our model, we simulate the earlier model estimation set-ups and provide realistic estimates of the model-parameters uncertainties for these cases. We then discuss the signiﬁcance of the difference between the existing fault models and demonstrate that both the inclusion of additional independent data as well as the covariance-based data weights improve the model-parameter estimation.

distance along the line-of-sight (LOS) of the radar, and thus provide only one component of the surface displacement vector. Due to the different viewing geometries of ascending-and descending-pass In-SAR measurements, one can observe two different projections of the displacement vector. Wright et al. (2004) demonstrate the advantage of simultaneously using data from both look directions in narrowing model parameter confidence intervals, with respect to the use of single interferograms. Beside InSAR and GPS data, offsets between two SAR images, measured parallel to the radar flight direction (azimuth offsets), can provide valuable information about azimuthal surface movements, which is orthogonal to the line-of-sight displacement vector imaged with InSAR. However, the applicability of this method is limited to large surface displacements as its noise level is much higher than for standard InSAR and GPS observations. In summary, while GPS data are spatially incomplete, they provide full 3-D surface displacement vectors, and while InSAR only give a 1-D projection of the surface displacements, they can provide spatially dense observations. Several authors have therefore taken advantage of these different data characteristics and combined the different data sets in source estimations. For example, Johanson et al. (2006) combined InSAR from different look directions and GPS in their source imaging of the 2004 Parkfield earthquake. In a study of the 1999 Hector Mine earthquake Jónsson et al. (2002) used, in addition to GPS data, ascending and descending InSAR data, azimuth offsets for the source imaging.
The way one combines different geodetic data sets, such as In-SAR and GPS, can significantly influence the result of source model estimations. A variety of data weighting approaches has been published, which include using no specified weights or equal weights (e.g. Wright et al. 2003;Funning et al. 2007), arbitrary weights for InSAR (e.g. Delouis et al. 2002;Árnadóttir et al. 2004), weights that minimize the misfit (e.g. Schmidt et al. 2005), weights that reflect the area each data point represents in subsampled InSAR data (e.g. Simons et al. 2002;Lasserre et al. 2005), weight factors based on the data error variance (e.g. Jónsson et al. 2002;Pedersen et al. 2003), and data weights based on both the data error variance and covariance (e.g. Fukushima & Cayol 2005;Lohman et al. 2005). Only the last two approaches quantitatively consider the data quality of the independent data sets, but since InSAR and GPS data exhibit correlated data errors (Hanssen 2001), using only the variance is not an adequate description of the data errors. Data weights that also account for the data error covariances do not only balance the independent data sets with respect to one another, but also each single data point consistently across the entire data set. Consistently, because the complete representation of data errors enables the propagation of the error structure through the different processing steps, such as data subsampling that commonly is carried out for InSAR data.
The weighting based on the data errors characteristics requires an empirical estimation of the apparent error structure. The error of GPS displacements is often estimated from the variance of repeated measurements, which is not applicable in case of InSAR. A few different methods have been used to retrieve estimates of the autocovariance of data errors in interferograms. Hanssen (2001) found the characteristics of InSAR errors in many interferograms to be similar, apart from a scale factor, and to be appropriately represented by power-law functions according to elementary turbulence theory. In his approach only this scale factor needs to be estimated from the interferogram. Knospe & Jónsson (2008) estimated 2-D covariance functions for InSAR data exhibiting strong anisotropic atmospheric errors. Lohman et al. (2005) discuss, in addition to atmospheric errors, the contribution of decorrelation and other ef-fects to InSAR data errors introduced in the processing, for example, from filtering and multilooking. In purely empirical autocovariance functions these effects are included while no assumptions on the error origin are required.
In this paper we demonstrate the effect of combining different geodetic data on source-model estimations and the corresponding model-parameter uncertainties, as well as addressing the role of data weights. In a case study we use data from the June 2000 Kleifarvatn earthquake in southwest Iceland. This choice is motivated by the fact that little is known from field and seismological data about this event near the city of Reykjavík and that the two existing geodetic source models of this event by Pagli et al. (2003Pagli et al. ( ) andÁrnadóttir et al. (2004 are somewhat different. In these two studies descending InSAR data alone , and together with GPS data (Árnadóttir et al. 2004), were used, resulting in two sets of model parameters that appear to be significantly different, although no model parameters uncertainties were reported in these studies. We complement the previously published data sets by adding ascending InSAR data with the intention of improving the source model of the Kleifarvatn earthquake. To find appropriate weights for the three independent data sets in the source model optimization we estimate the data error covariances in the interferograms. From these empirical covariances we also compute synthetic realizations of data errors and use them to estimate the uncertainties of the optimum model parameters. We then attempt to reproduce the two earlier source model results of the Kleifarvatn event Árnadóttir et al. 2004) to estimate the model parameter uncertainties for these two case studies. With these results we address several questions, such as: Are the differences in the source models significant? Can we improve the Kleifarvatn source model with the use of additional data? How do the data weights affect the source model estimations? At last, we also present the first variable slip and rake model for the Kleifarvatn earthquake.

T H E J U N E 0 0 0 K L E I FA RVAT N E A RT H Q UA K E
The Kleifarvatn earthquake occurred on 2000 June 17 on Reykjanes Peninsula in southwest Iceland, about 20 km south of the city of Reykjavík (Fig. 1). It ruptured a previously unknown fault east of Lake Kleifarvatn and was a part of a seismic sequence triggered by the magnitude 6.5 Holt earthquake ( Fig. 1), which occurred 77 km to the east in the South Icelandic Seismic Zone (SISZ). Four earthquakes of that sequence reached magnitudes of 5 or larger: the Second Holt event, the Hvalhnúkur event, the Kleifarvatn event, and the Núpshlí arháls event Clifton et al. 2003;Arnadóttir et al. 2004;Hjaltadóttir & Vogfjör 2005). From the timing of the Hvalhnúkur and Kleifarvatn events, 26 and 30 s after the Holt main shock, these events are thought to have been triggered dynamically, as the timing roughly corresponds to the shear wave traveltime from the main shock epicentre (Antonioli et al. 2006).
The Kleifarvatn earthquake was not felt as a separate event and remained undetected for about one year. The local seismic network recorded complex seismograms containing superimposed wave trains of three M > 5 events, which, along with an unfortunate coincidence of instrument malfunction and exceedance of the dynamic range at stations installed on the Reykjanes Peninsula, prohibited standard seismological analysis of the data (Clifton et al. 2003;K. Vogfjör personal communication, 2007). The Kleifarvatn event was eventually discovered in ERS interferograms, when clear coseismic deformation pattern was detected in the vicinity of in UTM coordinates (UTM zone 27V). Shown with stars are the M ∼ 5 earthquake epicentres of the June 2000 seismic sequence: the main shock ('J17'), the dynamically triggered Hvalhnúkur ('H') and Kleifarvatn ('K') events and the aftershocks: the Núpshlí arháls event ('N'), the Second Holt event ('H2') and the June 21 main shock ('J21'). Open circles mark locations of the campaign GPS sites used in this study. The inset shows the map location in southwest Iceland (solid box) and the main volcanic and tectonic structures. The on-land spreading segments (black lines), the Northern (NVZ), the Eastern (EVZ), and the Western Volcanic Zone (WVZ), are connected to the oceanic Reykjanes Ridge (RR) and the Kolbeinsey Ridge (KR) through the southern transform zones Reykjanes Peninsula (RP) and the South Icelandic Seismic Zone (SISZ) and the Tjörnes Fracture Zone (TFZ) in the north (grey lines). Lake Kleifarvatn by Pagli et al. (2003). The surface fault rupture appears to be similar to the other north-south striking, right-lateral strike-slip faults on both the Reykjanes Peninsula and in the South Icelandic Seismic Zone (Einarsson 1991;Clifton & Kattenhorn 2006).
The Kleifarvatn event has already been studied twice using geodetic data sets, resulting in two different fault models, which both have rectangular planar faults with uniform slip. The first model by Pagli et al. (2003) is based on a single descending ERS interferogram, while the second model byÁrnadóttir et al. (2004) was constrained using the same descending interferogram together with campaign GPS data. What is noteworthy about these models is their large discrepancy in fault dip (66 • to the east versus 78 • ), as well as that both models have a shallower fault dip to the east compared to the near-vertical (88 • ) distribution of relocated aftershocks (Hjaltadóttir & Vogfjör 2005) and the established understanding of the regional faulting style (Einarsson 1991). Other parameters of the two fault models are similar, except that the second fault model has somewhat more strike-slip. No estimates of the model parameter uncertainties were reported for these models, which makes it hard to judge whether or not the two fault models differ significantly from one another.

DATA
In the source optimization we use three independent data sets. Two InSAR images that provide spatially dense information about the line-of-sight (LOS) projection of the surface displacement field along both the ascending and descending imaging directions. We complement the InSAR data with the full coseismic surface dis-placement vectors at eleven locations measured using GPS and provided byÁrnadóttir et al. (2004).

InSAR data
At the high latitude of Iceland the ERS swaths are overlapping by ∼65 per cent so that the investigation area can potentially be imaged from three parallel ascending and another three parallel descending tracks. However, long winters, limited data acquisitions, and large orbital baselines limit the number of good image pairs. We formed a total of eleven coseismic interferograms; seven using ERS-2 data from ascending tracks and four using data from descending tracks. For each viewing geometry we chose the interferogram with the highest quality (Table 1).
The descending InSAR image was generated using radar scenes acquired on 2 October 1999 and 16 September 2000 and the perpendicular baseline B ⊥ between the two orbits is only 5 m (Table 1 and Fig. 2). The radar scenes of the ascending interferogram were recorded on 2 September 1999 and 17 August 2000 (Table 1 and Fig. 2). Despite its fairly small B ⊥ (38 m) and the similar time span, the resulting correlation of the phase signal is for some reason lower than in the descending interferogram. We therefore, multilooked (complex-value averaged) three adjacent pixels in the range and azimuth directions of the ascending interferogram, resulting in less white noise, but also lower resolution. The removal of the topographic phase and the transformation from radar to geographic coordinates (geocoding) are based on a digital elevation model with a resolution of about 25 m. In addition, we applied an adaptive filter (filter window size 32, filter exponent 0.8) to enhance the signalto-noise ratio (SNR) of the interferograms (Goldstein & Werner 1998). Finally, we used the snaphu software for phase unwrapping, the statistical-cost network-flow unwrapping algorithm by Chen & Zebker (2001).
The dominant feature in the interferograms is the deformation signature of the Kleifarvatn event around the lake (Fig. 2). In the descending interferogram we observe positive phase shifts (range increase) east of the lake and negative phase shifts to the southeast. The pattern in the ascending interferogram is almost antisymmetric with the largest phase shifts to the west and southwest of the lake. Such patterns are typical for a right-lateral strike-slip motion on a north-south striking fault, as the horizontal ground displacements, combined with alternating uplift and subsidence in the compressional and extensional quadrants around the fault, result in weak LOS displacement signals on the side of the fault that is farther away from the down-looking radar. Possible fault surface ruptures cannot be observed in the InSAR data, due to the near-field decorrelation (Fig. 2).
Both interferograms span time periods that include 2-3 months after the Kleifarvatn earthquake and could therefore be affected by post-seismic deformation. Rapid poro-elastic rebound was observed after the two June main shocks in the SISZ  and it is not unlikely that similar rebound occurred after the Kleifarvatn event. Post-seismic poro-elastic deformation tends to act against the vertical coseismic deformation in the near-field, although its amplitude is much smaller. Therefore, if some poroelastic deformation contaminated the coseismic interferograms it could lead to an underestimation of shallow fault slip in source modelling.
We attempted to measure the azimuth offset field by estimating the pointwise mismatch or offsets between two single-look radar images. The spatially variable offsets in the azimuth direction are determined by cross-correlating small image patches within the amplitude images (Michel 1999). No significant offsets or clear signs of the fault trace could be found using this method. With the maximum expected strike-slip on the fault of about 0.7 m (Árnadóttir et al. 2004) and an image resolution of 4.9 m in the azimuth direction, we are trying to detect a signal that is smaller than 15 per cent of the resolution cell. Michel (1999) estimated the error of such offset measurements to be somewhat lower, or around 10 per cent of the resolution. This, together with the unstable near-field conditions on the ground, implies that we are close to the limits of this technique in our case.

Estimating error statistics in InSAR data
Errors in interferograms arise from several different sources. There are noise sources at the radar instrument itself, on the path of the radar waves, at the reflecting surface and errors can be introduced in the processing of radar records as well (Hanssen 2001). Beside random white noise induced by phase decorrelation, InSAR data exhibit spatially correlated errors, due to smoothly varying atmospheric signal delays. Therefore, InSAR data errors are particular for each InSAR image depending on the state of the atmosphere and the ground surface at the time of the two radar acquisitions. As a consequence the processing steps and parameters for noise reduction, like multilooking and filtering, are usually adjusted for each interferogram.
The final power and structure of the data errors are quality measures of the processed InSAR images and the estimation of these InSAR error characteristics is important for two reasons in particular. With information about the quality of all the data we can assign meaningful and consistent weights to each data set. Furthermore, we can use the empirical covariances to generate multiple sets of synthetic data errors, which we can add to the original data, and through multiple optimizations obtain a distribution of model parameters. By this, we 'propagate' the data uncertainties to fault model uncertainties (Wright et al. 2003).
We estimate the error variances and autocovariances in an area of the interferograms where neither deformation signal is expected nor visible. We assume the error to be stationary, which implies that the error statistics estimated in non-deforming parts of the images are the same as in the adjacent deforming areas. The areas chosen for the error estimation are of about the same size as the investigation area to capture the whole bandwidth of periods present in the noise (Figs 3a and b).
We use sample semi-variogramsγ (h) to estimate the InSAR variances and sample covariogramsĈ(h) to estimate the spatial correlation in the data (Chilés & Delfiner 1999). The discrete sample semi-variogram value for distance class h c is: with N being the number of data-point pairs at locations r i and s i such that r i − s i h c . Thus, when assuming isotropic noise, the semi-variogram depends only on distance h between data points. Similarly, the sample covariogram is: (2) When calculating the sample semi-variograms and covariograms, we first remove an overall linear ramp from the error estimation image and then pick randomly data-point pairs d(r i ) and d(s i ) with distances h ranging from 40 m to 14 km. We then form the sample semi-variogramγ (h) and the sample covariogramĈ(h) by taking the average in 100 meter intervals containing about 700 single measurements. The data variance is estimated from the level at which the sample semi-variogramγ (h) forms a sill at distances larger than the correlation length. In presence of white noise, the covariance functions therefore have a step at a zero lag (Fig. 3c). For a continuous description of the covariances we fit functions to the sample covariograms. The covariance is by definition a positive-definite function and we therefore use function types ensuring positive-definiteness (Chilés & Delfiner 1999): an exponential decay of the type b · e −h/a to represent the descending covariance and an exponential decay complemented by a cosine term, b · e −h/a · cos( h c ), to account for the anti-correlation present in the error structure of the ascending interferogram. For the latter case, positive-definiteness is limited to parameter values a < c.
The ascending image has lower error variance (15 mm 2 ) than the descending image (25 mm 2 ), but shows higher autocovariance for distances smaller than 5 km (Fig. 3c). The reason for the relatively low ascending variance is that this image was multilooked in the processing. The estimated covariance functions C asc (h) and C desc (h) are We use these covariance functions to form the InSAR data covariance matrices and to generate multiple synthetic realizations of InSAR error images, see Section 4.3.

Subsampling of InSAR images
The InSAR data as shown in Fig. 2 consist of several hundred thousand data points. As the displacement field is varying smoothly, we can decimate the numerous phase values without losing important information. We subsample the unwrapped interferograms with a quadtree algorithm (Jónsson et al. 2002) to obtain a reasonable number of data points and good spatial representation of the LOS displacements. This algorithm subsequently divides the InSAR images into squares until the phase values within each box do not exceed a certain variance threshold. The average phase value of the contributing pixels is then assigned to their focal point. The algorithm is therefore sensitive to the variability of the phase values across the area and to possible data gaps. We used division-thresholds of 81 and 64 mm 2 for the ascending and descending images, respectively, which is significantly higher than the variance derived from the error estimation (Section 3.2), but leads to an appropriate representation of the deformation field with only 634 data points (Fig. 4). Data gaps result from decorrelation, masking of phase unwrapping errors, and layovers.
The subsampling procedure can be described as a linear operation. The operator A relates the full data vector d c to the subsampled data vector d: For the N i pixels averaged in the ith quadtree square, the corresponding N i values in the ith row of A are equal to 1/N i , but zero elsewhere.

GPS data
The Eleven campaign GPS sites are within our investigation area and in most cases the coseismic GPS data agree very well with the InSAR data when projected into the radar line-of-sight direction (Fig. 4). However, at two stations in the northeast the vertical components deviate strongly from the InSAR data, as well as from neighbouring GPS stations. This suggests a possible error in the measured antenna height and we therefore do not include the vertical components of these two sites in our analysis.

FAU LT M O D E L L I N G
With the combined use of ascending and descending InSAR data together with GPS data we aim to find an improved fault model of the Kleifarvatn earthquake. For the Kleifarvatn and Núpshlí arháls events we first assume simple uniform-slip faults and use a nonlinear optimization approach to find the optimum model parameters. We also estimate the corresponding model parameter uncertainties. To compare our results with earlier studies of Pagli et al. (2003) and Arnadóttir et al. (2004), we simulate these modelling calculations and estimate their corresponding model parameter uncertainties. We then discuss separately the influence of additional data and of different data weights. Finally, we take the obtained fault geometry and invert for variable slip and rake on the Kleifarvatn fault.

Data weighting
The data are a compound of data points from two InSAR images and of campaign GPS displacements and we weight the data based on their estimated uncertainty. In the fault parameter optimization we are seeking the minimum of the L2-norm: where d obs and d pred are the observed and predicted data vectors, respectively. The predicted data vector is calculated through forward modelling for a given set of fault model parameters. R is a weighting matrix and balances the data residual (d obs − d pred ) in the optimization, so that the influence of data points with high uncertainties and/or correlation is reduced. R is based on the data variance-covariance matrix with The weight of each data point is inversely proportional to the sum of its variance and its covariances with respect to all other data points.
The variance-covariance matrix of the subsampled InSAR data and GPS is formed from the variance-covariance matrix of the full data set c , which is defined by the covariance functions we estimated for the interferograms in Section 3.2. As the quadtree subsampling is an averaging procedure, the variance of each subsampled data point in d is lower than the variance of the full resolution data in d c and it depends on the quadtree square size. We propagate c to with the subsampling operator A (eq. 3): c is such a large matrix (740.000 × 740.000 values) that we calculate consecutively each element of from a subset of c , here C N i ×N j , containing the covariances between the pixels in each pair of quadtree squares, d i and d j : where N i and N j are the number of pixels averaged in d i and d j , respectively. The diagonal of (eq. 7, i = j) is therefore the average autocovariance of the pixels in the quadtree squares themselves and decreases with increasing number of pixel (or increasing area of the quadtree square). The off-diagonal values of (i = j) approximately follow C asc (h) or C desc (h) (Fig. 3c). Hence, the weights of data points in densely sampled areas are relatively low, because the quadtree squares sizes are small and they have many neighbouring squares at short distances (Fig. 5). Quadtree squares located at the border of the investigation area or near data gaps have relatively high weights because of fewer correlated data point neighbours. In summary, it is the area the subsampled data point is representing (square size) and its position with respect to other data points that control its resulting weight (Fig. 5).
In case of quadtree subsampling the point density is irregular and depends on the local gradient of the deformation signal. This gradient is usually higher in the near field of the fault and so is the point density. Disregarding the data autocovariance inherently increases the weight per area in densely sampled regions, which holds for all irregular InSAR subsampling procedures. Through the full propagation of the variance-covariance matrix (eqs 6 and 7), the assigned variance for a data point representing a small area remains high, balancing the weight per area across the irregularly sampled investigation area.

Nonlinear optimization/uniform slip models
We model the Kleifarvatn and the Núpshlí arháls events as two rectangular planar faults with uniform slip in a homogeneous, elastic half-space (Okada 1992), using a Poisson-ratio of 0.28 (Árnadóttir et al. 2004). Each model is defined by nine model parameters describing the fault location, orientation, dimensions and the mechanism. We also account for possible orbit errors and long-period atmospheric noise by adding plane parameters (3 parameters for each InSAR image) to the model vector.
To minimize the data misfit (eq. 4) we use a simulated annealing (SA) optimization approach (Cervelli et al. 2001) that samples the model parameter space within predefined bounds in search for regions with low model costs. The SA algorithm progressively samples these low-cost regions more exhaustively until identifying a set of model parameters that ideally are close to the global minimum of the misfit space. Since the algorithm uses a fixed number of model calls and has a predefined minimum step size in the parameter space, it can only find a set of model parameters that is near the absolute model-cost minimum. Therefore, we complement this procedure with a subsequent derivative-based Levenberg-Marquardt inversion, using the result of the SA optimization as a starting model and using smaller step sizes to improve the model further.
The performance of the SA algorithm in finding the global minimum in a misfit space is controlled by the so called critical temperature T c and by a time-evolving temperature function, that is, the cooling schedule. The cooling schedule usually begins with fast cooling, slowing down near T c , and then speeding up again after T c has been reached. Cervelli et al. (2001) found a Gaussian time-pertemperature function centred at T c to be appropriate for SA. In this case T c is equal to the mean temperature T . However, using the fast approach to determine T c (and thus T ) proposed by Cervelli et al. (2001), we find that a considerably lower T worked better in our case, in the sense that the scatter of the final model parameters decreases as well as their mean model cost (Fig. 6). However, some far off outliers appear in the model cost histograms, when T is too low, and in these cases the algorithm does not escape unfavourable local minima. We chose a value of 100 for T . The model cost of 500 solutions then has a narrow distribution without any outliers.
The SNR for the Kleifarvatn data is not very high, compared to most other InSAR studies of earthquakes. The deformation is caused by medium-sized earthquakes and the maximum observed LOS displacement is less than 15 cm (Fig. 4), while the maximum error amplitudes are about 3 cm (Fig. 3 top panel). According to Cervelli et al. (2001), a low SNR results in more minima in the model-misfit space. We attempted to account for this by deviating from the Gaussian cooling schedule and making the search more exhaustive, which we did by using slower cooling after T c , allowing for about 20 per cent of extra searching time. The final model-cost histograms indicate the finite precision of the optimization (Fig. 6) and they suggest that the global minimum of the misfit space is located in a broad valley with shallow slopes.
The model fault plane that best explains the measured surface displacement near Lake Kleifarvatn is NNE striking and almost vertical ( Fig. 7 and Table 2). The fault plane reaches the surface and extends 6.5 km southwards from the eastern lake shore, near the earthquake epicentre, with a strike of N5 • E. The modelled mechanism is oblique, with dextral strike slip of 0.67 m and a small amount of normal faulting (0.1 m). Based on the estimated fault parameters we calculate the moment M 0 = μAu and the moment magnitude M W = 2 3 log 10 M 0 − 6.03 of the Kleifarvatn earthquake from the fault plane area A, the static displacement u, and by assuming rock rigidity μ = 3 × 10 10 N m −2 the resulting moment is M 0 = 7.15 × 10 17 N m and the corresponding moment magnitude close to 5.9.
The majority of the residuals between observed and predicted data is smaller than 2 cm (Fig. 7), that is, comparable to the magnitude of the data errors (see Fig. 3). However, systematic data residuals are visible in the near field that likely result from slip variations along the fault and therefore cannot be explained by this simple uniform-slip model.
The deformation signal of the Núpshlí arháls event is small in comparison with the Kleifarvatn event and we cannot reliably determine its source parameters due to the superposition of the two deformation signals. Therefore, we put tight bounds on the location of the smaller event and fixed its strike to 12 • (K. Vogfjör personal communication, 2007) to stabilize the optimization. Similar to the Kleifarvatn event, the Núpshlí arháls is not well represented by a simple fault model as significant residuals remain. However, the lack of data in the near field hampers a more detailed analysis of this event. Moreover, seismological investigations in this case are pointing to a more complex source mechanism than a planar failure (K. Vogfjör personal communication, 2007).

Estimation of fault model parameter uncertainties
The uncertainties of the inferred model parameters are estimated by using the same optimization procedure and multiple data sets that have been modified by synthetic noise. First, 2500 realizations of data errors were generated, ε synth,i , based on the GPS and InSAR data covariance matrices. We then added each realization to the original data and eq. (4) becomes: From tests we know that in our case the difference between the empirical covariance and the covariance of two superposed realizations of correlated noise is negligible. Hence, the weighting of the modified data sets based on the empirical covariance functions is still meaningful. We carried out multiple fault parameter estimations of the Kleifarvatn event using all the modified data sets and formed histograms of the resulting model parameters (Fig. 8, top row). From the histograms we find the corresponding 95 per cent confidence bounds of the model parameters, or strictly speaking the 0.025-and 0.975quantiles of the samples. The fault parameter values of the optimal model correlate well with the histogram maxima (Fig. 8), except for the fault strike. For the fault width and the fault dip we observe relatively broad histogram peaks and some degree of asymmetry. The scatter plots in Fig. 8 show how one model parameter may depend on another parameter and thus enable visual detection of parameter trade-offs. Here, we observe parameter trade-offs between fault slip, fault width and fault dip.
Compared to the earlier fault models of the Kleifarvatn earthquake by Pagli et al. (2003Pagli et al. ( ) andÁrnadóttir et al. (2004, our optimal fault is significantly longer and therefore extends farther to the north. The fault length of both of the earlier models lies outside the 95 per cent confidence interval. In addition, our fault has a steeper dip and more slip than the earlier fault models.

Simulation of earlier fault model optimizations
We suspect model parameter trade-offs, for example, between fault slip and dip, to cause the discrepancy between the three source models of the Kleifarvatn event: Model 1 ), Model 2 (Árnadóttir et al. 2004), and our model ( Fig. 8 and Table 2). Model 1 is based on a single descending InSAR image, from which we can only measure the 1-D LOS projection of the 3-D coseismic displacement field. Model 2 is based on the same InSAR data as well as on several 3-D coseismic GPS displacements and it should therefore be less affected by model parameter trade-offs than Model 1. We expect our fault model to be even less sensitive to parameter trade-offs as it is also based on ascending InSAR data (Wright et al. 2004).
The final result of geodetic source optimizations can also be influenced by many methodological factors, such as by choices made in the data processing, InSAR data subsampling, data weighting, the particular optimization approach, and by assumptions made about the forward model. These assumptions and methods are similar in all the three studies of the Kleifarvatn event Arnadóttir et al. 2004, and this study), except for the data weighting. Therefore, the three fault models may possibly differ due to this methodological difference as well as due to the fact that different data sets were included in the fault model optimizations.
No model parameter uncertainties or trade-offs were reported for Model 1 and Model 2 Árnadóttir et al. 2004). Therefore, we here simulate these optimizations to estimate parameter uncertainties and dependencies of these models in the same way as we described in Section 4.3. The similarity of the applied methods and assumptions made in the optimizations enable the simulations where we separately investigate the influence of additional independent data (this section) and the impact of data weights accounting for the correlation of data errors (next Section 4.5).
Originally, Model 1 and Model 2 were estimated considering the data variances only Árnadóttir et al. 2004). To simulate these fault model optimizations we build the data variance matrix from the variance we estimated in the descending interferogram of 25 mm 2 and the variance of the GPS data. For comparison, we also apply this variance based weighting in our set-up and refer to it as simulation of Model 3. The low variance we obtained in our error analysis for the multilooked ascending image (Fig. 3) does not reflect its poor quality and we therefore estimate the full-resolution variance by extrapolating the autocovariance function resulting in 46 mm 2 . Apart from these changes in the data weighting we run the simulations of Model 1 and Model 2 with adjustments in the parametrization of the SA algorithm and with a simplified forward model. In the SA algorithm we adapt the mean temperature T for each optimization set-up, since T depends on the input data and the weighting matrix R. For simplicity, we fix the parameters (see Table 2) for the smaller Núpshlí arháls event (Figs 1 and 7), despite the fact it was a part of the optimization in all three studies, as its signal is not easy to explain with a simple fault model and the optimization calculations tend to become unstable .
In Fig. 9, we show the distributions of four selected model parameters resulting from the simulations when the autocorrelation of data errors in the data weighting is neglected. The Kleifarvatn fault model parameters published by Pagli et al. (2003Pagli et al. ( ) andÁrnadóttir et al. (2004 fall within the corresponding parameter histograms from the simulation of Model 1 and Model 2 (Fig. 9). The older fault model estimations thus seem to be reproducible, which is a requirement for drawing meaningful conclusions from this case study.
It is clear that increasing the number of independent data sets improves the precision of the estimated fault parameters. However, despite the fact that GPS data are typically an important complement to InSAR, because they provide full 3-D displacement vectors, we see here only a limited influence on the model parameter estimates (Fig. 9). For instance, the fault dip values for Model 1 show a broad distribution of more than 30 • that only slightly improves by adding the GPS data (Model 2). Including the ascending InSAR data, on the other hand, narrows the histogram by a factor of three (Model 3). The weight of each GPS data point is similar to the InSAR data point weights, so the InSAR data dominate the optimization due to their large number. The effect of adding the ascending data is therefore larger, contributing significantly to the reduction of fault model parameter uncertainties represented by the histograms (Fig. 9) and the effective suppression of model parameter trade-offs, for example, between fault dip and strike-slip.
The fault width appears to be poorly constrained in all three cases (Fig. 9) and the histograms are bimodal, revealing the existence of two almost equally likely fault model solutions. The second favourable fault model, with a relatively large amount of slip on a small-width fault, can not be ruled out by adding the ascending In-SAR data. The large width of most the model parameter histograms demonstrates the importance of estimating fault model uncertainties, particularly when the observations are of a low SNR and are limited in their number and/or spatial extent.

Model parameter uncertainties using different data weights
The data weighting has a strong influence in the fault model estimation and we therefore separately test the impact of different data weights on the fault model optimizations and the influence of additional input data. We compare the result from the simulation of the fault models using data weights based on the variance matrix  only (Model 1, Model 2 and Model 3, Section 4.4) with simulations considering also the data covariances in the data weights (section 4.2). The distribution of fault parameter values resulting from the latter simulations reflect the uncertainties of fault models we call Model 1 * , Model 2 * and Model 3 * .
When the full data autocovariances are included the simulation results show improved precision in the estimates of all model parameters (Fig. 10), compared to the model estimations using only variance matrices. Most remarkable is the suppression of one of the bimodal histogram peaks for fault width, now showing a clear preference for a fault model with a width of more than 5 km (Fig. 10).
Beside the more precise fault model estimations, we can also detect a shift of the parameter histograms for Model 2 * and Model 3 * , with respect to Model 2 and Model 3, for example, for the slip-components (Fig. 10). We find that this shift in fault slip is accompanied by shifts in the estimates of the orbital ramps also and is likely caused by the reduced InSAR data weights near the centre of the interferograms. This shift is also noteworthy because fault slip parameters of most Model 3 * solutions fall outside the 95 per cent confidence interval of Model 3. This reveals that the model parameter uncertainties outlined by the Model 3 histograms do not reflect the true uncertainty.
The introduction of autocorrelated data errors has a pronounced positive effect on the precision of the model parameter estimates for all simulations (Fig. 10). The influence of GPS is clearly visible and stronger than in the variance-weighted simulations. The data point weights of the GPS measurements (Fig. 5) are larger than those of the adjacent InSAR points. We find that taking autocorrelated data errors into account is not only a more consistent way to balance independent data sets in fault model optimizations, but also leads to better results. We link this obvious improvement to the way the data are subsampled. The quadtree subsampling is an irregular subsampling approach that samples certain regions of an image more densely to capture details of the deformation signal. When we do not balance the InSAR data against each other, but allocate equal weights to them, areas with a high point density will be overrepresented in the model cost evaluation. Densely sampled areas are usually near the fault, but there is no physical reason to focus on near fault areas, while putting lower weights on the surrounding regions. We show that the covariance weighting circumvents this problem  ( Fig. 10), because the weight per area becomes more even. The 'far field' information is thus better represented in the optimization and the results become more stable. The good agreement of the 95 per cent 'confidence ellipses' drawn from the 2500 original optimizations (Fig. 8) with the simulation of Model 3 * validates our simplified one-fault optimizations and the limited number of 200 runs in each simulation. The fault model parameters estimated by Pagli et al. (2003) lie outside these confidence regions in most cases, especially the fault dip. From our analysis we therefore reject the Kleifarvatn fault model estimated by Pagli et al. (2003). The fault model estimated byÁrnadóttir et al. (2004) also differs significantly from our optimal model, particularly in fault length. But considering that we only take data errors into account, we stand back from rejecting it based on the results given in Fig. 8 and Table 2.

Variable slip and rake model
Assuming uniform slip on fault planes is a simplification since heterogeneity in fault displacement and rake is observed on all scales of faults (Mai & Beroza 2002). Also, systematic near-field residuals of the uniform-slip model (Figs 7c and d) suggest that a uniform-slip model is too simple. Therefore, we use the fault plane estimated in Section 4.2, extend its length to 10 km, its width to 8 km, and divide it into 10 × 8 subfaults to estimate variations in slip on the fault.
Taking the locations and orientations of the Núpshlí arháls fault and the Kleifarvatn subfaults, we solve for the slip, which constitutes a linear inverse problem: where G is the Green's function relating the data vector d to the model vector m, and R is again the weighting matrix. In our specific case the two orthogonal components of fault slip m are related to the observed surface displacement d through G, defining the fault plane and the medium. The model parameters defining the geometry of the fault plane carry considerable uncertainty, as we showed in Section 4.3. Therefore, we solve eq. (9) not only for the observed data d obs and the optimum plane parameters incorporated in the Green's function, G 0 , to obtain our optimum distributed slip model m 0 : but we also solve the problem for the N data sets modified with synthetic data errors ε synth,i and their corresponding optimum model plane (Fig. 8) applied in the Green's functions G i : The distribution of m i will reflect the sensitivity according to variations of the two fault planes and the data noise ε synth,i . We apply the Fast Non-Negative Least Square (FNNLS) algorithm of Bro & de Jong (1997) to the problem in eq. (10). The non-negativity of the algorithm allows the rake to vary only 90 de- grees and we rotate the local in-plane coordinate axes to allow rake variations (±45 • ) about the estimated rake of the uniform source model (189 • ). We introduce correlations between slip on neighbouring subfaults through a smoothing operator D (Jónsson et al. 2002), with the balance between data misfit and fault slip smoothness controlled by factor κ. The coupled system of equations becomes: We compile a cost-roughness trade-off curve of our specific case by inverting eq. (12) multiple times using different values of κ (Fig. 11). One can select the appropriate value of κ from where the misfit stops decreasing strongly with increasing slip roughness. In our case we find that high roughness causes solution instability, such that some near-zero slip values become negative. We chose κ = 2.2 (Fig. 11) because it corresponds to the start of the trade-off curve bend and ensures a high degree in solution stability for all N + 1 inversion problems (eqs 10 and 11).
The optimum source model m 0 has a concentration of fault slip of more than 60 cm in the upper central part of the fault, extending from the surface down to 3 km depth (Fig. 12). In the southern part of the fault the total displacement decreases from the surface to very small displacements below 5 km depth, while in the central and northern parts of the fault the fault slip decreases only slightly below 4 km. In summary, significant fault displacement of above 10 cm is restricted to the area of the simple rectangular fault model.
The uncertainties of m 0 are outlined by the set of solutions m i (Fig. 12). At depths between 1 and 4 km the fault slip shows the smallest confidence intervals, indicating that it is constrained to within about 20-30 cm. The slip constraints for the shallowest part of the fault are slightly poorer than at the centre of the fault plane. The lateral extent of the rupture is well determined through sharply decreasing slip with little scatter in the solutions m i . For the deeper and central part of the fault the slip only decreases to a value about 10-20 cm, while the relative uncertainty of the displacement increases strongly with depth. At these subfaults the scatter of m i clearly outlines the parameter bounds of the FNNLS algorithm along the axes of the local rotated coordinate system, so there we have very little fault slip resolution.

D I S C U S S I O N A N D C O N C L U S I O N S
With the use of additional InSAR data from an ascending track, together with descending InSAR and GPS data, we have located the previously unknown fault east of lake Kleifarvatn more precisely than the two earlier geodetic source studies by Pagli et al. (2003) and byÁrnadóttir et al. (2004). We have shown that the differences between the three source models Árnadóttir et al. 2004, this study) primarily arise from model parameter trade-offs. These trade-offs are effectively reduced when using InSAR data from the two viewing directions. We distinguish this improvement of using additional data from the impact of using different data weighting in the source model estimation, for which we also find considerable fault parameter estimation improvement when considering correlated data errors.
The improvement of fault model precision is not the only good argument for a weighting approach based on the full data covariance. This implementation of empirical data error characteristics is also a consistent and a quantitative way to combine independent data sets. We further note that we reduce the bias of the InSAR data subsampling with the strict propagation of the full data covariance matrix to the covariance matrix of the subsampled data sets. Both these effects increase the precision of the model parameter estimates and the reproducibility of the modelling results, which justifies the effort estimating of the data-error characteristics even more.
Aftershock locations can provide independent information about the location of the rupture plane. Relocated aftershocks of the Kleifarvatn earthquake are concentrated at depths between 5.5 and 7 km east of the lake and form a 6 km long, slightly north-northeast striking lineament, extending from the hypocentre region to the south (Hjaltadóttir & Vogfjör 2005, Fig. 13). Together with a few shallow aftershocks they outline a steeply (88 • E) dipping plane. The accuracy of the hypocentre locations depends on the quality of the velocity model. The velocity model used for the aftershock relocalization is of limited resolution leading to large uncertainties, especially in hypocentre depths, which tend to be overestimated (Hjaltadóttir & Vogfjör 2005). Therefore, the dip of the estimated fault plane from aftershocks (88 • E) is rather too large as well and we consider the difference between these results and our estimation about the fault dip ( Fig. 8 and Table 2) insignificant. The distributed slip model shows a concentration of shallow slip at the southeastern lake shore where surface ruptures were found in the field (Clifton et al. 2003, Figs 12 and 13). Some of these rupture traces form small pressure ridges, pointing to a localized component of thrust faulting in the shallow part of the fault, which is indicated in our geodetic source model (Figs 12  and 13). Furthermore, shallow fault slip extends beneath the lake to its northern shore near the epicentre (Fig. 13). Clifton et al. (2003) found here another concentration of surface ruptures and hypocentres of relocated aftershocks also show an increased activity in the north (Hjaltadóttir & Vogfjör 2005). The fault slip of our distributed slip model is somewhat small (0.3 m) in the northern The multicoloured stripe shows the trace of our model plane, where the colour represents the total slip on the shallow subfaults. The dashed circles outline the area of surface ruptures observed by Clifton et al. (2003). The dashed dark-blue line gives the outline of aftershock locations projected to the surface (Hjaltadóttir & Vogfjör 2005), and the star the epicentre of the Kleifarvatn event.
part and the uniform slip model does not extend to the northern lake shore. The source models by Pagli et al. (2003) andbyÁrnadóttir et al. (2004) have even smaller fault extensions and the match to observations in the field and from seismology is worse than for our new model. From our estimated parameter uncertainties we can reject the source model by Pagli et al. (2003) as being a realistic representation of the Kleifarvatn fault. The dominant discrepancies are the fault dip, the fault length and the slip on the fault (Fig. 8). The difference between our source model and the model byÁrnadóttir et al. (2004) is only significant for the fault length (Fig. 8). Since our estimation of the model parameter uncertainties is solely based on the estimated data errors, we stand back from rejecting this model.

A C K N O W L E D G M E N T S
We thank Steffen Knospe for sharing his knowledge on statistical analysis and for many fruitful discussions on the error estimation in InSAR data. We thank Martin Mai for reviewing an early version of the manuscript. We also thank ThóraÁrnadóttir for providing the GPS data and Kristín Vogfjör for seismological data of the Núpshlí arháls event and informative discussions. Comments from Gilles Peltzer and an anonymous reviewer improved the manuscript. The data were provided by ESA through Category-1 project #3639.