Journal of Volcanology and Geothermal Research 194 (2010) 139–149 Contents lists available at ScienceDirect Journal of Volcanology and Geothermal Research j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / j vo l g e o r e s Plate boundary deformation and man-made subsidence around geothermal fields on the Reykjanes Peninsula, Iceland M. Keiding a,⁎, T. Árnadóttir a, S. Jónsson b, J. Decriem a, A. Hooper c a b c Nordic Volcanological Centre, Institute of Earth Sciences, University of Iceland, Iceland Division of Physical Sciences and Engineering, King Abdullah University of Science and Technology (KAUST), Saudi Arabia Delft Institute of Earth Observation and Space Systems, Delft University of Technology, The Netherlands a r t i c l e i n f o Article history: Received 3 November 2009 Accepted 21 April 2010 Available online 3 June 2010 Keywords: plate boundary geothermal fields man-made subsidence triggered earthquakes Interferometric Synthetic Aperture Radar (InSAR) a b s t r a c t We present Interferometric Synthetic Aperture Radar (InSAR) data from 1992–1999 and 2003–2008 as well as GPS data from 2000–2009 for the active plate boundary on the Reykjanes Peninsula, southwest Iceland. The geodetic data reveal deformation mainly due to plate spreading, anthropogenic subsidence caused by geothermal fluid extraction and, possibly, increasing pressure in a geothermal system. Subsidence of around 10 cm is observed during the first 2 years of production at the Reykjanes geothermal power plant, which started operating in May 2006. We model the surface subsidence around the new power plant using point and ellipsoidal pressure sources in an elastic halfspace. Short-lived swarms of micro-earthquakes as well as aseismic fault movement are observed near the geothermal field following the start of production, possibly triggered by the stresses induced by geothermal fluid extraction. © 2010 Elsevier B.V. All rights reserved. 1. Introduction Many different natural and man-made processes associated with fluid migration at depth cause measurable deformation at the surface. The fluid-related processes are often so large that they locally obscure the deformation due to tectonic processes such as plate boundary deformation. Examples of processes involving fluid migration are ground-water extraction (e.g. Amelung et al., 1999; Hoffmann et al., 2001; Anderssohn et al., 2008), mining (e.g. Donnelly, 2009), geothermal or hydrocarbon production (Grasso and Wittlinger, 1990; Mossop and Segall, 1997; Fialko and Simons, 2000), naturally occurring fluctuations in geothermal and magmatic systems (Wicks et al., 1998; Peltier et al., 2009), or transient post-seismic processes (e.g. Jónsson et al., 2003). Probably the most prominent example of man-made subsidence around a geothermal reservoir is the Wairakei geothermal field in New Zealand, where 50 years of geothermal fluid extraction has resulted in a total of 15 m subsidence (Allis et al., 2009). The host rock deformation associated with geothermal fluid extraction can provide important insight in the extent, morphology and dynamics of the subsurface fluid reservoirs (e.g. Glowacka et al., 1999; Fialko and Simons, 2000; Vasco et al., 2002). The fluid flow in reservoirs is often highly anisotropic due to variations in permeability related to geological structures such as faults or sediment composition (Amelung et al., 1999), hence spatially dense ⁎ Corresponding author. E-mail address: keiding@gfz-potsdam.de (M. Keiding). 0377-0273/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.jvolgeores.2010.04.011 observations are needed in order to fully map the resulting ground deformation. InSAR offers excellent possibilities for this. Whereas ground-based observations, such as levelling and GPS data, are usually sparse, the radar technique can provide very dense spatial sampling of the ground deformation. In one example, (Fialko and Simons, 2000) examined InSAR data showing the subsidence around the Coso geothermal field in California, and modelled the subsidence using multiple ellipsoidal sources in an elastic halfspace. They also showed that clusters of micro-earthquakes associated with the geothermal fluid extraction may result from perturbations in the pore fluid pressure, as well as normal and shear stresses caused by the contraction of the geothermal reservoir. In this paper we examine the ground deformation on the Reykjanes Peninsula in southwest Iceland, using a combination of descending and ascending InSAR, as well as GPS data. The MidAtlantic plate boundary comes onshore on the Reykjanes Peninsula, where it forms a diffuse transtensional plate boundary zone characterised by high seismicity and recent volcanism (Fig. 1). The main tectonic features on the peninsula are a large number of NEtrending eruptive fissures and fractures, grouped into four volcanic fissure swarms (Sæmundsson, 1978; Clifton and Kattenhorn, 2006). The volcanic fissure swarms are intersected by a series of N–S oriented right-lateral strike-slip faults, which are the surface expressions of the left-lateral E–W shear at depth. Several high-temperature geothermal fields are present on the peninsula, located primarily at the intersections of the eruptive fissures and the strike-slip faults (Amy Clifton, personal communication, 2009). Following the start of geothermal energy production in the Reykjanes field in 2006, a 140 M. Keiding et al. / Journal of Volcanology and Geothermal Research 194 (2010) 139–149 Fig. 1. Tectonic map of the Reykjanes Peninsula, with fracture locations from (Clifton and Kattenhorn, 2006). The fractures are mainly NE–SW trending normal faults and tension fractures located within four volcanic fissure swarms. The hatched areas show the locations of high-temperature geothermal fields, labelled as R: Reykjanes, E: Eldvörp, S: Svartsengi, K: Krísuvík, B: Brennisteinsfjöll and H: Hengill. The Iceland inset shows the neovolcanic systems (grey shades) and the location of the study area. The arrows show the direction of the 2 cm/yr spreading across the peninsula between North America and Eurasia (DeMets et al., 1994). marked zone of subsidence of several cm/yr has evolved around the power plant. We examine the observed subsidence in some detail, to gain insight into the reservoir dynamics and the effect of the reservoir contraction on the surrounding crust. 2. Utilisation of geothermal energy on the Reykjanes Peninsula The utilisation of geothermal waters has been an integral part of people's life since the settlement of Iceland in the 9th century. The capital Reykjavík bears in its name a clear reference to geothermal springs (Reykjavík literally means “Smoky Bay”), and historical records describe how the springs were used for washing and bathing in past centuries. As of 2009, geothermal energy accounts for around 25% of the electricity production and almost all domestic heating in Iceland. Four geothermal power plants are in operation on the Reykjanes Peninsula and in the Hengill area. In 1976, the Svartsengi geothermal power plant was taken into use, and it has been progressively expanded since then. The most recent expansion of Svartsengi took place in early 2008, resulting in a production capacity of 75 MW electricity and 150 MW heat. The Nesjavellir power plant started operating in the Hengill area in 1990, and is today the largest geothermal power plant in Iceland with a capacity of 120 MW electricity and 300 MW heat. In May 2006, production started in the Reykjanes power plant on the tip of the Reykjanes Peninsula with a capacity of 100 MW electricity. Later that year, the Hellisheidi power plant started operating in the southern Hengill area. The Hellisheidi power plant holds a production capacity of 210 MW electricity, as of February 2009, but there are plans to expand its capacity to a total of 300 MW electricity and 400 MW heat. As the pressure drawdown may diminish the well field productivity, waste fluids are typically reinjected into the geothermal reservoirs. In Svartsengi, reinjection has been carried out at intermittent rates since 1984, but injection was increased progressively during 2002–2008 so that around 50% of the volume of extracted water was reinjected in 2008 (Vatnaskil, 2009). Systematic reinjection in the Reykjanes field has not started, as of summer 2009. In Hellisheidi, all waste fluids from the production are reinjected into the reservoir. The reinjected waste fluids, however, never make up the volume of the extracted fluids due to the loss of steam to the air. Therefore, a pressure decrease occurs, and results in contraction of the rock matrix within the reservoir, which in turn causes subsidence above the reservoir. Subsidence in the Svartsengi field was first documented by a levelling and gravity study based on repeated measurements during 1975–1999 (Eysteinsson, 2000). The results of the levelling showed subsidence rates between 7 and 14 mm/yr, with the highest rates during the first years of production. The study also demonstrated that the subsidence at Svartsengi varies linearly with the pressure decrease observed at 900 m depth in boreholes. The subsidence around Svartsengi was later confirmed by an InSAR study (Vadon and Sigmundsson, 1997), as well as GPS studies (Hreinsdóttir et al., 2001; Magnússon and Thorbergsson, 2004; Árnadóttir et al., 2006; Keiding et al., 2008). 3. Data and methods 3.1. GPS data analysis We report GPS data from a network of around 60 campaign stations and 8 continuous stations on the Reykjanes Peninsula and the Hengill area. Annual surveys of selected campaign stations have been carried out since 2000. Each campaign measurement lasted at least two days during 2000–2006 and three days during 2007–2009. The GPS data analysis was done in two steps. First, we calculate daily solutions using the Bernese v5.0 software (Dach et al., 2007), with orbit and Earth rotational information from the International GPS Service (IGS) (Dow et al., 2005). Six international IGS stations were included in the processing to aid the stabilisation in the International Terrestrial Reference Frame (ITRF). Second, we combine the daily campaign and Icelandic CGPS solutions with three IGS global solutions (IGS1, IGS3, and EURA) using the GLOBK software (Herring et al., 2006). We use a regional stabilisation approach (McClusky et al., 2000) to estimate the station positions and velocities in a reference frame that is approximately aligned with the ITRF2005 (Altamimi et al., 2007). Discontinuities in the time series are caused by earthquakes in 2003 and 2008. On 23 August 2003 a MW5.0 earthquake occurred on the central Reykjanes Peninsula, causing horizontal offsets of up to 1 cm at eight stations (Keiding et al., 2008). The velocities at these stations are computed using only the GPS data after the 2003 earthquake. Two main shocks with a total moment release equivalent of MW6.2 occurred immediately east of the study area on 29 May 2008, causing coseismic offsets of up to several cm on the eastern part of the peninsula (Hreinsdóttir et al., 2009; Decriem et al., 2010). Hence the M. Keiding et al. / Journal of Volcanology and Geothermal Research 194 (2010) 139–149 velocities at all stations east of 22°W are estimated using data before the 2008 earthquake sequence. Furthermore, the start of production in the Reykjanes and Hellisheidi geothermal fields in 2006 affected the velocities at nearby stations, as seen in the time series of the campaign stations RNES and HH04 in Fig. 2. We therefore divide the time series at affected stations into two time periods before and after 2006, using the GPS campaign measurements in March–April 2006 as respectively the last and the first observation of each. In the following, we mainly focus on these two time periods, before and after the start of production in the Reykjanes and Hellisheidi geothermal fields. However, for stations that show steady velocities we choose to include all available data during 2001–2009, in order to obtain the best constraints on the linear regression estimates. On average, five observations are used for the regression of constant velocities at the campaign stations. For the continuous stations we compute the velocities using observations at the time of each campaign. 3.2. InSAR data analysis The Reykjanes Peninsula is well suited for a radar based study because its surface mainly consists of young and sparsely vegetated lava fields, hence the surface reflectivity is sufficiently high and changes little with time. We form InSAR images from data collected by the ERS and Envisat satellites, operated by the European Space Agency. The ERS and Envisat data sets are from the descending track 138, and comprise 18 images spanning 11 May 1992–16 Sep 1999 and 22 images spanning 25 Sep 2003–18 Sep 2008, respectively. All images were acquired during May–October to avoid decorrelation due to snow. The SAR system operates with a side-looking geometry, illuminating an approximately 100 km wide swath of the ground. The incidence angle at ground is 20–26°, and the descending line-of-sight (LOS) unit vector from ground to satellite is approximately [east north up] = [0.4 −0.1 0.9], hence the observations are most sensitive to vertical ground motion, less sensitive to east–west motion and least sensitive to north–south motion. The SAR image has ground resolution elements, or pixels, of approximately 5 × 20 m. The data are processed with the StaMPS/MTI software (Hooper et al., 2007; Hooper, 2008), which applies both a Persistent Scatterer (PS) and a Small Baseline (SB) approach. This multi-temporal InSAR method involves the processing of multiple acquisitions and addresses the problems of decorrelation caused by differences in position and orientation of the master and slave sensor or by physical 141 changes at the surface (e.g. Zebker and Villasenor, 1992). The PS method identifies pixels that are dominated by the echo from a single bright scatterer and therefore have little decorrelation due to changes in satellite geometry or relative movement of scatterers within the pixel. The SB method, on the other hand, minimises the decorrelation by computing interferograms with short temporal and spatial baselines. The combination of the two methods has the potential for improving the spatial sampling considerably, and thereby increase the resolution of deformation signals and aid a more reliable phase unwrapping. For the PS processing a master is chosen that minimises the sum of decorrelation due to the time interval, the perpendicular baseline and the difference in Doppler frequency. For the SB processing we form 42 ERS and 63 Envisat multiple-master interferograms. The Doris software (Kampes et al., 2003) is used for the interferometric processing. Each slave image is resampled to the master geometry and corrected for the difference in position of the master and slave sensor, using the WGS84 reference ellipsoid and a 25 m digital elevation model from the National Land Survey of Iceland. The pixels selected by the two methods are then combined and the phase unwrapped using a statistical cost flow algorithm applicable to singleor multiple-master time series (Hooper et al., 2009). Finally, the unwrapped phase is corrected for atmospheric delay plus errors in orbits and the elevation model, using a combination of temporal and spatial filtering (Hooper et al., 2007). 4. Results 4.1. GPS velocities Fig. 3 shows the GPS velocities relative to stable North America, computed using the ITRF2005 absolute rotation pole for the North American plate (Altamimi et al., 2007). The horizontal GPS velocities on the Reykjanes Peninsula mainly reflect the plate motion, that is, leftlateral shear in the E–W direction as well as some N–S extension. The station velocities are close to zero on the northern part of the peninsula, and gradually increase in magnitude moving south across the plate boundary zone. The stations along the southern shore of the peninsula are moving toward ESE with horizontal rates of 18–19 mm/yr relative to Reykjavík, close to the full spreading rate between North America and Eurasia (DeMets et al., 1994). The vertical GPS velocities are close to zero on the central part of the peninsula, but uplift is observed along the SE part of the peninsula. During 2001–2006, local subsidence is observed in the Svartsengi field, with a subsidence rate of 5–10 mm/yr, while no clear signal of subsidence is observed around the Nesjavellir power plant in the northern Hengill area. During 2006–2009, marked signals of subsidence appear around the Reykjanes and Hellisheidi fields. The Reykjanes subsidence is confined to a small area on the tip of the peninsula, with a maximum subsidence rate of around 40 mm/yr at the nearest campaign GPS stations RNES (see Fig. 2). The Hellisheidi subsidence bowl covers a larger area, in accordance with the larger extent of the well field. The maximum subsidence rates are smaller than in the Reykjanes field, however, the rates are not well-constrained due to the short time span of these station velocities from 2006 to 2008 (until the May 2008 earthquake). Interestingly, the station DRAU which is located less than 1 km from the Hellisheidi power plant, has close to zero vertical rate, probably because it is located near one of the areas where waste fluids are being reinjected (see Fig. 3). 4.2. InSAR time series Fig. 2. North, east and vertical time series at the campaign stations RNES and HH04 (see station locations in Fig. 3), in the ITRF2005 reference frame. Error bars indicate 1σ uncertainties. The solid vertical lines show the approximate start of production at the Reykjanes and Hellisheidi geothermal power plants. The stippled line shows the time of the earthquake sequence on 29 May 2008. The analysis of the InSAR data results in ERS and Envisat time series that each comprises around 2 million pixels, providing an exceptionally good spatial sampling of the ground deformation on the Reykjanes Peninsula. Small areas of decorrelation are observed at 142 M. Keiding et al. / Journal of Volcanology and Geothermal Research 194 (2010) 139–149 Fig. 3. GPS velocities based on annual measurements during 2001–2009. The arrows show horizontal velocities with 1σ confidence ellipses, while the contour colours in the background show interpolated vertical velocities. Black dots show locations of production boreholes at the Svartsengi (S), Nesjavellir (N), Hellisheidi (H) and Reykjanes (R) geothermal power plants, while red dots show the location of reinjection boreholes in Hellisheidi. The velocities at stations near the Reykjanes and Hellisheidi power plants have been computed for 2001–2006 (upper panel) and 2006–2009 (lower panel, western peninsula) or 2006–2008 (lower panel, eastern peninsula). Note that only the velocity vectors coloured in red are from time series that have been split up into two periods, the other velocities are based on all available data. Some variations in the velocity vectors between the two time periods are due to changes in the network configuration, that is, abandonment of old stations or installation of new ones. steep slopes and in areas with lakes or relatively dense vegetation. The InSAR time series show relatively steady regional deformation, but the rates and spatial extent of deformation vary locally around the geothermal fields. The main change in the deformation rates occurs in 2006, due to the start of production in the Reykjanes and Hellisheidi fields. Thus we divide the InSAR time series into three intervals, and compute the mean LOS velocity fields for the time periods 11 May 1992–16 Sep 1999 (ERS), 25 Sep 2003–29 Sep 2005 (Envisat) and 6 Jul 2006–1 May 2008 (Envisat). Several of the ERS and Envisat interferograms have a bilinear phase ramp, showing increasing range change from east to west. The ramps may be partly due to orbital errors, but they may also reflect a regional subsidence of the western part of the peninsula, relative to its eastern part, as has been documented by a countrywide GPS study (Árnadóttir et al., 2009). In order to better display the local signals of deformation along the Reykjanes Peninsula, we estimate a bilinear ramp for each image by linear regression and remove the ramps before computing the mean LOS velocities. The 2003–2005 and 2006–2008 LOS velocity fields have some noise, seen as areas of patchy LOS rates, reflecting that they are estimated from only nine images for each period. The standard deviation of the mean LOS velocity for individual pixels is around 1 mm/yr for the 1992–1999 rates and 3–7 mm/yr for the 2003–2005 and 2006–2008 rates. Although the standard deviation on individual pixels is sometimes high, the pattern of deformation is generally smooth and the signal-to-noise ratio is improved by multilooking the radar data. Fig. 4. Residual mean LOS velocity fields after removal of bilinear ramps, relative to the mean value during each period in the area near Reykjavík (shown with the box). The time spans of the images are 11 May 1992–16 Sep 1999 (18 ERS images), 25 Sep 2003– 29 Sep 2005 (9 Envisat images) and 6 Jul 2006–1 May 2008 (9 Envisat images). The profiles in the bottom panel show moving averages of the mean LOS rates and 1σ standard deviations along the line shown on the maps. The LOS velocity fields in Fig. 4 reveal both regional and local deformation. Positive LOS rates indicate motion toward the satellite (primarily uplift), while negative LOS rates show motion away from the satellite (primarily subsidence). All three images show increasing LOS rates moving from north to south across the peninsula, as illustrated with the profiles in the lower panel of Fig. 4. This increase in LOS rates must in part reflect the increase in eastward velocities across the plate boundary zone, but it is also possible that they reflect some uplift along the SE part of the peninsula, as indicated by the GPS velocities in Fig. 3. A subtle zone of negative LOS rates are observed along the central part of the peninsula during 2003–2005 and 2006– 2008, with negative rates of 0–4 mm/yr relative to Reykjavík. The negative rates most likely reflect subsidence, caused by the extension 143 M. Keiding et al. / Journal of Volcanology and Geothermal Research 194 (2010) 139–149 across the plate boundary. The subsidence may indicate that the extension is not completely balanced by inflow of material from below, as suggested by (Vadon and Sigmundsson, 1997). During 1992–1999 a marked zone of positive LOS rates is observed along the eastern margin of the image. A persistent earthquake swarm and uplift of the Hengill area was observed during 1994–1998, culminating in 1998 with two earthquake swarms including moderate sized events on the eastern part of the peninsula (Sigmundsson et al., 1997; Feigl et al., 2000; Clifton et al., 2002). The positive LOS rates during 1992–1999 thus reflect uplift and possibly coseismic deformation or some widening of the Hengill fissure swarm, as suggested by a GPS study (Keiding et al., 2008). The LOS velocity fields also show local deformation around the geothermal fields. Negative LOS rates are observed around the Svartsengi field in all three images. During 1992–1999, the maximum negative rates average 5 mm/yr, relative to Reykjavík, but a considerable higher negative rate of 20 mm/yr was in fact observed during 1992–1993, decreasing to 4 mm/yr after 1993. The varying subsidence rate is most likely due to changes in the reinjection of waste fluids at Svartsengi, as reinjection was taken up again in 1993, after a break during 1991–1992 (Vatnaskil, 2009). The subsidence around Svartsengi is elongated in the NE–SW direction, and includes the Eldvörp geothermal field located 5 km further SW (see Fig. 1). The Eldvörp field has not yet been directly utilised, but a pressure connection between the Svartsengi and Eldvörp fields indicates that fluids are also withdrawn from Eldvörp during production in Svartsengi (Eysteinsson, 2000). During 2003–2005 and 2006–2008, negative rates of around 10 mm/yr are observed in the Svartsengi field, but the subsidence signal is less localised than during 1992– 1999, and part of the negative rates seems to be related to the regional zone of subsidence along the central plate boundary zone. After 2006, a marked signal of negative LOS rates appears around the Reykjanes geothermal field, reflecting the newly formed subsidence bowl due to geothermal fluid extraction. The subsidence signal in the Reykjanes field shows up almost immediately after the power plant was put into operation in May 2006, and the maximum negative LOS rate during 2006–2008 is more than 30 mm/yr. Negative LOS rates are also observed around the Hellisheidi geothermal field in the east. Finally, an area of positive LOS rates is observed in the Krísuvík area on the central part of the peninsula. Local anomalies are also observed at some of the GPS stations in this area, which show accelerating uplift and SE motion during 2008–2009. The most likely explanation for this signal is that the geothermal system in the Krísuvík area is inflating due to overpressure. Another possibility, however, is that we are seeing the signs of a slow magmatic intrusion (Keiding et al., 2009a). The ascending and descending radar LOS displacements are shown in Fig. 5a and b, with displacements relative to the mean values in the NW part of the area (as shown with the box in Fig. 5a). The ascending and descending displacements show a similar deformation pattern, mainly the subsidence bowls around the Reykjanes and Svartsengi geothermal fields. The ascending radar data show a larger maximum negative displacement (∼12 cm) than the descending radar data (∼8 cm), but this is probably due to the horizontal shear along the plate boundary rather than a discrepancy in the vertical displacements, as described below. The average LOS unit vectors for the ascending track 173 and descending track 138 in the area shown in Fig. 5 are ([east north up]) nasc ndes = = ½−0:32 −0:10 0:94Š ½0:40 −0:11 0:91Š ð1Þ Adding and subtracting the mean LOS displacements gives the following sensitivity vectors nþ = nasc + ndes n− = ndes −nasc = = ½0:09 ½0:72 −0:21 1:85Š −0:01 −0:04Š ð2Þ Hence, adding the ascending and descending radar data provides near-vertical deformation and subtracting the ascending from the descending data will show deformation approximately in the east– 4.3. The western Reykjanes Peninsula We examine the subsidence around the Reykjanes geothermal field on the western Reykjanes Peninsula in more detail, using a combination of the GPS and InSAR data. The GPS data show the 3-dimensional deformation but are spatially sparse, while the InSAR data provide a dense spatial sampling but only a 1-dimensional observation of the ground deformation. Ascending radar data can add another 1-dimensional observation, because the ascending LOS unit vector of approximately [east north up]= [−0.4 −0.1 0.9] differs from the descending LOS unit vector. Too few ascending radar data are available to be processed using the multi-temporal StaMPS/MTI method. We therefore process ascending ERS data from track 173, using the GAMMA software (Werner et al., 2000), and include one such interferogram, spanning June 2005–May 2008. We compute the descending LOS displacement during same time period, from 13 images. The rate of deformation during this time period is not constant as it includes a full year before the start of production in the Reykjanes field, thus we show the total displacement during this time period rather than annual rates. Fig. 5. Linear combinations of ascending and descending radar data and comparison with GPS data. a) Ascending LOS displacements during 18 June 2005–3 May 2008 from two Envisat track 173 acquisitions. b) Descending LOS displacements during 16 June 2005–1 May 2008, estimated from 13 Envisat track 138 acquisitions. c) Near-vertical radar displacements from addition of ascending and descending LOS displacements. d) Approximately east–west radar displacements obtained by subtracting the ascending from the descending LOS displacements. The coloured circles show the magnitudes of the GPS displacements projected onto the direction of the radar displacements. The radar data in all four panels are shown relative to the mean value within the area shown with the box in panel a, and the GPS data are shown relative to the continuous GPS station located within the box. The arrows in panels a and b show the line-of-sight direction from ground toward the ascending and descending satellites. 144 M. Keiding et al. / Journal of Volcanology and Geothermal Research 194 (2010) 139–149 west direction. After forming the linear combinations of the InSAR data, we normalise the results to get the displacements in cm. The near-vertical radar displacements in Fig. 5c indicate that the maximum subsidence in the Reykjanes geothermal field is around 10 cm, during June 2005–May 2008. The Reykjanes subsidence is clearly elongated in the NE–SW direction. The elongation thus aligns with the fractures in the area, suggesting that the permeability in the reservoir is controlled by the fractures. The largest subsidence is observed within an ellipse of approximately 3 × 5 km, but there is a NNE-ward extension to the subsidence bowl, making it slightly curved. The boundary toward NW is particularly sharp, but the boundary toward east is also quite sharp, indicating that there is little pressure connection between the Reykjanes field and the adjacent Eldvörp field. A comparison of the near-vertical displacements with the ascending and descending LOS displacements (Fig. 5a and b) shows that the maximum LOS displacements are offset by 1–2 km from the centre of subsidence in the near-vertical deformation field, because the horizontal motion is also mapped onto the LOS directions. The east–west radar displacements in Fig. 5d show horizontal motion on the order of several cm toward the centre of subsidence around the Reykjanes field. Furthermore, there is a clear regional pattern showing westward displacements (negative rates) on the northern part of the peninsula and eastward displacements (positive rates) on the southern part, reflecting the left-lateral shear along the plate boundary. In order to compare the radar and GPS data, we compute the GPS velocities assuming linear rates during 2005–2009, scale them to the time period of the radar data to get the displacements, and project the displacements onto the radar unit vectors. The values of the projected GPS displacements are shown with the coloured circles on top of the radar displacements in Fig. 5, and generally agree well with the radar data. A quantitative comparison can be obtained if we estimate the values of the radar data as the mean and standard deviation of the pixel values within a small area centred at each GPS station (here we use circular areas with 600 m diameter). Doing so, we find that the differences between the GPS and radar data are for the most part less than 1 cm and within the 1σ uncertainties of the data. around active volcanos (e.g. Battaglia and Hill, 2009) or geothermal fields (Fialko and Simons, 2000), than does the simple point source. The ascending and descending radar LOS displacements and the GPS displacements (in east, north and vertical) are included as three separate data sets in the optimisations. The data reflect both the local deformation due to production in the Reykjanes and Svartsengi fields and the plate motion, as demonstrated in Fig. 5d and described above. Thus, in order to model the subsidence around the geothermal fields, we first correct for deformation due to plate motion in both GPS and radar data. For the GPS data, we assume that the plate motion is constant during 2001–2006, and subtract this background signal from the velocities estimated for 2005–2009. The four stations located nearest to Svartsengi are affected by the subsidence around the geothermal field during 2001–2009, hence we correct for the plate motion signal at these stations using the average velocities at nearby GPS stations. We finally scale the velocities during 2005–2009 to match the time period of the radar data. In the radar displacement fields, the east–west shear along the plate boundary causes linear ramps perpendicular to the plate boundary. We estimate such ramps in areas where there is negligible subsidence (here we use data east of 22.35°W or north of 63.96°N), and subtract them from the LOS displacements. The InSAR displacements are calculated relative to an area in the northwestern part of the peninsula (shown with a box in Fig. 5a). In the optimisation, we therefore estimate a constant offset for the ascending and descending radar displacements to account for possible biases due to this choice of reference area. We sub-sample the radar data using a quadtree algorithm that reduces the number of data points while maintaining a good spatial representation of the deformation (e.g. Jónsson et al., 2002; Decriem et al., 2010). We start by dividing the radar displacement fields into squares of 2 × 2 km, which corresponds approximately to the correlation distance within the data. The squares that have a variance larger than a certain threshold are then recursively subdivided into quadrants until the variance within each quadrant does not exceed the threshold. The weights of the sub-sampled radar data are estimated from the full data variance–covariance matrix, as described by (Sudhaus and Jónsson, 2008). The sample covariogram for a discrete distance class hc is given by 5. Modelling the subsidence around the Reykjanes field 1 N ∑ dðr Þ⋅dðsi Þ 2N i = 1 i 5.1. Methodology Ĉðhc Þ = We estimate source models that may describe the subsidence around the Reykjanes and Svartsengi fields by joint optimisation of the GPS and InSAR data. The simplest source model relating ground deformation to volume change at depth is an isotropic point pressure source (Mogi, 1958), defined by four parameters describing its location (latitude, longitude, and depth) and volume change. The point source gives a good approximation to roughly equi-dimensional bodies undergoing uniform volume change and has been widely used to model observed deformation of active volcanoes (e.g. Lu et al., 2002). The point source has also been applied to model subsidence due to geothermal fluid extraction (Mossop and Segall, 1997; Fialko and Simons, 2000), although more than a single point pressure source is usually required in order to fit the spatial irregularities observed around geothermal fields. Another commonly used elastic halfspace source model that relates ground deformation to pressure change at depth is the finite prolate spheroidal source derived by (Yang et al., 1988). The ellipsoidal source is defined by eight parameters describing its location (latitude, longitude, and depth), geometry (length of semimajor and semiminor axes), orientation (strike and plunge of the semimajor axis) and uniform pressure at the ellipsoidal surface. The ellipsoidal source with its slightly more complicated finite geometry, often provides a better fit to the deformation observed where d(ri) and d(si) are the radar LOS displacements at the locations ri and si separated by distances of approximately hc, and N is the number of data pairs in the distance class. We estimate the sample covariograms in the northwestern part of the peninsula (west of 22.5°W and north of 63.95°N), where there is negligible deformation from the geothermal fields and plate boundary shear. The sample covariograms are fitted by exponential decay functions, which are then used to compute the covariance matrix for the sub-sampled data sets. The two-pass ascending interferogram has more noise and thus a covariance that is approximately three times higher than the covariance of the averaged descending interferogram. As a result, the ascending data are given less weight in the optimisation. The weights of the GPS data are simply based on the variance of the estimated velocities. After sub-sampling, we have 185 ascending and 274 descending data points, while the GPS data set includes 16 displacement vectors in east, north and vertical. We apply a non-linear optimisation scheme to find the set of model parameters that minimises the weighted residual sum of squares, WRSS = rT∑− 1 r, where r is the difference between the observed and predicted displacements, and ∑ is the data covariance matrix. During optimisation we use a simulated annealing algorithm, followed by a derivative-based algorithm (e.g. Cervelli et al., 2001). The simulated annealing performs a random ð3Þ 145 M. Keiding et al. / Journal of Volcanology and Geothermal Research 194 (2010) 139–149 search through a predefined model parameter space and identifies the region near the absolute model-cost minimum. The algorithm is able to escape local minima due to the randomness of the search, but it does not necessarily find the minimum itself. The result of the simulated annealing is therefore passed to a derivative-based algorithm to improve the model further. The mean and confidence intervals of the model parameters are estimated using a bootstrap algorithm (Efron and Tibshirani, 1986), that performs the optimisations on a large number of randomly resampled data sets and computes the confidence intervals from the range of the estimated parameters. We report the goodness of fit to the data using the WRSS , where N is the number reduced chi-squared, calculated as χ2ν = N−m of data and m is the number of unknown model parameters. Table 2 Estimated ellipsoidal source parameters for the subsidence around the Reykjanes field. A Poisson's ratio of 0.25 and a shear modulus of 10 GPa are used (Fialko and Simons, 2000). The volume decrease is computed using the relation for cigar-shaped ellipsoidal sources ΔV = VΔP/μ, where V is the ellipsoid volume, ΔP is pressure change and μ is the shear modulus (Amoruso and Crescentini, 2009). 5.2. Results of modelling 6. Discussion The simplest model with a single point source for each of the Reykjanes and Svartsengi subsidence bowls has a χ2ν = 4.2. The fit is improved if we increase the number of point sources in the Reykjanes field to two (χ2ν = 3.7) or three (χ2ν = 1.8). In the two- and three-point source models, the sources align along the NE–SW trend of the subsidence signal, at 1.5–2.8 km depth (Tables 1). We then test a model with an ellipsoidal source for the Reykjanes subsidence and a point source for the Svartsengi subsidence. When the optimisation is performed with unconstrained strike and plunge and loose constrains on the other parameters, the resulting ellipsoidal parameters show a bimodal distribution: one group of the estimated ellipsoids are centred at 2 km depth and plunge shallowly toward NE, while another group of ellipsoids are deeper (∼ 4 km) and have steep plunge (∼70°). A deep and steep source seems unrealistic, given that the depth of production boreholes in the Reykjanes field is only 1–2.5 km (Jónsson et al., 2009). A sensitivity study by (Shirzaei and Walter, 2009) has, in fact, demonstrated that the plunge and downward depth of the ellipsoidal source may not be well-constrained. If we constrain the plunge to be less than 45°, we get consistent results with a χ2ν = 3.0. The results of the bootstrap optimisations indicate an ellipsoidal source plunging 10° toward N53°E and centred at 2.2 km depth. The ellipsoid is clearly elongated with a semimajor to semiminor axis ratio of 5 (Tables 2). Fig. 6 shows the data that were used as input in the optimisations and the predicted displacements from the ellipsoid model. The GPS and radar observations are reproduced to the first order. The residuals of the ascending data are generally larger than for the descending data, and some of the GPS stations show large residuals. This reflects that the descending data are given most weight during the optimisations due to their large number of data points and relatively low covariance. We note that the bootstrap confidence limits in Tables 1 and 2 seem small and may be underestimated due to the correlations in the InSAR data. The Svartsengi point source is located at 3.2 km depth and has a volume decrease of around 1.0 × 10− 3 km3. As for the Reykjanes subsidence, the single point source is a little deeper than the depth of the production bore holes, which in Svartsengi are located at 1–2 km 6.1. Subsidence and pressure changes Table 1 Estimated point source parameters for the subsidence around the Reykjanes field. No. is the number of point sources in the Reykjanes field. The confidence limits are 68% percentiles from the bootstrap model parameters. No. Lon (°W) Lat (°N) Depth (km) Volume decrease (×10− 3 km3) 1 2 22.670 22.697 22.645 22.711 22.665 22.615 63.823 63.818 63.836 63.816 63.826 63.855 3.4 ± 0.4 0.6 1.9+ − 0.5 2.6 ± 0.3 1.5 ± 0.3 2.0 ± 0.2 0.2 2.8+ − 0.3 0.8 3.6+ − 0.7 1.0 1.0+ − 0.6 1.6 ± 0.4 0.5 ± 0.2 0.3 1.2+ − 0.2 1.2 ± 0.2 3 Lon (°W) Lat (°N) Depth (km) Strike Plunge Semimajor Semiminor Volume decrease (°) (°) (km) (km) (×10− 3 km3) 22.657 63.832 2.2 ± 0.2 53 ± 1 10 ± 1 5.5 ± 0.5 0.5 0.9+ − 0.4 2.1 ± 0.2 (Jónsson et al., 2009). We also test the ellipsoidal source for Svartsengi, but the results are scattered, with a roughly equi-dimensional source without a preferred orientation. In the previous section we found that the observed surface subsidence around the Reykjanes geothermal field can be fitted well using point sources or a finite ellipsoidal source in an elastic halfspace. While the three-point source model provides the best fit to the data, the model with an ellipsoidal source is probably more physical, since it mimics the reservoir as a finite volume within the crust. Our elastic halfspace model does not consider the poroelastic processes related to the flow of the interstitial fluid and the deformation of the porous rock. However, the pore-pressure changes induced by fluid extraction diminish outside the reservoir, thus the deformation of the crust surrounding the reservoir can be assumed elastic (Segall and Fitzgerald, 1998). Our models, therefore, provide a reasonable simulation of the deformation of the crust surrounding the reservoir. The estimated volume changes from our ellipsoid and multiplepoint source models are in the range (2–3) × 10− 3 km3. During 2005– 2008, 57.9 Mton of fluids were extracted from the Reykjanes field (Vatnaskil, 2009). Assuming a density of 820 kg/m3 of the geothermal water (Eysteinsson, 2000), this corresponds to a volume withdrawal of 70.6 × 10− 3 km3, that is, ∼ 25 times the estimated volume change from our point and ellipsoid models. Although our models with simplified source geometries can only provide rough estimates of the actual volume decrease in the reservoir, this difference still seems large enough to indicate that it may be significant. A similar difference was estimated in the Svartsengi field during1975–1999 (Eysteinsson, 2000). The extracted fluids are mostly replaced by the natural recharge into the system. However, the natural recharge does not completely make up for the extracted fluids, hence the pore pressure within the reservoir decreases and the ground-water level drops. The decrease in pore pressure, in turn, results in a small contraction of the rock matrix. The apparent difference between the volume of the extracted fluids and the estimated volume change of the host rock indicates that the rock matrix is relatively strong and that the permeability is high and not considerably reduced by the fluid extraction from the reservoir. Fig. 7 shows a comparison of the vertical GPS displacements at the station RNES, the maximum descending LOS displacements in the Reykjanes field and the pressure observed at 1500 m depth in three boreholes (Jónsson et al., 2009). The boreholes RN12 and RN23 are production holes located within in the main well field, while RN16 is located some hundred meters further NW and is only used for monitoring and research. The geodetic data and the pressure observations show the same pattern: slow subsidence and slowly decreasing pressure during 2003–2006, followed by an abrupt change to higher rates at the time of the start of production in the Reykjanes field in 2006. The pressure observations clearly show that the pressure decrease starts to tail out during 2007, indicating that it is stabilising at a slower rate. The stabilisation of the pressure decrease is expected as the recharge into the system usually increases during the 146 M. Keiding et al. / Journal of Volcanology and Geothermal Research 194 (2010) 139–149 Fig. 6. Ascending and descending radar LOS displacements and east, north and vertical GPS displacements included in optimisation for the Reykjanes and Svartsengi source parameters. The predicted surface displacements from the model with an ellipsoidal source for Reykjanes and a point source for Svartsengi are shown, as well as the residuals between the observations and predictions. The coloured circles in the GPS figures show the vertical displacements. The colour scale for the ascending LOS, descending LOS and vertical GPS displacements is identical for all figures. first years of production, and the observed pressure changes are in good agreement with numerical simulations of the reservoir fluid and heat flow (Björnsson et al., 2008). The surface subsidence observed by GPS and InSAR also appears to tail off, albeit at a slower rate than the pressure observations, indicating that the crustal response is somewhat delayed from the pressure changes. The pressure decrease in RN16 is smaller than that observed in the production boreholes located within the main well field, however it shows a very similar pattern. The total pressure decrease in the production boreholes is around 3 MPa. 6.2. Induced seismicity Fig. 7. Time series spanning 2003–2008 with vertical GPS and descending radar LOS displacements, as well as pressure observed at 1500 m depth in three boreholes in the Reykjanes field. The GPS and radar displacements are in the same reference frame as in Fig. 5. Extraction of fluids from geothermal and hydrocarbon reservoirs are known to trigger seismicity in several cases, such as in the Geysers geothermal field (Eberhart-Phillips and Oppenheimer, 1984), the Coso geothermal field in USA (Fialko and Simons, 2000) and the Lacq gas field in southwestern France (Grasso and Wittlinger, 1990). Near M. Keiding et al. / Journal of Volcanology and Geothermal Research 194 (2010) 139–149 the Lacq gas field, earthquakes with magnitudes up to 4.2 have been recorded, in an area that is clearly separated from the tectonic seismicity in the Pyrenees (Grasso and Wittlinger, 1990). The seismicity began when the pressure had decreased by 30 MPa, 10 years after the start of gas extraction. The subsidence around the Lacq gas field varies linearly with the average reservoir pressure drop, indicating that the contraction of the reservoir can be explained by a linear poroelastic model (Segall et al., 1994). The poroelastic stressing due to reservoir contraction has been examined in detail by (Segall, 1989) and (Segall and Fitzgerald, 1998). They found through analytical modelling of a circular disk-shaped reservoir, that the least compressive horizontal stress within the reservoir decreases with decreasing pore pressures, thus enhancing tensional fracturing. Outside the reservoir production does not directly decrease the pore pressure so that the tendency for normal faulting is considerably larger. In extensional environments, normal faulting will be promoted near the edge of the reservoir, or anywhere there is a steep gradient in pore-pressure reduction (Segall and Fitzgerald, 1998). A change in the pattern of seismicity is observed following the start of production at the Reykjanes power plant (Fig. 8). Earthquakes on the Reykjanes Peninsula usually occur along a narrow ENE trending zone extending from the tip of the peninsula to the Hengill area (Tryggvason, 1973). An intense swarm of several thousand micro-earthquakes were recorded on the tip of the peninsula in 1972 by a local seismic network (Klein et al., 1977). Accurate epicentre locations of a subset of the swarm events showed that most of the earthquakes occurred within a less than 2 km wide and approximately 12 km long zone, as outlined in Fig. 8. Since the early 1990s, earthquakes in Iceland have been recorded by the SIL seismic network operated by the Icelandic Meteorological Office (Bödvarsson et al., 1999; Jakobsdóttir, 2008). Earthquake locations and magnitudes as well as focal mechanism solutions are determined as part of the routine SIL analysis. During the first months after the start of production in the Reykjanes field, the SIL network recorded three short-lived earthquake swarms SE of the tip of the peninsula, along the periphery of the subsidence bowl. The swarms occurred on 31 May–1 June, 9–10 July and 27–28 September 2006 (orange dots in Fig. 8), and each of them had 40–80 recorded events with a maximum local magnitude of ML3. A swarm occurred in the same area on 6–11 January 2008, and another swarm NW of the tip on 8 July 2008 (blue dots). The focal mechanisms, as determined by the SIL analysis, showed that the largest events were typically consistent with normal faulting on NEtrending planes (see Fig. 8), although there is some uncertainty on the 147 mechanisms as the earthquakes are located outside the seismic network. Very few earthquakes have previously been recorded in these areas by the SIL network, and the seismicity is clearly separated from the area of the intense 1972 swarm. We investigate if the unusual earthquake activity could be triggered by the crustal stresses caused by the contraction within the geothermal field. From our ellipsoid model for the subsidence around the Reykjanes geothermal field, we compute the change in Coulomb failure stress (ΔCFS) for normal slip on NE–SW trending fault planes dipping 60° toward NW, along a vertical profile AA′ shown in Fig. 8. The maximum stress change is close to 0.3 MPa, which may be enough to trigger earthquakes (e.g. King et al., 1994). The earthquakes recorded during the swarms in 2006 and 2008 were generally located a little deeper than the area of maximum ΔCFS. However, the uncertainty of the hypocentre depths in this area reported by the SIL catalogue are typically 3–6 km, and recent work on the velocity model on the Reykjanes Peninsula indicates that the earthquake depths here may be a little too deep (K. Vogfjörd, personal communication, 2008). Short-lived swarms also occurred SE of the Svartsengi field during 2–3 July 2007 (red dots) and 22–24 January 2008 (blue dots in Fig. 8). The swarms only had 40–60 recorded events, but they both included a number of ML N 3 events, and the January 2008 swarm included two ML4 events. Very little seismicity has been observed in the area since the start of production in the Svartsengi field in 1976, indicating that the pressure drawdown in the Svartsengi reservoir of around 3 MPa (Vatnaskil, 2009) have raised the fracture limit and thus reduced the micro-earthquake activity temporarily (Brandsdóttir et al., 2002; Keiding et al., 2009b). The small swarms along the periphery of the Svartsengi subsidence bowl may occur in response to the increased reinjection of waste fluids during recent years. 6.3. Aseismic faulting Interestingly, the radar data reveal subtle discontinuities that probably reflect fault movement, as can be seen in the near-vertical displacements along profile BB′ in Fig. 8. Near the NW end of the profile, we see a 3 km long NE-trending discontinuity, consistent with subsidence in a graben-like structure. This discontinuity becomes visible between two acquisitions from 14 September 2006 to 19 October 2006, a few months after the start of production in the Reykjanes field. No fractures have been mapped in this area, but the amplitude images from the radar data clearly show linear structures aligning with the edges of the discontinuity. The total offset in the Fig. 8. Close-up on the near-vertical radar displacement field during June 2005–May 2008 (same as in Fig. 5c). Earthquake locations and focal mechanisms from the SIL seismic catalogue are shown as background events (small black dots), and distinct swarm events in 2006 (orange), 2007 (red) and 2008 (blue). Also shown are focal mechanisms for some of the largest swarm events with local magnitudes ranging 2.9–4.1. The stippled outline shows the approximate location of the 1972 swarm activity (redrawn from Klein et al., 1977, their Fig. 5). Profile AA′ shows the predicted change in Coulomb failure stress, for normal slip on NE–SW trending fault planes, computed from the elastic halfspace ellipsoidal source model for the subsidence around the Reykjanes geothermal field. Profile BB′ shows the observed near-vertical radar displacement across the Reykjanes subsidence bowl. 148 M. Keiding et al. / Journal of Volcanology and Geothermal Research 194 (2010) 139–149 near-vertical radar image is approximately 1 cm. No earthquakes were associated with this fault motion. A similar, albeit less clear, discontinuity can be seen near the SE end of the profile. The discontinuities are located at some distance from the Reykjanes subsidence bowl, however they are within the areas of increased tensile stress, and the appearance of the largest discontinuity in fall 2006 suggests that these faulting events are also induced by the subsidence due to geothermal fluid extraction. 7. Conclusions We have examined the crustal deformation observed on the Reykjanes Peninsula during 1992–2009, using InSAR and GPS data. The plate boundary is transtensional with both left-lateral motion and extension, and a number of geothermal fields are located along the central plate boundary zone. The geodetic data show deformation due to plate motion, anthropogenic subsidence around the Reykjanes, Svartsengi and Hellisheidi geothermal power plants and, possibly, pressure increase in the Krísuvík geothermal system. We investigate the subsidence around the Reykjanes field in more detail, and estimate a maximum subsidence of around 10 cm during the first 2 years of production. The subsidence bowl around the Reykjanes field is clearly elongated in the NE–SW direction, thus aligning with the main trend of fractures in the area. The observed surface subsidence can be modelled to the first order using point or ellipsoidal pressure sources in an elastic halfspace. Following the start of production at the Reykjanes power plant, short-lived swarms of micro-earthquakes occurred along the SE and NW periphery of the subsidence bowl, in areas where few earthquakes have previously been recorded. The earthquake swarms, as well as aseismic fault movement revealed by the radar data, may be triggered by stress changes caused by geothermal fluid extraction at Reykjanes. Acknowledgements The ERS and Envisat data were provided by the European Space Agency. We thank Halldór Geirsson for providing the continuous GPS data, and Halldór Ólafsson for skilled and cheerful assistance during numerous GPS campaigns. The earthquake locations, magnitudes and focal mechanisms included in this study are from the SIL seismic catalogue courtesy of the Icelandic Meteorological Office. Páll Jónsson and Gudmundur Ómar Fridleifsson provided the pressure data from the Reykjanes geothermal field. We thank Maurizio Battaglia and Yuri Fialko for the codes for computing displacements and stresses due to an ellipsoidal source. Páll Einarsson, Grímur Björnsson, Ingvar Thór Magnússon and Ómar Sigurdsson are thanked for insightful comments. We are also grateful to Thomas R. Walter and an anonymous reviewer for constructive reviews that helped improving the paper. The figures were prepared using the GMT software (Wessel and Smith, 1998). This work is supported by a grant from the Eimskip Fund of the University of Iceland. Funding for GPS equipment used in this study came from the Icelandic Research Fund, the University of Arizona and NSF. References Allis, R., Bromley, C., Currie, S., 2009. Update on subsidence at the Wairakei–Tauhara geothermal system, New Zealand. Geothermics 38, 169–180. Altamimi, Z., Collilieux, X., Legrand, J., Garayt, B., Boucher, C., 2007. ITRF2005: a new release of the International Terrestrial Reference Frame based on time series of station positions and Earth Orientation Parameters. J. Geophys. Res. 112. Amelung, F., Galloway, D.L., Bell, J.W., Zebker, H.A., Laczniak, R.J., 1999. Sensing the ups and downs of Las Vegas: InSAR reveals structural control of land subsidence and aquifer-system deformation. Geology 27, 483–486. Amoruso, A., Crescentini, L., 2009. Shape and volume change of pressurized ellipsoidal cavities from deformation and seismic data. J. Geophys. Res. 114. Anderssohn, J., Wetzel, H.U., Walter, T.R., Motagh, M., Djamour, Y., Kaufmann, H., 2008. Land subsidence pattern controlled by old alpine basement faults in the Kashmar Valley, northeast Iran: results from InSAR, levelling and imaging spectrometry. Geophys. J. Int. 174, 287–294. Árnadóttir, T., Jiang, W., Feigl, K.L., Geirsson, H., Sturkell, E., 2006. Kinematic models of plate boundary deformation in southwest Iceland derived from GPS observations. J. Geophys. Res. 111. Árnadóttir, T., Lund, B., Jiang, W., Geirsson, H., Björnsson, H., Einarsson, P., Sigurdsson, T., 2009. Glacial rebound and plate spreading: Results from the first countrywide GPS observations in Iceland. Geophys. J. Int. 177, 691–716. Battaglia, M., Hill, D., 2009. Analytical modeling of gravity changes and crustal deformation at volcanoes: the Long Valley caldera, California, case study. Tectonophysics 471, 45–57. Björnsson, H., Thorgilsson, G., Halldórsdóttir, S., 2008. Numerical modelling of the geothermal system in Reykjanes, predictions for 50 MWe production increase from year 2011 (in Icelandic). Report ÍSOR-08053. Iceland Geosurvey, Reykjavík, Iceland. Bödvarsson, R., Rögnvaldsson, S.T., Slunga, R., Kjartansson, E., 1999. The SIL data acquisition system — at present and beyond year 2000. Phys. Earth Planet. Inter. 113, 89–101. Brandsdóttir, B., Franzson, H., Einarsson, P., Árnason, K., Kristmannsdóttir, H., 2002. Seismic monitoring during an injection experiment in the Svartsengi geothermal field, Iceland. Jökull 51, 43–52. Cervelli, P., Murray, M.H., Segall, P., Aoki, Y., Kato, T., 2001. Estimating source parameters from deformation data, with an application to the March 1997 earthquake swarm off the Izu Peninsula, Japan. J. Geophys. Res. 106 (11), 217–11238. Clifton, A.E., Kattenhorn, S.A., 2006. Structural architecture of a highly oblique divergent plate boundary segment. Tectonophysics 419, 27–40. Clifton, A.E., Sigmundsson, F., Feigl, K.L., Gudmundsson, G., Árnadóttir, T., 2002. Surface effects of faulting and deformation resulting from magma accumulation at the Hengill triple junction, SW Iceland, 1994–1998. J. Volcanol. Geotherm. Res. 115, 233–255. Dach, R., Hugentobler, U., Fridez, P., Meindl, M., 2007. Bernese GPS software version 5.0. Switzerland. Decriem, J., Árnadóttir, T., Hooper, A., Geirsson, H., Sigmundsson, F., Keiding, M., Ófeigsson, B.G., Hreinsdóttir, S., Einarsson, P., LaFemina, P., Bennett, R.A., 2010. The 29 May 2008 earthquake doublet in SW Iceland. Geophys. J. Int. DeMets, C., Gordon, R.G., Argus, D.F., Stein, S., 1994. Effect of recent revisions to the geomagnetic reversal time scale on estimates of current plate motions. Geophys. Res. Lett. 21, 2191–2194. Donnelly, L.J., 2009. A review of international cases of fault reactivation during mining subsidence and fluid abstraction. Quart. J. Eng. Geol. Hydrol. 42, 73–94. Dow, J.M., Neilan, R.E., Gendt, G., 2005. The International GPS Service (IGS): celebrating the 10th anniversary and looking to the next decade. Adv. Space Res. 36, 320–326. Eberhart-Phillips, D., Oppenheimer, D.H., 1984. Induced seismicity in the Geysers geothermal area, California. J. Geophys. Res. 89, 1191–1207. Efron, B., Tibshirani, R., 1986. Bootstrap methods for standard errors, confidence intervals, and other measures of statistical accuracy. Stat. Sci. 1, 54–75. Eysteinsson, H., 2000. Elevation and gravity changes at geothermal fields on the Reykjanes Peninsula, SW Iceland. Proceedings World Geothermal Congress 2000, pp. 559–564. Feigl, K.L., Gasperi, J., Sigmundsson, F., Rigo, A., 2000. Crustal deformation near Hengill volcano, Iceland 1993–1998: coupling between magmatic activity and faulting inferred from elastic modeling of satellite radar interferograms. J. Geophys. Res. 105, 25655–25670. Fialko, Y., Simons, M., 2000. Deformation and seismicity in the Coso geothermal area, Inyo County, California: observations and modeling using satellite radar interferometry. J. Geophys. Res. 105, 21781–21793. Glowacka, E., Gonzáles, J., Fabriol, H., 1999. Recent vertical deformation in Mexicali Valley and its relationship with tectonics, seismicity, and the exploitation of the Cerro Prieto geothermal field. Mexico. Pure Appl. Geophys 156, 591–614. Grasso, J.R., Wittlinger, G., 1990. Ten years of seismic monitoring over a gas field. Bull. Seis. Soc. Am. 80, 450–473. Herring, T.A., King, R.W., McClusky, S.C., 2006. GLOBK Reference Manual, Global Kalman filter VLBI and GPS analysis program, Release 10.3. Technical Report. Massachusetts Institute of Technology, USA. Hoffmann, J., Zebker, H.A., Galloway, D.L., Amelung, F., 2001. Seasonal subsidence and rebound in Las Vegas Valley, Nevada, observed by synthetic aperture radar interferometry. Water Resour. Res. 37, 1551–1566. Hooper, A., 2008. A multi-temporal InSAR method incorporating both persistent scatterer and small baseline approaches. Geophys. Res. Lett. 35. Hooper, A., Pedersen, R., Sigmundsson, F., 2009. Constraints on magma intrusion at Eyjafjallajökull and Katla volcanoes in Iceland, from time series SAR interferometry. In: B., C.J., et al. (Ed.), The VOLUME Project, VOLcanoes: Understanding subsurface mass moveMEnt, pp. 13–24. Hooper, A., Segall, P., Zebker, H., 2007. Persistent scatterer interferometric synthetic aperture radar for crustal deformation analysis, with application to Volcán Alcedo, Galápagos. J. Geophys. Res. 112. Hreinsdóttir, S., Árnadóttir, T., Decriem, J., Geirsson, H., Tryggvason, A., Bennett, R.A., LaFemina, P., 2009. A complex earthquake sequence captured by the continuous GPS network in SW Iceland. Geophys. Res. Lett. 36. Hreinsdóttir, S., Einarsson, P., Sigmundsson, F., 2001. Crustal deformation at the oblique spreading Reykjanes Peninsula, SW Iceland: GPS measurements from 1993 to 1998. J. Geophys. Res. 106, 13803–13816. Jakobsdóttir, S.S., 2008. Seismicity in Iceland: 1994–2007. Jökull 58, 75–100. Jónsson, P., Halldórsdóttir, S., Björnsson, H., 2009. Svartsengi–Reykjanes, temperature and pressure measurements 2008 (in Icelandic). Report ÍSOR-2009/036. Iceland Geosurvey, Reykjavík, Iceland. Jónsson, S., Segall, P., Pedersen, R., Björnsson, G., 2003. Post-earthquake ground movements correlated to pore-pressure transients. Nature 424, 179–183. M. Keiding et al. / Journal of Volcanology and Geothermal Research 194 (2010) 139–149 Jónsson, S., Zebker, H., Segall, P., Amelung, F., 2002. Fault slip distribution of the 1999 Mw7.1 Hector Mine earthquake, California, estimated from satellite radar and GPS measurements. Bull. Seis. Soc. Am. 92, 1377–1389. Kampes, B., Hanssen, R., Perski, Z., 2003. Radar Interferometry with Public Domain Tools, in: Proceedings of FRINGE 2003. ESA, Frascati, Italy. Keiding, M., Árnadóttir, T., Sturkell, E., Geirsson, H., Lund, B., 2008. Strain accumulation along an oblique plate boundary: the Reykjanes Peninsula, southwest Iceland. Geophys. J. Int. 172, 861–872. Keiding, M., Hooper, A., Árnadóttir, T., Jónsson, S., Decriem, J., 2009a. Natural and manmade deformation around geothermal fields on the Reykjanes Peninsula, SW Iceland, in: Proceedings of FRINGE 2009. ESA, Frascati, Italy. Keiding, M., Lund, B., Árnadóttir, T., 2009b. Earthquakes, stress and strain along an obliquely divergent plate boundary: the Reykjanes Peninsula, southwest Iceland. J. Geophys. Res. 114. King, G., Stein, R.S., Lin, J., 1994. Static stress changes and the triggering of earthquakes. Bull. Seis. Soc. Am. 84, 935–953. Klein, F.W., Einarsson, P., Wyss, M., 1977. The Reykjanes Peninsula, Iceland, earthquake swarm of September 1972 and its tectonic significance. J. Geophys. Res. 82, 865–888. Lu, Z., Masterlark, T., Power, J., Dzurisin, D., Wicks, C., 2002. Subsidence at Kiska Volcano, Western Aleutians, detected by satellite radar interferometry. Geophys. Res. Lett. 29. Magnússon, I.T., Thorbergsson, G., 2004. GPS measurements on the outer Reykjanes Peninsula 2004 (in Icelandic). Report IV Hluti. Iceland GeoSurvey, Reykjavík. McClusky, et al., 2000. Global Positioning System constraints on plate kinematics and dynamics on the eastern Mediterranean and Caucasus. J. Geophys. Res. 105, 5695–5719. Mogi, K., 1958. Relations between the eruptions of various volcanoes and the deformations of the ground surfaces around them. Bull. Earthquake Res. Inst. Univ. Tokyo 36, 99–134. Mossop, A., Segall, P., 1997. Subsidence at The Geysers geothermal field, N. California from a comparison of GPS and leveling surveys. Geophys. Res. Lett. 24, 1839–1842. Peltier, A., Hurst, T., Scott, B., Cayol, V., 2009. Structures involved in the vertical deformation at Lake Taupo (New Zealand) between 1979 and 2007: new insights from numerical modelling. J. Volc. Geotherm. Res. 181, 173–184. Sæmundsson, K., 1978. Fissure swarms and central volcanoes of the neovolcanic zones of Iceland. Geol. J. Special Issue 10, 415–432. 149 Segall, P., 1989. Earthquakes triggered by fluid extraction. Geology 17, 942–946. Segall, P., Fitzgerald, S.D., 1998. A note on induced stress changes in hydrocarbon and geothermal reservoirs. Tectonophysics 289, 117–128. Segall, P., Grasso, J.R., Mossop, A., 1994. Poroelastic stressing and induced seismicity near the Lacq gas field, southwestern France. J. Geophys. Res. 99, 15423–15438. Shirzaei, M., Walter, T.R., 2009. Randomly iterated search and statistical competency as powerful inversion tools for deformation source modeling: application to volcano interferometric synthetic aperture radar data. J. Geophys. Res. 114. Sigmundsson, F., Einarsson, P., Rognvaldsson, S., Foulger, G., Hodgkinson, K., Thorbergsson, G., 1997. The 1994–1995 seismicity and deformation at the Hengill triple junction, Iceland: triggering of earthquakes by minor magma injection in a zone of horizontal shear stress. J. Geophys. Res. 102, 15151–15161. Sudhaus, H., Jónsson, S., 2008. Improved source modelling through combined use of InSAR and GPS under consideration of correlated data errors: application to the June 2000 Kleifarvatn earthquake, Iceland. Geophys. J. Int. 176, 389–404. Tryggvason, E., 1973. Seismicity, earthquake swarms, and plate boundaries in the Iceland region. Bull. Seis. Soc. Am. 63, 1327–1348. Vadon, H., Sigmundsson, F., 1997. Crustal deformation from 1992 to 1995 at the midAtlantic ridge, southwest Iceland, mapped by satellite radar interferometry. Science 275, 194–197. Vasco, D.W., Wicks, C., Karasaki, K., Marques, O., 2002. Geodetic imaging: reservoir monitoring using satellite interferometry. Geophys. J. Int. 149, 555–571. Vatnaskil, 2009. Svartsengi–Reykjanes, production overview for 2008 (In Icelandic). Report 09.04. Hitaveita Sudurnesja. Reykjavík, Iceland. Werner, C., Wegmüller, U., Strozzi, T., Wiesmann, A., 2000. GAMMA SAR and interferometric processing software. Proc. ERS-ENVISAT Symp. European Space Agency, Gothenburg, Sweden. Wessel, P., Smith, W.H.F., 1998. New, improved version of the Generic Mapping Tools released. Eos Trans. AGU 79, 579. Wicks, C., Thatcher, W., Dzurisin, D., 1998. Migration of fluids beneath Yellowstone Caldera inferred from satellite radar interferometry. Science 282, 458–462. Yang, X.M., Davis, P.M., Dieterich, J.H., 1988. Deformation from inflation of a dipping finite prolate spheroid in an elastic half-space as a model for volcanic stressing. J. Geophys. Res. 93, 4249–4257. Zebker, H.A., Villasenor, J., 1992. Decorrelation in interferometric radar echoes. IEEE Transactions on Geoscience and Remote Sensing 30, 950–959.