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

Amoroso, Ortensia; Napolitano, Ferdinando; Hersir, Gylfi Pall; Agustsdottir, Thorbjorg; Convertito, Vincenzo; De Matteis, Raffaella; Gunnarsdóttir, Sveinborg Hlíf; Hjörleifsdóttir, Vala; Capuano, Paolo

doi:10.3389/feart.2022.994280

ORIGINAL RESEARCH article

Front. Earth Sci., 12 October 2022

Sec. Solid Earth Geophysics

Volume 10 - 2022 | https://doi.org/10.3389/feart.2022.994280

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

Ortensia Amoroso¹*

Ferdinando Napolitano¹

Gylfi Pall Hersir^2,3

Thorbjorg Agustsdottir²

Vincenzo Convertito⁴

Raffaella De Matteis⁵

Sveinborg Hlíf Gunnarsdóttir²

Vala Hjörleifsdóttir⁶

Paolo Capuano¹

¹Dipartimento di Fisica “E.R.Caianiello”, Università Degli Studi Di Salerno, Fisciano, SA, Italy
²Iceland GeoSurvey (ÍSOR), Reykjavík, Iceland
³Independent Researcher, Reykjavík, Iceland
⁴Istituto Nazionale di Geofisica e Vulcanologia, Osservatorio Vesuviano, Napoli, Italy
⁵Dipartimento di Scienze e Tecnologie, Università Degli Studi Del Sannio, Benevento, Italy
⁶Reykjavík Energy, Reykjavík, Iceland

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 $v_{P}$ , $v_{S}$ , and $\frac{v_{P}}{v_{S}}$ 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 $v_{P}$ and $v_{S}$ , and high $\frac{v_{P}}{v_{S}}$ 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 $\frac{v_{P}}{v_{S}}$ 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 $\frac{v_{P}}{v_{S}}$ 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 $\frac{v_{P}}{v_{S}}$ 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 $\frac{v_{P}}{v_{S}}$ 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 $\frac{v_{P}}{v_{S}}$ 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

FIGURE 1. The Hengill study area: seismic stations and earthquake locations. (A) shows the location of the study area (red box) and the volcanic centers in SW-Iceland. (B) Seismic events extracted from ÍSOR’s catalog (2016–2020) and relocated using the 3D velocity model based on this study. Blue, green, magenta, and red diamonds represent the seismic stations of the ÍSOR/ON network, the COSEISMIQ and S4CE project, and the SIL seismic network deployed in the study area, respectively. AA′ and BB′ profiles, shown as white dotted lines, are discussed in the results paragraph. Relocated earthquakes (2016–2020) are color coded by depth. Solid blue lines are the three volcanic centers of the region with their respective names in white on a blue background. The white square marked as “Ne” is the Nesjavellir power plant. Black lines are the roads. Red lines represent the eruptive fissure and craters, while the fractures are shown using gray lines (Sæmundsson et al., 2016). (C) is a zoomed-in of the part inside the gray box of the figure to the left.

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 ( $v_{P}$ ), S-wave velocity ( $v_{S}$ ) and of the $\frac{v_{P}}{v_{S}}$ 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 $v_{P}$ , $v_{S}$ , and $\frac{v_{P}}{v_{S}}$ 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 $\frac{v_{P}}{v_{S}}$ ratio equal to $1.77 \pm 0.08$ (for the extent of Figure 1B). This value is fully consistent with the $\frac{v_{P}}{v_{S}}$ ratio computed on a larger scale in the wider Hengill region: Foulger (1995) and Miller et al. (1998) found $\frac{v_{P}}{v_{S}} = 1.77$ , Tryggvason et al. (2002) obtained $\frac{v_{P}}{v_{S}} = 1.78$ , Jousset et al. (2011) found $\frac{v_{P}}{v_{S}} = 1.75 \pm 0.12$ .

FIGURE 2

FIGURE 2. Wadati diagram. 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 the 6,906 events recorded at the 14 available stations. The best fit line corresponds to $\frac{v_{P}}{v_{S}} = 1.77 \pm 0.08 .$

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 $\frac{v_{P}}{v_{S}}$ 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

FIGURE 3. 1D-velocity model analysis. (A) Velocity models available in literature for the Hengill area (Stefánsson et al., 1993; Tryggvason et al., 2002; Jousset et al., 2011) used as starting models in VELEST. A fourth synthetic gradient model has been used to test the stability of the results. (B) The 1D velocity models (whole lines) obtained using 473 events from the study area. A final 1D velocity model has been obtained (dashed black line) averaging the velocity in each layer up to a depth of 10 km. This model has been used as input in VELEST. (C) Minimum P-wave 1D velocity model specific for the study area, the Nesjavellir geothermal area and surroundings.

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 $\frac{v_{P}}{v_{S}}$ 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, $v_{P}$ , $v_{S}$ , 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).

TABLE 1

TABLE 1. Minimum 1D P-wave velocity model obtained by VELEST (see Section 3.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 $v_{P}$ and the $v_{S}$ models are the direct output of the tomographic code, the $\frac{v_{P}}{v_{S}}$ model is derived by dividing the values of the $v_{P}$ and $v_{S}$ 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 $v_{P}$ and the $v_{S}$ 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 $v_{P}$ anomaly (labeled as PL in Figure 4), initially emerges in the NW part of the map and extends down to 4 km depth. The $v_{S}$ 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 $\frac{v_{P}}{v_{S}}$ model, the well-resolved area features values less than 1.77 down to 2.5 km depth (Figure 6).

FIGURE 4

FIGURE 4. 3D tomographic model: $v_{P}$ horizontal slices, depth below sea level. Black dots are earthquakes within 1 km around each slice from the 2,358 high quality events used for the tomography. The regions not covered by the ray paths are white. The dark gray contour line delimitates the well-resolved area inferred by the Derivative Weight Sum (DWS). The uppermost slide to the left (z = 0 km) shows for orientation, the seismic stations, roads in black, fissures, and craters in red, fractures in gray and blue lines are the three volcanic centers of the region. The Nesjavellir power plant is located in the white square, labeled Ne. PL is a low $v_{P}$ anomaly, PH is a high $v_{P}$ anomaly (discussed in Section 4,Section5,).

FIGURE 5

FIGURE 5. 3D tomographic model: $v_{S}$ horizontal slices. Black dots are earthquakes within 1 km around each slice from the 2,358 high quality events used for the tomography. The regions not covered by the ray paths are shown in white. The dark gray contour line delimitates the well-resolved area inferred by the Derivative Weight Sum (DWS). The uppermost slide to the left (z = 0 km) shows for orientation the seismic stations, roads in black, fissures and craters in red, fractures in gray and blue lines are the three volcanic centers of the region. The Nesjavellir power plant is located in the white square, labeled Ne. SL is a low $v_{S}$ anomaly (discussed in Section 4,Section5,).

FIGURE 6

FIGURE 6. 3D tomographic model: $\frac{v_{P}}{v_{S}}$ horizontal slices. Black dots are earthquakes within 1 km around each slice from the 2,358 high quality events used for the tomography. The regions not covered by the ray paths are white. The dark gray contour line delimitates the well-resolved area inferred by the Derivative Weight Sum (DWS). The uppermost slide to the left (z = 0 km) shows for orientation the seismic stations, roads in black fissures and craters in red, fractures in gray and blue lines are the three volcanic centers of the region. The Nesjavellir power plant is located in the yellow square, labeled Ne. PSL is a low $\frac{v_{P}}{v_{S}}$ anomaly, PSH is a high $\frac{v_{P}}{v_{S}}$ anomaly (discussed in Section 4,Section5,).

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 $v_{S}$ does not change significantly (Figure 5). The anomalies in the $v_{P}$ model produce a notable complex pattern in the $\frac{v_{P}}{v_{S}}$ model down to 4 km. In particular, at 3 km depth, we note two localized areas of high $\frac{v_{P}}{v_{S}}$ ratio (PSH) and a low $\frac{v_{P}}{v_{S}}$ 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 $\frac{v_{P}}{v_{S}}$ 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 $v_{P}$ , $v_{S}$ and $\frac{v_{P}}{v_{S}}$ ratio in depth view are displayed on Figure 7.

FIGURE 7

FIGURE 7. Tomographic cross-sections across (A) A′ and along (B) B′ the Nesjavellir valley. Profile locations are given on Figure 1. Black dots represent earthquakes located within 250 m on each side from the vertical cross-section (A) A′ (left panels) and within 500 m on each side from the vertical cross-section (B) B′ (right panels). From top to bottom: P-wave velocity ( $v_{P}),$ S-wave velocity ( $v_{S})$ and $\frac{v_{P}}{v_{S}}$ ratio. The dark gray contour line delimitates the well-resolved area inferred by the Derivative Weight Sum (DWS). The vertical black line represents the intersection of the two cross-section. Nesjavellir power plant is marked with a black square. The production area is indicated with two vertical dashed lines, the well trajectories are shown on Figure 9.

The distribution of the relocated seismicity using the new detailed 3D $v_{P}$ and $v_{S}$ 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 $v_{P}$ value (PL), low $v_{S}$ value (SL), and low $\frac{v_{P}}{v_{S}}$ ratio (PSL) (Figures 4–7). The seismicity at depths greater than about 3 km is located in a high $v_{P}$ (PH) and high $\frac{v_{P}}{v_{S}}$ ratios (PSH) area.

FIGURE 8

FIGURE 8. Relocated seismicity used in the 3D velocity models. Map and cross-section of the relocated seismicity of the Nesjavellir geothermal field and surroundings. The inverted orange triangles are the well heads and the diamonds show the seismic stations. Orange lines in the cross-sections are borehole trajectories. Red lines represent the eruptive fissure and craters, while the fractures are shown using gray lines (Sæmundsson et al., 2016). Roads are blue and the Nesjavellir power plant is in the yellow square, labeled Ne. The deep injection well NJ-18 is labeled, and its trajectory shown in black. The shallow wells northeast of NJ-18, are also injection wells (note some are shallower than 400 m and do, therefore, not plot clearly on the figure). Earthquakes are color coded by depth.

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 $v_{P}$ , $v_{S}$ , and $\frac{v_{P}}{v_{S}}$ 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 $v_{P}$ (4.0 km/s, Figure 4) and the still lower $v_{S}$ (2.3 km/s, Figure 5), which results in high $\frac{v_{P}}{v_{S}}$ ratio values (1.9, Figure 6), suggest the presence of porous material characterized by the presence of fluid. A similar $\frac{v_{P}}{v_{S}}$ 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 $\frac{v_{P}}{v_{S}}$ ratios (1.64–1.70), due to low $v_{P}$ (5.2 km/s) and low-to-average $v_{S}$ (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 $\frac{v_{P}}{v_{S}}$ 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 $v_{P}$ (PH) and low-to-average $v_{S}$ values, resulting in a high $\frac{v_{P}}{v_{S}}$ ratio anomaly (labeled PSH). The highest value of $\frac{v_{P}}{v_{S}}$ is larger than 1.9 and is reached at 4.5 km (the retrieved $\frac{v_{P}}{v_{S}}$ 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

FIGURE 9. Cross-sections (A) A′ and (B) B′ (locations are given on Figure 1). Earthquake’s hypocenters are shown on all four figures as black dots [500 m to each side for (A) A′ and 250 m for (B) B′] and contour lines for the $\frac{v_{P}}{v_{S}}$ ratio. Left panel: Resistivity model (Vilhjálmsson et al., 2019); Right panel: Formation temperature (Gunnarsson et al., 2011). Well heads and trajectories at 300 m distance from each profile have been added to the profiles. Nesjavellir power plant is marked with a black square. Red box on the left panel shows the extent of the temperature data on the right panel in relation to the resistivity data.

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, $\frac{v_{P}}{v_{S}}$ 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 $\frac{v_{P}}{v_{S}}$ 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 $\frac{v_{P}}{v_{S}}$ 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 $\frac{v_{P}}{v_{S}}$ 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 $\frac{v_{P}}{v_{S}}$ 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 $\frac{v_{P}}{v_{S}}$ 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 $\frac{v_{P}}{v_{S}}$ ratio is due to an increase in $v_{P}$ that dominates with respect to the decrease in $v_{S}$ (Figure 7). This anomalously high $\frac{v_{P}}{v_{S}}$ ratio (>1.9) correlates well with the deep-seated low resistivity on both profiles as shown on Figure 9 and discussed above. A high $\frac{v_{P}}{v_{S}}$ 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 $\frac{v_{P}}{v_{S}}$ 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 $\frac{v_{P}}{v_{S}}$ 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 $\frac{v_{P}}{v_{S}}$ ratio. It is most interesting, that the NNE-SSW trending seismicity is located outside, on the border, of the high $\frac{v_{P}}{v_{S}}$ 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

FIGURE 10. A horizontal slice at 4 km depth through the resistivity model (Vilhjálmsson et al., 2019) and the tomographic model showing $\frac{v_{P}}{v_{S}}$ ratio as contours. Black dots are earthquakes, at 500 m distance from each side of the slide.

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 $\frac{v_{P}}{v_{S}}$ 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 $\frac{v_{P}}{v_{S}}$ 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 $\frac{v_{P}}{v_{S}}$ 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 $\frac{v_{P}}{v_{S}}$ 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.

Data availability statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author contributions

OA and PC conceived the work. FN, OA, VC, and RD performed the 1D velocity model, earthquake relocation and 3D tomographic analysis. GH, TÁ, SG, OA, and FN worked on the joint interpretation of 3D seismic tomography, resistivity and well temperature data. VH and TÁ provided the seismic catalogue data and related metadata. OA, FN, GH, and TÁ worked on draft manuscript preparation. All authors discussed the results and contributed to the final manuscript.

Funding

This work has been supported by the S4CE (“Science for Clean Energy”) project, funded by the European Union’s Horizon 2020—R&I Framework Programme, under grant agreement No 764810 and by PRIN-2017 MATISSE project, No 20177EPPN2, funded by the Italian Ministry of Education and Research.

Acknowledgments

We thank the Geothermica ERA-net funded COSEISMIQ project (project No. 170167-4401) for data from six of their stations, and Icelandic Meteorological Office (IMO) for data from one station. We thank ÍSOR and ON Power for the use of their extensive dataset on manually refined phase-picks.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2022.994280/full#supplementary-material

References

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

Google Scholar

Ágústsson, K., and 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

Ágústsson, K., Kristjánsdóttir, S., Flóvenz, Ó. G., and 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

Amoroso, O., Ascione, A., Mazzoli, S., Virieux, J., and Zollo, A. (2014). Seismic imaging of a fluid storage in the actively extending Apennine mountain belt, southern Italy. Geophys. Res. Lett. 41, 3802–3809. doi:10.1002/2014GL060070