3D seismic imaging of the Nesjavellir geothermal field, SW-Iceland

Abstract

We present a detailed seismic imaging of the harnessed Nesjavellir geothermal area, SW-Iceland, which is one of several geothermal fields on the flanks of the Hengill volcano. We map the , , and ratio using seismic data recorded in 2016–2020 and compare them with both a resistivity model of the same area and the rock temperature as measured in boreholes. The results show that the shallower crust (depth less than 1 km) is characterized by low and , and high ratio (around 1.9). Shallow low resistivity values at similar depths in the same area have been interpreted as the smectite clay cap of the system. At depths between 1 and 3 km the observed low ratio of 1.64–1.70 correlates with high resistivity values. In this area, characterized by temperatures larger than 240°C, the seismicity appears to be sparse and located close to the production wells. This seismicity has been interpreted as induced by the production in combination with naturally occurring earthquakes. At depths greater than 4 km, high ratio of 1.9 correlates well with low resistivity values. In the valley of Nesjavellir, a deep-seated conductive body, domes up at about 4.500 m b.sl. and coincides spatially with a significant high ratio anomaly (>1.9). Above these anomalies an elevated temperature is registered according to borehole temperature data. This is proposed here to be caused by hot 600°C–900°C cooling intrusives, close to the brittle ductile transition—probably the heat source(s) of the geothermal field above. These anomalies are at the same location as the last fissure eruption in Hengill almost 2,000 years ago. The NNE-SSW trending, deeper seismic cluster at 3–6 km depth is located at the edge of this high anomaly. The heat source of the Nesjavellir geothermal field is most likely connected to this most recent volcanism as reflected by the deep-seated low resistivity body and high ratio, located beneath the deep fault that connects the flow path of the high temperature geothermal fluid, resulting in an actively producing reservoir.

1 Introduction

The Hengill volcano, with its surrounding geothermal fields, is the most productive harnessed high temperature geothermal region in Iceland. The Nesjavellir power plant (Figure 1), to the north of the Hengill volcano (120 MWe and 350 MWt), was commissioned in 1990 and the Hellisheiði power plant, to the southwest of the volcano (303 MWe and 210 MWt), was commissioned in 2006. The area is located in the southwest region of Iceland, about 30 km east of the capital Reykjavík at a triple junction formed by the Reykjanes Peninsula (RP), the Western Volcanic Zone (WVZ) and the South Iceland Seismic Zone (SISZ) (Stefánsson et al., 1993; Einarsson, 2008).

FIGURE 1

The geology of the Hengill area is dominated by a succession of alternating layers made of hyaloclastite and basalt deposits formed during sub-glacial and subaerial volcanic eruptions, respectively (Sæmundsson, 1967; Árnason et al., 1986; Sæmundsson, 1992). Below 1.5 km depth, the subsurface is dominated by basaltic intrusions, dykes, and sills (Gunnarsdóttir et al., 2020 and references therein). Hyaloclastites have formed ridge-shaped mountains and tuyas, including the Hengill volcano itself (Sæmundsson et al., 2016). The focus of this study is the Nesjavellir geothermal sub-field (Figure 1). It lies within the N30°E trending Hengill fissure swarm, which is in total about 40 km long, parallel to the hyaloclastite ridges, faults and fissures of the Hengill fissure swarm (Sæmundsson, 1992; Sæmundsson et al., 2016).

The geothermal activity of the Hengill area is closely connected to the three volcanic systems in the region: Grensdalur (extinct), Hrómundartindur (last eruption 11,000 years ago) and Hengill (last eruption 2,000 years ago) (see, Árnason et al., 1987; Sæmundsson, 1992 and reference therein). The conceptual model of the Hengill geothermal system was assumed to consist of a single central up-flow underneath the Hengill volcano (e.g., Björnsson et al., 2003). However, more recent studies indicate a more complex picture with multiple local heat sources with cooler regions in between (Gunnarsson et al., 2011; Gunnarsson & Aradóttir, 2015). The geothermal fluids migrate through the porous basaltic crust, fractures and fissures in the area where the two most productive power plants in Iceland are located: the southern Hellisheiði and the northern Nesjavellir power plants, both operated by ON Power (https://www.on.is/), a subsidiary of OR-Reykjavík Energy. The reservoir temperature for the Hengill area ranges from 200°C to 340°C (Gunnarsson & Aradóttir, 2015). The energy is used to provide space heating and for other direct uses in Reykjavík and neighboring communities, as well as for electricity generation. In the Nesjavellir sub-field a total of 30 production wells have been drilled. Part of the surplus geothermal water from the plant goes into injection wells.

The seismic activity in the Hengill area has been monitored since the middle of 1970s (e.g., Foulger, 1988a, b). Earthquake detection was much improved in 1990 through the installation of the denser South Iceland Lowland seismic network (SIL) of the Icelandic Meteorological Office (IMO) (Stefánsson et al., 1993; Jakobsdóttir, 2008). The permanent seismic network in the Hengill area now consists of 10 stations operated by the Iceland GeoSurvey (ÍSOR) for ON Power, and 5 stations operated by the IMO. The seismicity is characterized by persistent tectonic background activity (Bodvarsson et al., 1996; Rögnvaldsson et al., 1996; Brandsdóttir et al., 2010; Li et al., 2019), and earthquake swarm activity in connection with magma intrusion and geothermal fluid migration in hot rocks (Sigmundsson et al., 1997; Jakobsdóttir, 2008).

Drilling and fluid injection operations can induce seismicity and thus increase the already high rate of seismicity, particularly in the production and injection areas of the Hellisheiði and Nesjavellir power plants (e.g., Flóvenz et al., 2015; Ágústsson et al., 2015; Kristjánsdóttir et al., 2019; Hjörleifsdóttir et al., 2021a, b). Close to the Hellisheiði geothermal plant (64°02′14″N 21°24′03″W), seismicity occurred already during drilling of the injection-wells in the Húsmúli formation, when the wells intersected permeable formations or the bottom of the feed zones with maximum magnitude of M_L 2.5 (Ágústsson et al., 2015). When injection began in Húsmúli in September 2011 the rate of seismicity greatly increased, including swarms and two large events, M_L 4.0 and 3.8, in the middle of October 2011 (Bessason et al., 2012). The reinjection also caused a local uplift of 2 cm during the first year (Juncu et al., 2017). In the southernmost Hengill sub-field, Hverahlíð, seismicity increased significantly in 2016 when drilling and production started, and has persisted to date (Kristjánsdóttir et al., 2019; Hjörleifsdóttir et al., 2021b).

Detailed analysis of the elastic properties of the crustal media in relation to the local seismicity are important to better understand the crustal structure as well as the current status and permeability of geothermal systems and volcanoes. Variations of P-wave velocity (), S-wave velocity () and of the ratio are mainly thought to depend on: lithology, rock porosity, hydrothermal alteration, fluid composition, pore fluid pressure and saturation, and temperature (e.g., Tryggvason et al., 2002; Mavko et al., 2009). Monitoring of harnessed geothermal areas with a dense seismic network provides the opportunity to record changes in the seismicity rate and the application of local earthquake seismic tomography. This imaging method of the subsurface can complement and better constrain other independent geological and geophysical analysis, such as gravity and resistivity modeling widely applied in geothermal areas (e.g., Ussher et al., 2000; Árnason et al., 2010; Hersir et al., 2022).

Several 1D and 3D velocity models have been made on a regional scale for the greater Hengill area, starting with Foulger (1984) and Foulger and Toomey (1989). Tryggvason et al. (2002) used about 7,500 local earthquakes recorded at 20 seismic stations of the SIL network in a 224 km² × 112 km² wide area, including the South Iceland Seismic Zone, the Hengill volcanic system and the Reykjanes Peninsula, to model the 3D velocity structure down to 10–15 km depth. More recently, Jousset et al. (2011) performed a joint tomographic inversion for the estimation of 3D P- and S-wave velocity structure and local earthquakes in the whole Hengill geothermal area including its subareas. They used 250 local earthquakes recorded from June to October 2006 by 8 seismic stations. Wagner (2019) did local-earthquake tomography of Hengill. His , , and models show structures that can be associated with the local geothermal activity. Finally, Sánchez-Pastor et al. (2021) performed ambient seismic noise tomography of the entire Hengill area, with a good resolution in the uppermost 4 km of the crust.

The results presented here aim to provide a more accurate and detailed 3D velocity structure of the Nesjavellir geothermal sub-field using high-resolution seismic tomography and relocation of the local earthquakes recorded from 18 October 2016 to 10 June 2020. The results are then compared with a resistivity model and rock temperature data from boreholes, taking the local geology into account. This study is an important contribution to a better understanding of the crustal structure in and around the Nesjavellir geothermal area, as well as addressing the possible mechanisms generating the observed seismicity.

2 Data

Seismicity in the Hengill area has been monitored by several independent permanent and temporary seismic networks. For almost 3 decades, the activity was monitored by the regional seismic network (the SIL network) operated by the Icelandic Meteorological Office (IMO) (Figure 1, red diamonds). In 2016, Iceland GeoSurvey (ÍSOR) installed a seismic network consisting of 10 stations (Figure 1, blue diamonds), for ON Power (ON) in the Hengill area, which drastically improved the earthquake detection threshold of the pre-existing SIL network. In 2018–2021, the Geothermica EU funded COSEISMIQ project (COntrol SEISmicity and Manage Induced earthQuakes), operated a dense network of 23 stations in the area (Figure 1, green diamonds). In 2019–2020, two additional stations were set up in Nesjavellir by the European Union’s Horizon 2020 funded Science4CleanEnergy (S4CE) project (Figure 1, magenta diamonds). This adds up to a total of 38, three component seismic stations monitoring the Hengill area in 2019–2020. For the purpose of this study, we used 14 permanent and temporary stations in the vicinity of Nesjavellir, comprising 5 ON/ÍSOR stations (Lennartz LE-3DliteMkII 1 s), 2 S4CE stations (Lennartz LE-3DliteMkII 1 s), 1 IMO station (Lennartz 5 s) and 6 COSEISMIQ stations (3 STS2s 120 s and 3 Guralp 6D 30 s) (Figure 1).

The available dataset consists of 6,906 earthquakes extracted from ÍSOR’s catalog, of local magnitude −0.8 ≤ ML ≤ 3.8, recorded between 18 October 2016 and 10 June 2020, in a 120 km² area around the Nesjavellir geothermal field (Figure 1). ÍSOR’s manually refined catalog consists of 31,320 P- and 25,624 S-waves travel times recorded at the 14 permanent and temporary seismic stations used in our study. Travel time errors associated with the P- and S-wave manual picks are 0.01 and 0.02 s, respectively. In order to evaluate the P- and S-wave travel time consistency, we applied the Wadati diagram analysis (Wadati and Oki, 1933). Specifically, we plotted the difference T_S-T_P between S-waves (T_S) and P-waves (T_P) travel times, as a function of T_P -T₀ (T₀ is the origin time) for each event recorded at the available stations (Figure 2). The resulting Wadati diagram shows a linear trend, revealing the good consistency of the data. From the best fitting line, we estimated a local ratio equal to (for the extent of Figure 1B). This value is fully consistent with the ratio computed on a larger scale in the wider Hengill region: Foulger (1995) and Miller et al. (1998) found , Tryggvason et al. (2002) obtained , Jousset et al. (2011) found .

FIGURE 2

3 Method and inversion strategy

In order to seismically image the study area we performed a 3D local earthquake tomography. We first selected the best constrained catalog events (Section 3.1), based on the quality of the preliminary location, to calculate a minimum 1D velocity model. Using the new velocity model and the estimated value of 1.77, we performed a first relocation of all the events contained in the catalog. Based on the quality of the relocations, we selected the highest quality events and used them to perform the 3D tomographic inversion. In the following subsections we describe the details of the methodologies and the related inversion strategies

3.1 1D velocity model and earthquake relocation

Reliable initial hypocenter locations and a reference seismic velocity model are crucial to successfully perform local earthquake tomography (Kissling, 1988; Kissling et al., 1994). For this reason, a reliable minimum 1D P-wave velocity model for the study area has been obtained solving the coupled hypocenter-velocity problem by implementing the VELEST software (Kissling et al., 1995). For this purpose, we used only the P-wave arrival times of the best constrained 473 earthquakes, selected from the starting dataset by requiring a minimum number of 5 P-phases (and at least 2 S-phases), and an azimuthal gap smaller than 180°. Several different velocity models have previously been proposed for Southwest Iceland using older datasets (Figure 3A): 1) the SIL model (Bjarnason et al., 1993; Stefánsson et al., 1993) is routinely used for earthquake locations in Iceland; 2) the velocity model by Tryggvason et al. (2002), obtained by using data recorded in Southwest Iceland, comprising the Reykjanes Peninsula, the Hengill area and the SISZ; 3) the Jousset et al. (2011) model, which includes seismicity of the Hengill area (including the Nesjavellir and Hellisheiði subareas). These models, together with our gradient velocity model were all tested in the present study as starting velocity models in VELEST (Figure 3A) and used to produce four final models (solid coloured lines in Figure 3B). Additionally, the average of these four models (dashed black line in Figure 3B) is used as a starting model for an additional inversion step whose solution finally represents the best minimum 1D P-wave velocity model (indicates as “Final Model” in Figure 3C).

FIGURE 3

We used a non-linear global search location method, implemented in the NonLinLoc software (Lomax et al., 2000), to perform the relocation of the collected 6,906 events using the best minimum 1D P-wave velocity model (Figure 3C) and the equal to 1.77 as inferred from the Wadati diagram (Figure 2). We retained all those events whose relocation was obtained by using at least 6 P-wave and 2 S-wave arrival time readings, maximum gap of 180°, root mean square residual lower than 0.3 s, and horizontal and vertical location error lower than 0.6 km. The result of this selection is a dataset of 2,358 seismic events shown in Figure 1, which is then used as input for the 3D tomographic inversion.

3.2 3D P and S tomography

We used an iterative linearized delay-time inversion to estimate both the 3D P- and S-wave velocity models and earthquake locations (Latorre et al., 2004). The tomographic code has been used for several applications in tectonic (Amoroso et al., 2014, 2017; Napolitano et al., 2021), volcanic (Serlenga et al., 2016; De Landro et al., 2017), and geothermal areas (Amoroso et al., 2018), as well as in feasibility study in off-shore exploration areas (De Landro et al., 2020). First arrival times are computed through a finite-difference solution of the eikonal equation (Podvin and Lecomte, 1991) on a fine grid of 100 m³ × 100 m³ × 100 m³. Travel times are recalculated by numerical integration of the slowness along the rays. The parameters, that is, , , and earthquake location, are inverted by using the least square root (LSQR) method of Paige and Saunders (1982). The number of inversion steps was set to a maximum of 10 and the regularization of the inversion is achieved through a damping factor, which was calibrated using the L-curve empirical approach (Supplementary Figure S1). The spatial distribution of stations and events allowed us to investigate a volume of 12 km³ × 10 km³ × 13 km³ (Figure 1B). The velocity distribution in a continuous medium is described by a trilinear interpolating function based on a grid of regularly spaced nodes. The optimal grid mesh has been chosen taking into account the dimension of the investigated volume, the source/station geometry and corresponds to 0.5 km³ × 0.5 km³ × 0.25 km³. The root-mean-square function (RMS), defined as the weighted sum of the squared time delay, is analyzed a-posteriori. It has been verified that the convergence was reached after 5 iterations with RMS reduction of 80% (Supplementary Figure S2). To determine the well-resolved area, we computed the Derivative Weight Sum (DWS) and performed checkerboard tests (Supplementary Figure S3). The DWS threshold was conservatively chosen according to the checkerboard result images.

We have estimated the P- and the S-phase pick errors to be equal to 0.01 and 0.02 s, respectively. The reliability of the results has also been tested with respect to the quality of the initial dataset. The test consisted of performing a set of tomographic inversions considering different perturbed datasets. The results indicate that the main features identified in the tomographic images do not change when considering perturbed dataset using standard deviation less or equal to 0.05 s for P travel-time and 0.1 s for S, respectively (Supplementary Figure S4).

4 Results

The minimum 1D P-wave velocity model for the Nesjavellir geothermal area and surroundings (Figure 3C) shows two large velocity steps at 1 and 3 km depth of more than 1 km/s each, but an overall almost linear increase from 3.32 km/s at surface to 7.80 km/s at 10 km depth (Table 1).

The final 3D velocity model is represented through horizontal slices at several depths (b.s.l.) in Figures 4–6 and in cross-sections in Figure 7. Note that, while the and the models are the direct output of the tomographic code, the model is derived by dividing the values of the and models at each node of the inversion grid. As indicated by the DWS analysis (outlined with a gray line), the best resolution is attained at depths ranging between 1 and 5 km, in agreement with the earthquake locations. The and the values in the 3D model, considering only the well-resolved area inferred by the DWS, vary from 3.2 to 6.7 km/s and from 2.3 to 3.6 km/s, respectively. From surface to 3 km depth, a low anomaly (labeled as PL in Figure 4), initially emerges in the NW part of the map and extends down to 4 km depth. The model shows a low-velocity anomaly from 1 km down to 2.5 km depth (labeled as SL in Figure 5) in the northern part of the model. As for the model, the well-resolved area features values less than 1.77 down to 2.5 km depth (Figure 6).

FIGURE 4

FIGURE 5

FIGURE 6

For depths greater than 3 km, the presence of a high P-wave velocity anomaly (PH) is detected in the northwest sector in addition to PL (Figure 4). In the same depth range does not change significantly (Figure 5). The anomalies in the model produce a notable complex pattern in the model down to 4 km. In particular, at 3 km depth, we note two localized areas of high ratio (PSH) and a low ratio (PSL) (Figure 6). At depth greater than 4 km, the well-resolved area has a limited extent in accordance with the earthquake distribution, where the ratio reaches the highest value of 2.1 in the Nesjavellir valley, slightly southwest of the power plant (labeled Ne). Cross-sections AA′ and BB’ on Figure 1 for the , and ratio in depth view are displayed on Figure 7.

FIGURE 7

The distribution of the relocated seismicity using the new detailed 3D and models, obtained applying the NonLinLoc software (Lomax et al., 2000), is shown on Figures 1, 8. The seismicity is clustered in two different groups (Figure 1): one to the south, between seismic stations SKA10, OLK26, KAT03, and BIT06, located between the Hengill and Hrómundartindur volcanic centers; and the second one to the north in the Nesjavellir valley, between stations SKEGG, GRAFN, NESJV, and KOLDU. The events belonging to the first group are part of the seismicity affecting a larger region outside the study area and not well covered by the seismic stations used in this study. Earthquakes belonging to the second group (Figure 8), are located in the Nesjavellir production area and are characterized by a well-defined NNE-SSW direction in agreement with the NNE-trending system of normal faults (Helgadóttir et al., 2010; Sæmundsson et al., 2016). The seismicity in Nesjavellir (second group) has a clear dual depth distribution (Figure 8). The shallow seismicity, at 1.0–3.0 km depth, is located to the northwest of the Nesjavellir valley. The deeper seismicity, at 3–6 km depth, is located to the southeast of the valley and aligned with the NNE-SSW trending surface faults. The deeper events seem to be located on at least one west-dipping structure. We note that the shallower seismicity occurs in an area characterized by low value (PL), low value (SL), and low ratio (PSL) (Figures 4–7). The seismicity at depths greater than about 3 km is located in a high (PH) and high ratios (PSH) area.

FIGURE 8

5 Discussion

In the following section we will discuss and interpret the results of the high resolution 3D tomographic images of the Nesjavellir geothermal area we provided, in terms of , , and and of the detailed earthquake relocations obtained by using these new 3D models. Then, we will compare velocity, resistivity and well temperature data in order to achieve the most reliable interpretation of the crustal structure surrounding the Nesjavellir geothermal field.

5.1 Seismic tomography

At shallow depths (<1 km), the observed low (4.0 km/s, Figure 4) and the still lower (2.3 km/s, Figure 5), which results in high ratio values (1.9, Figure 6), suggest the presence of porous material characterized by the presence of fluid. A similar value has been found also by Duran (2021). Although this depth range is at the edge of the well-resolved area, the presence of porous hyaloclastites that have been confirmed with drilling, is reflected in low velocities (Sæmundsson, 1967; Árnason et al., 1986; Sæmundsson, 1992).

Between 1.0 and 3.0 km depth, the production area (Figures 1C, 7) is marked by well constrained low ratios (1.64–1.70), due to low (5.2 km/s) and low-to-average (2.9 km/s) values compared to the 1D starting velocity model. This result has been also found by Wagner (2019). We assume that these velocity anomalies are related to effects of the geothermal system on the shallow crustal structure, possibly in combination with field operations. The geothermal system is characterized by higher than average temperature, clay minerals and over pressurized fluids. From lab experiments on Nesjavellir chlorite altered hyaloclastite, Jaya et al., 2009 found out that Vp decreases with increasing temperature.

The extensive fluid extraction at Nesjavellir lowers the pore-pressure in the field and consequently increases the steam dominated zone, leading to lower ratios. The P-wave velocity decreases as the fluid changes from water to steam, while the S-velocities are less affected. This is because changes in the compressibility of pore contents (fluid/gas) alter the rock’s compressibility while having no effect on its shear resistance (Tryggvason et al., 2002). This is most evident at 1.5–3.0 km depth in the cross-sections. The presence of the vapor zone is also corroborated by surface manifestations, such as fumaroles and steaming ground. A vapor steam cap is likely to be located in the production area above 1,500 ml. Additionally, vapor may be present in extremely small volumes, very locally around each well in their feed zones.

At depth greater than 3 km, down to the bottom of the well resolved area (5 km depth), the tomographic images show high (PH) and low-to-average values, resulting in a high ratio anomaly (labeled PSH). The highest value of is larger than 1.9 and is reached at 4.5 km (the retrieved value is also in accordance with the results obtained by Duran, 2021), at the tip of the deeper cluster of seismicity (Figure 7). These values indicate high temperatures or melt fractions (e.g., Tryggvason et al., 2002). A more thorough discussion of these values will be taken up in Section 5.3.

5.2 Relocated seismicity

Figure 8 displays the 2,358 relocated seismic events in the 3D final P- and S-wave models. The vertical and horizontal error modal value is 200 m. The seismicity is confined to the Nesjavellir valley and surrounding ridges, which contains the Nesjavellir power plant, its production and injection wells. The bimodal depth distribution of the seismicity is striking, with the shallower cluster located at 0.5–3.0 km depth and the deeper cluster at 3–6 km depth (Figure 8). The shallowest part of the deeper cluster is located slightly to the NE of the power plant (labeled Ne, on Figure 8) reaching under the power plant and production area at greater depth. We hypothesize that this NNE-SSW trending, west (∼40°) dipping, deeper cluster (3–6 km) is of tectonic nature, as it is highly unlikely that the production would influence stresses down to 6 km depth. It is more likely the seismicity is taking place on the NNE trending normal faults mapped in the area, which belong to the fissure swarm (on-land rift segment) of the Hengill volcano.

It is well established that the reinjection of cold wastewater in hot rocks can yield in vapor creation, and giving rise to cracks opening due to thermal contraction, which together with a fluid pore pressure increase may be the cause of the shallower seismicity. This mechanism has been invoked for other geothermal areas such as The Geysers in southern California (e.g., Majer & Peterson, 2007; Lin & Wu, 2018), Krafla NE-Iceland (Ágústsson and Blanck, 2019) and on Hellisheiði (Ágústsson et al., 2015). However, prior to 2020, all of the wastewater in Nesjavellir was injected into wells shallower than 500 m, and the injection fluid is, therefore, not thought to be in direct contact with the main underlying geothermal system. The re-injection of wastewater into well NJ-18 (the deepest well; for location, see Figure 8) only started in 2020 and no change in seismicity has been observed in association with that injection. It is, therefore, highly unlikely that the shallow seismicity observed during our study period in the Nesjavellir valley is induced by the reinjection. On the other hand, seismicity can also be induced by production, as demonstrated by Kristjánsdóttir et al. (2019) in Hverahlíð, another sub-field of the Hengill geothermal area. We suggest that the rather sparse shallower seismicity is more likely related to the production, in combination with naturally occurring earthquakes, and taking place within the geothermal reservoir. Well temperature data (Figure 9) show that the shallow seismicity is observed in areas of temperature >240°C. It is an interesting observation how little shallow seismicity is observed with a dense seismic network despite the intense production of 120 MWe and 350 MWt.

FIGURE 9

5.3 Joint interpretation of the 3D seismic tomography, resistivity and well temperature data

The combination of local high resolution seismic and resistivity tomographic images provides a powerful tool for mapping the structure of the subsurface and further understanding of the role of the fluids in the production area of the Nesjavellir geothermal field. Notably, for the first time a joint interpretation or comparison of resistivity data, ratio and well temperature data is applied (Figure 9). In this respect, we compare the results obtained in this study to a resistivity model, based on 1D joint inversion of Transient Electromagnetic (TEM) and Magnetotelluric (MT) data (Vilhjálmsson et al., 2019) and formation temperature estimated from down-hole measurements (Gunnarsson et al., 2011) (Figure 9). The formation temperature data are available only in the depth range of the wells in the area, or down to 1.500–2.000 m, with measured temperatures of up to >380°C.

Resistivity methods are commonly used in geothermal studies as a temperature proxy in the exploration of high-temperature geothermal fields, keeping in mind that the system might well be fossil and because of that or other reasons, the temperature may not be in equilibrium with the formation temperature (see e.g., Hersir et al., 2022 and references therein). Low resistivity (<10 Ωm), in the uppermost 1 km of the crust, is assumed to be related to smectite, which is a low temperature hydrothermal alteration mineral forming at 100°C–230°C (Árnason et al., 1986, 1987; Hersir et al., 1990; Árnason et al., 2010; Benediktsdóttir et al., 2021). Beneath this conductive clay cap, a high-resistivity core is associated with high-temperature alteration minerals (chlorite, epidote; formed at temperature >230°C). In the Hengill high-temperature geothermal field the up-doming low resistivity at depth (<10 Ωm) imaged by MT data (Árnason et al., 2010; Benediktsdóttir et al., 2021) is thought to represent the core of the geothermal field and presumably reflecting hot, solidified intrusions that are heat sources for the geothermal system.

Cross-section A-A′ on Figure 9, which lies across the Nesjavellir valley (for location, see Figure 1), shows both up-doming of the shallow low resistivity layer (the smectite clay cap) and the up-doming of the deep-seated low resistivity layer (between well NJ-12 and NJ-24), reaching up to about 4.500 m b.s.l. (upper left panel on Figure 9). The deep-seated low resistivity anomaly coincides spatially quite well with the high ratio anomaly (>1.9). This is also the location of the NNE-SSW trending, west (∼40°) dipping seismic cluster, discussed above. The formation temperature on cross-section A-A’ domes up right above the deep-seated resistivity low and high anomalies (Figure 9, upper right panel). These anomalies are at the same location as the last fissure eruption in Hengill almost 2,000 years ago (Sæmundsson, 1962; Sinton et al., 2005).

Similar correlations are seen on cross-section B-B’ (Figure 9), which lies along the Nesjavellir valley (for location, see Figure 1). The deep-seated low resistivity which extends from the Hengill volcano to the north domes up locally (lower left panel on Figure 9) coinciding with the high ratio anomaly (>1.9). The formation temperature domes up to shallower depths above these two anomalies (lower right panel on Figure 9). At 1.000–3.000 m b.s.l., a relatively low ratio anomaly (<1.72) is seen on both cross-sections. This could partly be due to the extensive fluid extraction at Nesjavellir which lowers the pore-pressure in the field and consequently increases the steam dominated zone, leading to lower ratios. Or alternatively, by the high temperatures and alteration.

At depths greater than 3 km, down to 5 km depth, our results indicate that the high ratio is due to an increase in that dominates with respect to the decrease in (Figure 7). This anomalously high ratio (>1.9) correlates well with the deep-seated low resistivity on both profiles as shown on Figure 9 and discussed above. A high ratio is in general to be expected when the rock temperature approaches its melting point, the shear strength reduces, and S-waves are attenuated far more than the P-waves. High ratio indicates high temperature, but very little, if any, partial melt as we still record S-waves. When temperature reaches more than 600°C, resistivity decreases significantly as the mineral conduction starts to play a big role and pore fluid and surface condition become of less importance (see e.g., Hersir et al., 2022 and references therein). The brittle/ductile transition in the northern Hengill area and Nesjavellir has previously been established at around 6–7 km (Tryggvason et al., 2002; Li et al., 2019). Our observations (absence of seismicity) are in agreement with these previous studies, marking the brittle ductile transition at 6 km depth (Figure 8). The brittle ductile transition for basaltic rocks is thought to lie at temperatures of 600°C–700°C (Ágústsson and Flóvenz, 2005; Violay et al., 2012). Seismic velocities will probably not change much until close to 900°C. In Krafla high-temperature area, NE-Iceland, magma temperature was estimated to be 950°C (Zierenberg et al., 2013). Laboratory experiments indicate that the temperature there seems to increase very fast from 600°C at the brittle ductile transition to the magma temperature of 950°C. Therefore, it is proposed here that the high ratio and deep-seated conductivity anomalies most likely reflect hot 600°C–900°C cooling intrusions, close to the brittle ductile transition—probably the heat source(s) of the geothermal field above. This is also in accordance with the elevated temperature right above these two anomalies—an upflow of geothermal heat.

Figure 10 is a horizontal slice at 4 km depth through the resistivity model showing earthquake locations and the ratio. It is most interesting, that the NNE-SSW trending seismicity is located outside, on the border, of the high ratio >1.85 anomaly, where the low resistivity (10 Ωm) is doming up. This is due to the proximity to the ductile part of the crust.

FIGURE 10

6 Conclusion

We present the first 3D seismic tomographic study specifically focused on the Nesjavellir geothermal field, and its associated natural and induced seismicity.

The relocated seismicity using the new 3D velocity model highlights the dual earthquakes depth distribution that was not discussed by previous large scale tomographic studies: 1) rather sparse seismicity within the production area of Nesjavellir at 1–3 km depth in relation to the intense production of 120 MWe and 350 MWt. This seismicity is located around the bottom of the wells and is likely induced by the production in combination with natural events taking place in the shallow part of the geothermal system. The shallow seismicity is observed in areas of temperature >240°C. We speculate that the seismicity is triggered within the part of the reservoir which is being extracted from. If that is the case then the seismicity outlines the harnessed reservoir. The joint interpretation for shallower depths (1–3 km) shows that the low of <1.72 correlates with the high resistivity values around the bottom of the production wells at 2 km depth, perhaps suggesting the presence of steam, probably partly because of the production due to lowered pore-pressure, or alternatively high temperatures and alteration. This can be supported by the observation that at this depth the formation temperature is measured 280°–380°. The deeper earthquake cluster at 3–6 km depth, located in the eastern part of the Nesjavellir valley is most likely of tectonic origin on one or more NNE-SSW trending faults, dipping W, within the Hengill fissure swarm.

The comparison of variations in crustal velocities in relation to resistivity and well temperature data in the valley of Nesjavellir, shows a deep-seated conductive body, which domes up at about 4.500 m b.sl. and spatially coincides with a significant high ratio anomaly (>1.9). Right above these two anomalies an elevated temperature is seen in borehole temperature data. It is proposed here that this is caused by hot 600°C–900°C cooling intrusions, close to the brittle ductile transition—probably the heat source(s) of the geothermal field above. These anomalies are at the same location as the last fissure eruption in Hengill almost 2,000 years ago. The NNE-SW trending, deeper seismic cluster at 3–6 km depth is located at the edge of this high anomaly. We propose a simple conceptual model that the heat source of the Nesjavellir geothermal field is connected to this most recent volcanism as reflected by the deep-seated low resistivity body and high ratio, located beneath the deep fault that connects the flow path of the high temperature geothermal fluid, resulting in an actively producing reservoir.

The 3D Nesjavellir optimized velocity model presented in this study can be of relevant importance for future applications. This is the case, for example, of the precise earthquake location aimed at monitoring the space and time evolution of the seismicity that can involve active faults able to generate large magnitude earthquakes. Moreover, the availability of a 3D model represents a starting point for a future 4D tomography study which will allow us to track changes in crustal properties over time. Eventually, the estimation of fault mechanisms and kinematic source parameters can of course benefit from the proposed 3D model.

References

• 1

  ÁgústssonK.BlanckH. (2019). Krafla, jarðskjálftar og niðurdæling (Krafla, earthquakes and reinjection). In Icelandic. ÍSOR-2019/022, LV-2019-032.

  • Google Scholar

• 2

  ÁgústssonK.FlóvenzÓ. G. (2005). “The thickness of the seismogenic crust in Iceland and its implications for geothermal systems in Iceland,” in Proceedings of the 2005 world geothermal congress (Turkey: Antalya).

  • Google Scholar

• 3

  ÁgústssonK.KristjánsdóttirS.FlóvenzÓ. G.GuðmundssonÓ. (2015). “Induced seismic activity during drilling of injection wells at the Hellisheiði power plant, SW Iceland,” in Proceedings world geothermal congress 2015 melbourne (Australia, 19–25.

  • Google Scholar

• 4

  AmorosoO.AscioneA.MazzoliS.VirieuxJ.ZolloA. (2014). Seismic imaging of a fluid storage in the actively extending Apennine mountain belt, southern Italy. Geophys. Res. Lett.41, 3802–3809. 10.1002/2014GL060070

  • CrossRef
  • Google Scholar

• 5

  AmorosoO.FestaG.BrunoP. P.D’AuriaL.De LandroG.Di FioreV.et al (2018). Integrated tomographic methods for seismic imaging and monitoring of volcanic caldera structures and geothermal areas. J. Appl. Geophy.156, 16–30. 10.1016/j.jappgeo.2017.11.012

  • CrossRef
  • Google Scholar

• 6

  AmorosoO.RussoG.De LandroG.GaramboisS.MazzoliS.ParenteM.et al (2017). From velocity and attenuation tomography to rock physical modeling: Inferences on fluid driven earthquake processes at the irpinia fault system in southern Italy. Geophys. Res. Lett.44, 6752–6760. 10.1002/2016gl072346

  • CrossRef
  • Google Scholar

• 7

  ÁrnasonK.EysteinssonH.HersirG. P. (2010). Joint 1D inversion of TEM and MT data and 3D inversion of MT data in the Hengill area, SW Iceland. Geothermics39 (1), 13–34. 10.1016/j.geothermics.2010.01.002

  • CrossRef
  • Google Scholar

• 8

  ÁrnasonK.HaraldssonG. I.JohnsenG. V.ThorbergssonG.HersirG. P.SaemundssonK.et al (1987). Nesjavellir - Ölkelduháls. Surface investigation in 1986. Orkustofnun report, in Icelandic, 112. OS-87018/JHD-02.

  • Google Scholar

• 9

  ÁrnasonK.HaraldssonG. I.JohnsenG. V.ThorbergssonG.HersirG. P.SaemundssonK.et al (1986). Nesjavellir. Geological and geophysical surveys in 1985 Orkustofnun report, in Icelandic, 125. OS-86014/JHD-02.

  • Google Scholar

• 10

  BenediktsdóttirÁ.ÁrnasonK.KarlsdóttirR.VilhjálmssonA. M.KristjánssonB. R. (2021). “Detailed resistivity models of sub-areas of the Hengill geothermal field, Iceland,” in Proceedings world geothermal congress.

  • Google Scholar

• 11

  BessasonB.ÓlafssonE. H.GunnarssonG.FlóvenzÓ. G.JakobsdóttirS. S.BjörnssonS.et al (2012). Verklag vegna örvaðrar skjálftavirkni í jarðhitakerfum. Reykjavik, Iceland: Orkuveita Reykjavíkur.

  • Google Scholar

• 12

  BjarnasonI. T.MenkeW.FlóvenzÓ. G.CaressD. (1993). Tomographic image of the mid‐Atlantic plate boundary in southwestern Iceland. J. Geophys. Res.98 (B4), 6607–6622. 10.1029/92JB02412

  • CrossRef
  • Google Scholar

• 13

  BjörnssonG.HjartarsonA.BöðvarssonG.SteingrimssonB. (2003). “Development of a 3-D geothermal reservoir model for the greater Hengill volcano in SW-Iceland,” in Proceedings of TOUGH symposium 2003 (Berkeley, CA: Lawrence Berkeley National Laboratory), 12–14.

  • Google Scholar

• 14

  BodvarssonR.RognvaldssonS. T.JakobsdottirS. S.SlungaR.StefanssonR. (1996). The SIL data acquisition and monitoring system. Seismol. Res. Lett.67 (5), 35–46. 10.1785/gssrl.67.5.35

  • CrossRef
  • Google Scholar

• 15

  BrandsdóttirB.ParsonsM.WhiteR. S.GuðmundssonO.DrewJ.ThorbjarnadóttirB. (2010). The May 29 2008 earthquake aftershock sequence within the South Iceland Seismic Zone: Fault locations and source parameters of aftershocks. Jokull J. Glaciol. Geol. Soc. Icel.60, 23–46.

  • Google Scholar

• 16

  De LandroG.AmorosoO.RussoG.ZolloA. (2020). 4D travel-time tomography as a tool for tracking fluid-driven medium changes in offshore oil–gas exploitation areas. Energies1322, 5878. 10.3390/en13225878

  • CrossRef
  • Google Scholar

• 17

  De LandroG.SerlengaV.RussoG.AmorosoO.FestaG.BrunoP. P.et al (2017). 3-D ultra-high resolution seismic imaging of shallow solfatara crater in campi flegrei (Italy): New insights on deep hydrothermal fluid circulation processes. Sci. Rep.7, 3412. 10.1038/s41598-017-03604-0

  • CrossRef
  • Google Scholar

• 18

  DuranA. (2021). Numerical sensitivity kernels in elastic media for imaging purposes and seismic tomography in the Hengill geothermal field. Dr. Diss. ETH-Zurich. 10.3929/ethz-b-000518037

  • CrossRef
  • Google Scholar

• 19

  EinarssonP. (2008). Plate boundaries, rifts and transforms in Iceland. Jokull58, 35–58.

  • Google Scholar

• 20

  FlóvenzÓ. G.ÁgústssonK.GuðnasonE. Á.KristjánsdóttirS. (2015). “Reinjection and induced seismicity in geothermal fields in Iceland,” in Proceedings world geothermal congress 2015 (Melbourne, Australia.

  • Google Scholar

• 21

  FoulgerG. R. (1988a). Hengill triple junction, SW Iceland 1. Tectonic structure and the spatial and temporal distribution of local earthquakes. J. Geophys. Res.93 (B11), 13493–13506. 10.1029/JB093iB11p13493

  • CrossRef
  • Google Scholar

• 22

  FoulgerG. R. (1988b). Hengill triple junction, SW Iceland 2. Anomalous earthquake focal mechanisms and implications for process within the geothermal reservoir and at accretionary plate boundaries. J. Geophys. Res.93 (B11), 13507–13523. 10.1029/JB093iB11p13507

  • CrossRef
  • Google Scholar

• 23

  FoulgerG. R. (1984). Seismological studies at the Hengill geothermal area, S W Iceland, 3. England: Univ.of Durham, 13. PhD thesis.

  • Google Scholar

• 24

  FoulgerG. R. (1995). The Hengill geothermal area, Iceland: Variation of temperature gradients deduced from the maximum depth of seismogenesis. J. Volcanol. Geotherm. Res.65 (1-2), 119–133. 10.1016/0377-0273(94)00088-X

  • CrossRef
  • Google Scholar

• 25

  FoulgerG. R.ToomeyD. R. (1989). Structure and evolution of the hengill‐grensdalur volcanic complex, Iceland: Geology, geophysics, and seismic tomography. J. Geophys. Res.94 (B12), 17511–17522. 10.1029/JB094iB12p17511

  • CrossRef
  • Google Scholar

• 26

  GunnarsdóttirS. H.HelgadóttirH. M.SnaebjornsdottirS. (2020). Geco: Geological properties, permeability and porosity of the Nesjavellir high temperature area in relation to the Re-injection of geothermal CO2 and H2S gases. Proc. World Geotherm. Congr.

  • Google Scholar

• 27

  GunnarssonG.AradóttirE. S. P. (2015). The deep roots of geothermal systems in volcanic areas: Boundary conditions and heat sources in reservoir modeling. Transp. Porous Media108, 43–59. 10.1007/s11242-014-0328-1

  • CrossRef
  • Google Scholar

• 28

  GunnarssonG.ArnaldssonA.OddsdóttirA. L. (2011). Model simulations of the Hengill area, southwestern Iceland. Transp. Porous Media90, 3–22. 10.1007/s11242-010-9629-1

  • CrossRef
  • Google Scholar

• 29

  HelgadóttirH. M.SnæbjörnsdóttirS. Ó.NíelssonS.GunnarsdóttirS. H.MatthíasdóttirT.HarðarsonB. S.et al (2010). “Geology and hydrothermal alteration in the reservoir of the Hellisheiði high temperature system, SW-Iceland,” in Proceedings world geothermal congress 2010 bali (Indonesia, 25–29.

  • Google Scholar

• 30

  HersirG. P.BjörnssonG.BjörnssonA. (1990). Eldstöðvar og jarðhiti á Hengilssvæði. Jarðeðlisfræðileg könnun. Orkustofnun OS-90031/JHD-06 93 (in Icelandic).

  • Google Scholar

• 31

  HersirG. P.GuðnasonE. Á.FlóvenzÓ. G. (2022). “Geophysical exploration techniques,” in Comprehensive renewable energy. Editor LetcherT. M.. 2nd edition (Oxford: Elsevier), 7, 26–79. 10.1016/B978-0-12-819727-1.00128-X

  • CrossRef
  • Google Scholar

• 32

  HjörleifsdóttirV.IngvarssonG.RatouisT.GunnarssonG.SnæbjörnsdóttirS. Ó.SigfússonB. (2021a). “Ten years of induced earthquakes in the Húsmúli CO2 injection site, Hellisheiði, Iceland,” in Sixth international conference on engineering Geophysics, virtual, 25–28. 10.1190/iceg2021-027.1

  • CrossRef
  • Google Scholar

• 33

  HjörleifsdóttirV.SnæbjörnsdóttirS. Ó.GunnarssonG.KristjánssonB. R.VogfjordK.JónsdóttirK.et al (2021b). “Induced earthquakes in the Hellisheiði geothermal field, Iceland,” in Proceedings world geothermal congress 2020+1 reykjavik (Iceland.

  • Google Scholar

• 34

  JakobsdóttirS. S. (2008). Seismicity in Iceland: 1994–2007. Jökull58, 75–100.

  • Google Scholar

• 35

  JoussetP.HaberlandC.BauerK.ÁrnasonK. (2011). Hengill geothermal volcanic complex (Iceland) characterized by integrated geophysical observations. Geothermics40 (1), 1–24. 10.1016/j.geothermics.2010.12.008

  • CrossRef
  • Google Scholar

• 36

  JuncuD.ÁrnadóttirT.HooperA.GunnarssonG. (2017). Anthropogenic and natural ground deformation in the Hengill geothermal area, Iceland. J. Geophys. Res. Solid Earth122, 692–709. 10.1002/2016jb013626

  • CrossRef
  • Google Scholar

• 37

  KisslingE.EllsworthW. L.Eberhart‐PhillipsD.KradolferU. (1994). Initial reference models in local earthquake tomography. J. Geophys. Res.99 (B10), 19635–19646. 10.1029/93JB03138

  • CrossRef
  • Google Scholar

• 38

  KisslingE. (1988). Geotomography with local earthquake data. Rev. Geophys.26 (4), 659–698. 10.1029/rg026i004p00659

  • CrossRef
  • Google Scholar

• 39

  KisslingE.KradolferU.MaurerH. (1995). Program VELEST user’s guide-Short Introduction. ETH Zurich: Institute of Geophysics.

  • Google Scholar

• 40

  KristjánsdóttirS.GudnasonE. Á.ÁgústssonK.ÁgústsdóttirTh. (2019). Hverahlíð, Hengill area: Detailed analysis of seismic activity from december 2016 to december 2019, 54. ÍSOR - Iceland GeoSurvey, Reykjavík, Report, ÍSOR-2019/051.

  • Google Scholar

• 41

  LatorreD.VirieuxJ.MonfretT.MonteillerV.VanorioT.GotJ. L.et al (2004). A new seismic tomography of Aigion area (Gulf of Corinth, Greece) from the 1991 data set. Geophys. J. Int.159 (3), 1013–1031. 10.1111/j.1365-246X.2004.02412.x

  • CrossRef
  • Google Scholar

• 42

  LiK. L.AbrilC.GudmundssonO.GudmundssonG. B. (2019). Seismicity of the Hengill area, SW Iceland: Details revealed by catalog relocation and collapsing. J. Volcanol. Geotherm. Res.376, 15–26. 10.1016/j.jvolgeores.2019.03.008

  • CrossRef
  • Google Scholar

• 43

  LinG.WuB. (2018). Seismic velocity structure and characteristics of induced seismicity at the Geysers Geothermal Field, eastern California. Geothermics71, 225–233. 10.1016/j.geothermics.2017.10.003

  • CrossRef
  • Google Scholar

• 44

  LomaxA.VirieuxJ.VolantP.Berge-ThierryC. (2000). “Probabilistic earthquake location in 3D and layered models,” in Advances in seismic event location (Dordrecht: Springer), 101–134.

  • Google Scholar

• 45

  MajerE. L.PetersonE. (2007). The impact of injection on seismicity at the Geysers, California Geothermal Field. Int. J. Rock Mech. Min. Sci.44 (8), 1079–1090. 10.1016/j.ijrmms.2007.07.023

  • CrossRef
  • Google Scholar

• 46

  MavkoG.MukerjiT.DvorkinJ. (2009). The rock physics handbook. Second Edition. New York: Cambridge University Press. Published in the United States of America by.

  • Google Scholar

• 47

  MillerA. D.JulianB. R.FoulgerG. R. (1998). Three-dimensional seismic structure and moment tensors of non-double-couple earthquakes at the Hengill–Grensdalur volcanic complex, Iceland. Geophys. J. Int.133 (2), 309–325. 10.1046/j.1365-246X.1998.00492.x

  • CrossRef
  • Google Scholar

• 48

  NapolitanoF.AmorosoO.La RoccaM.GervasiA.GabrielliS.CapuanoP. (2021). Crustal structure of the seismogenic volume of the 2010–2014 pollino (Italy) seismic sequence from 3D P- and S-wave tomographic images. Front. Earth Sci.9. 10.3389/feart.2021.735340

  • CrossRef
  • Google Scholar

• 49

  PaigeC. C.SaundersM. A. (1982). Lsqr: An algorithm for sparse linear equations and sparse least squares. ACM Trans. Math. Softw.8 (1), 43–71. 10.1145/355984.355989

  • CrossRef
  • Google Scholar

• 50

  PodvinP.LecomteI. (1991). Finite difference computation of traveltimes in very contrasted velocity models: A massively parallel approach and its associated tools. Geophys. J. Int.105 (1), 271–284. 10.1111/j.1365-246X.1991.tb03461.x

  • CrossRef
  • Google Scholar

• 51

  RögnvaldssonS. Th.GudmundssonG.ÁgústssonK.StefánssonR.JakobsdóttirS. (1996). “Recent seismicity near the Hengill triple-junction, SW Iceland,” in Proceedings of the XXV general assembly of the eur. Seismol. Comm., Icelandic meteorological Office, Reykjavík, Iceland. Editor ThorkelsonB., 461–466. Seismology in Europe.

  • Google Scholar

• 52

  SæmundssonK. (1962). Das alter der Nesja-lava (SW-island). Neues Jahrb. Geol. Paleontol. Monatsh.12, 650.

  • Google Scholar

• 53

  SæmundssonK.SaemundssonK. (1992). Geology of the thingvallavatn area. Oikos64 (1/2), 40–68. 10.2307/3545042

  • CrossRef
  • Google Scholar

• 54

  SæmundssonK.SigurgeirssonM. Á.HjartarsonÁ.KaldalI.KristinssonS. G.VíkingssonS. (2016). Geological map of SW Iceland, 1:100 000 (2. edition). Reykjavík: Iceland GeoSurvey.

  • Google Scholar

• 55

  SæmundssonK. (1967). Vulkanismus und Tektonik des Hengill-Gebietes in Sudwest-Island. Acta Nat. Isl.II-7, 195. (in German).

  • Google Scholar

• 56

  Sánchez-PastorP.ObermannA.ReinschT.ÁgústsdóttirT.GunnarssonG.TómasdóttirS.et al (2021). Imaging high-temperature geothermal reservoirs with ambient seismic noise tomography, a case study of the Hengill geothermal field, SW Iceland. Geothermics96, 102207. 10.1016/j.geothermics.2021.102207

  • CrossRef
  • Google Scholar

• 57

  SerlengaV.de LorenzoS.RussoG.AmorosoO.GaramboisS.VirieuxJ.et al (2016). A three‐dimensional Q_P imaging of the shallowest subsurface of Campi Flegrei offshore caldera, southern Italy. Geophys. Res. Lett.43, 218. 10.1002/2016GL071140

  • CrossRef
  • Google Scholar

• 58

  SigmundssonF.EinarssonP.RögnvaldssonS. T.FoulgerG. R.HodgkinsonK. M.ThorbergssonG. (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 (B7), 15151–15161. 10.1029/97JB00892

  • CrossRef
  • Google Scholar

• 59

  SintonJ.GrönvoldK.SæmundssonK. (2005). Postglacial eruptive history of the western volcanic zone, Iceland. Geochem. Geophys. Geosyst.6, Q12009. 10.1029/2005GC001021

  • CrossRef
  • Google Scholar

• 60

  StefánssonR.BöđvarssonR.SlungaR.EinarssonP.JakobsdóttirS.BungumH.et al (1993). Earthquake prediction research in the south Iceland seismic zone and the SIL project. Bull. Seismol. Soc. Am.83 (3), 696–716.

  • Google Scholar

• 61

  TryggvasonA.RögnvaldssonS. T.FlóvenzO. G. (2002). Three-dimensional imaging of the P-and S-wave velocity structure and earthquake locations beneath Southwest Iceland. Geophys. J. Int.151 (3), 848–866. 10.1046/j.1365-246X.2002.01812.x

  • CrossRef
  • Google Scholar

• 62

  UssherG.HarveyC.JohnstoneR.AndersonE. (2000). “Understanding the resistivities observed in geothermal systems,” in Proceedings world geothermal congressKyushu, 1915–1920.

  • Google Scholar

• 63

  VilhjálmssonA. M.BenediktsdóttirÁ.ÁrnasonK. (2019). Nesjavellir og Þverárdalur Einvíð túlkun viðnámsgagna og þrívíð fyrir Nesjavallasvæðið. ÍSOR report, in Icelandic, ÍSOR-2019/062,

  • Google Scholar

• 64

  ViolayM.GibertB.MainpriceD.EvansB.DautriaJ. M.AzaisP.et al (2012). An experimental study of the brittle-ductile transition of basalt at oceanic crust pressure and temperature conditions. J. Geophys. Res.117. 10.1029/2011JB008884

  • CrossRef
  • Google Scholar

• 65

  WadatiK.OkiS. (1933). On the travel time of earthquake waves.(Part II). J. Meteorological Soc. Jpn.11 (1), 14–28. 10.2151/jmsj1923.11.1_14

  • CrossRef
  • Google Scholar

• 66

  WagnerF. (2019). “Toward fully automatic earthquake detection and processing for tomography in the Hengill area,” in Digital comprehensive summaries of uppsala dissertations from the faculty of science and technology 1774. Uppsala: Acta Universitatis Upsaliensis, 53. ISBN 978-91-513-0574-5. Available at: http://uu.diva-portal.org/smash/get/diva2:1286849/FULLTEXT01.pdf.

  • Google Scholar

• 67

  ZierenbergR. A.SchiffmanP.BarfodG. H.LesherC. E.MarksN. E.LowensternJ. B.et al (2013). Composition and origin of rhyolite melt intersected by drilling in the Krafla geothermal field, Iceland. Contrib. Mineral. Petrol.165 (2), 327–347. 10.1007/s00410-012-0811-z

  • CrossRef
  • Google Scholar