3D S-wave velocity imaging of Reykjanes Peninsula high-enthalpy geothermal fields with ambient-noise tomography

Tomographic imaging based on ambient seismic noise measurements has shown to be a powerful tool, especially in areas like Iceland, where the microseism illumination is excellent. In this paper, we produce a 3D S-wave tomographic image over the western Reykjanes Peninsula high-enthalpy geothermal ﬁelds and evaluate the reliability of the tomographic results for different resolutions through simulated and real data. We use 30 broadband stations operating for approximately one-and-a-half year and apply ambient noise seismic interferometry for each station-pair. This results in empirical Green’s functions in which especially the ballistic surface waves (BSW) are well resolved. The retrieved BSW exhibit a high signal-to-noise ratio between 0.1 and 0.5 Hz, and the beamforming analysis indicates an apparent surface-wave velocity of 3 km/s over a broad azimuthal range. For the tomographic inversion, we invert the estimated phase velocities between all station pairs to frequency-dependent phase velocity maps in four different resolutions (1, 2, 3, and 4 km) using a Tikhonov regularisation. With the estimated regularisation parameter per frequency per resolution, we invert simulated data for checkerboard sensitivity tests per frequency for different combinations of velocity anomaly sizes and resolutions. Finally, after the inversion to depth, we detect S-wave velocity anomalies with variations between − 15% and 15% with reference to an estimated average velocity using 1 km and 3 km of lateral resolutions and 1 km of vertical resolution. This study shows the potential of ambient noise tomography as complementary seismological tool for reservoir characterization.


Introduction
The Western Reykjanes Peninsula (WRP), located at the southwestern tip of Iceland (Fig. 1) is a transitional zone between a tectonic spreading centre (the extension over Iceland of the Mid-Atlantic Ridge) and a transform zone (Pálmason and Saemundsson, 1974;Sigmundsson et al., 2018). At this transitional zone, an extensional component with northeast-southwest trending normal faulting, and a transform component with strike-slip faulting oriented north-south can explain observed displacement components (Klein et al., 1973;Einarsson, 2008). This tectonic setting facilitates the intrusion of magma in dikes along the NE-SW trending area of densely spaced fissures and faults. The presence of magmatic intru-sions and fault swarms, in turn, promotes cooler "fresh" ground water inflow and thermal upflow (Franzson, 1987). The resultant geothermal manifestations have made the WRP subject of geothermal power production for over 30 years.
Currently, attempts are made to explore deep geothermal sources, the supercritical part of hydrothermal systems, which offers different potential for harnessing geothermal energy. One of the most noticeable advantages of supercritical systems is that the increase of pressure (>221 bars or 22.06 MPa) and temperature (>374 • C higher for sea water) can potentially induce an estimated power output ∼10 times higher than traditional Icelandic geothermal wells (Friðleifsson and Albertsson, 2000;Albertsson et al., 2003). In Iceland, while customary geothermal systems usually reach down to 3 km depth with steam temperatures between 290 • C and 320 • C, predictions by Friðleifsson et al. (2014a) report that deep wells reaching depths between 4 km and 5 km may have temperatures between 400 • C and 600 • C. These results were tested  The green triangles denote the onshore broadband stations. The red dots identify wellsheads in Reykjanes geothermal fields, and the black lines the identified faults and fractures in the vicinity of the geothermal production areas. The squares locate high temperature areas (mapped by the Icelandic GeoSurvey ISOR (Guðnason et al., 2014)) from west to east, Reykjanes (R), Eldvörp (E), Svartsengi (S), of which the red squares indicate the location of the two existing power plants. b) Icelandic coastal boundary with Reykjanes peninsula within the red rectangle. Red dashed line locates the plate boundary and black lines locate the transform fault system. The neovolcanic zones are filled in light shaded grey along (Einarsson and Saemundsson, 1987;Erlendsson and Einarsson, 1996;Einarsson et al., 2002;Einarsson, 2008;Sigmundsson et al., 2018). and partially confirmed with IDDP2 in Reykjanes (drilled between 2016 and 2017) with measured temperature 427 • C and fluid pressure of 340 bars at ∼4.5 km depth (Friðleifsson et al., 2017a). In 2009, within the Icelandic Deep Drilling Project (IDDP), the first IDDP well (IDDP-1) drilled in Iceland in the Krafla area encountered rhyolitic magma at a depth of 2104 m (Hólmgeirsson et al., 2010). Magnetotelluric (MT) electromagnetic surveys performed at Krafla wrongly estimated a magma source at 4.5 km deep (Elders et al., 2011;Gasperikova et al., 2015), and the drilling encountered shallower magma. Even though MT surveys are the preferred geophysical measurement procedure to locate the water reservoirs required for the closed loop of geothermal production (e.g. Newman et al., 2008;Árnason et al., 2010;Hersir et al., 2018) and have been utilized frequently in Iceland (e.g. Hersir et al., 1984;Eysteinsson and Hermance, 1985;Oskooi et al., 2005), the magma pocket that was drilled into at the Krafla IDDP-1 well location was possibly below the level of resolution of the area surveyed using MT (Árnason et al., 2007). The IDDP-2 well drilled in December 2016 at the WRP saline geothermal system in southwest Iceland successfully reached approximately 4.5 km (vertical) depth (Friðleifsson et al., 2017b;Elders et al., 2014;Friðleifsson et al., 2018). For the location of the well, the operators could consider results from resistivity measurements (Karlsdóttir et al., 2018;Friðleifsson et al., 2014b;Darnet et al., 2018), as well as travel time seismic tomography (Jousset et al., 2017). The seismic tomography was estimated from a seismic campaign developed under the European-funded program Integrated Methods for Advanced Geothermal Exploration (IMAGE). The Jousset et al. (2017) tomographic study confirms the previous results of Bjarnason et al. (1993) and of Tryggvason et al. (2002) over the same area with enhanced details around well locations. The authors interpret the low-ratio anomaly of compressional-over shear-velocity as being due to the absence of a sizeable magmatic body at the tip of Reykjanes Peninsula, which was confirmed by Friðleifsson et al. (2018).
The unexpected drilling into a magma source in Krafla in 2009 highlighted the need to explore high-resolution imaging techniques as a complement to current measurement methods and to improve the existing ones. Seismic tomographic techniques and recent advances using ambient noise-based methodologies can play a role in assessing the necessary depth and resolution information and in constraining other geophysical estimations. In this regard, ambient noise seismic interferometry (ANSI) techniques can offer additional advantages by avoiding the cost of active seismic methods (therefore, making ANSI techniques economically more attractive), and circumventing limitations such as the limited number of earthquakes and irregular earthquake distribution. The straightforward data acquisition and the theoretical concepts, extended from 1D media (Claerbout, 1968) to arbitrarily heterogeneous 3D media by Wapenaar (2004), make ANSI attractive for tomography applications. The ANSI concept relies on a virtual source that is generated at the location of one of the two receivers by cross-correlation and summation of (noise) recordings from surrounding ambient noise sources. The tomographic results are subsequently derived from Rayleigh (or Love) waves retrieved between the virtual sources and the receivers. The number of applications of ambient noise tomography (ANT) studies has increased in recent years, especially in Iceland Benediktsdóttir et al., 2017;Jeddi et al., 2017;Green et al., 2017;Martins et al., 2019). Along with the direct advantage of characterising the subsurface within the depth range of Icelandic geothermal operations, ANT can contribute to constraining other geophysical measurements or interpretations that have previously been acquired over the same area and vice versa. ANT can be used to improve a subsurface image with complementary seismic studies, Blanck et al., 2019;Jousset et al., 2017).
In this study, we derive a 3D S-wave velocity tomographic image of the WRP's subsurface by applying ANT to the seismic survey deployed under the IMAGE project framework. On top of assessing the reliability of the retrieved dispersion curves, we devote particular attention to the model resolution given the deployed network configuration and to obtain results that allow to constrain other geophysical measurements.

Data and methodology
Within the IMAGE project framework, the German Research Center for Geosciences (GFZ Potsdam) and Iceland Geosurvey (ISOR) deployed 54 broadband seismometers on and around the  (Blanck et al., 2019, Jousset et al., 2017. The onshore instruments were operating from March 2014, and the OBS were placed in August 2014. All the equipment was collected in August 2015. The OBS deployed during the IMAGE project were not used in this study due to a phase shift in the instruments . We used the 30 onshore (Table 1 3 components seismometers (20 broadband Trillium compact sensors and 10 short-period MarkSensors) with a corner frequency as low as 0.005 Hz and a sampling rate of 200 Hz for a duration of almost one year and five months. We only use the vertical-component displacements and to reduce computation time, we down-sampled the records to 25 samples per second (Nyquist of 12.5 Hz).
The applied methodology follows the processing approach of Martins et al. (2019) integrally and the processing chain is depicted in Fig. 2 with a division between data-processing and inversion schemes (identified in Fig. 2 by the orange and green dashed squares, respectively). We divide the processing chain into two steps: data processing of the retrieved surface-waves and the inversion procedure. The pre-processing includes deconvolution of the instrument responses, spectral whitening, temporal averaging (Bensen et al., 2007) and frequency filtering. The next step is Empirical Green's Functions (EGF) retrieval with seismic interferometry. The tomographic inversion scheme makes use of the estimated Ballistic Surface Waves (BSW) arrival-times from the EGF to estimate frequency dependent spatial velocity anomalies. Finally, we estimate the depth-dependent 3D S-wave velocity anomalies.
In the rest of this section, we describe the data pre-processing, EGF retrieval, BSW arrival time picking, tomography and inversion to S-wave-velocity. Between these steps, we perform quality checks to ensure that: 1) The signal-to-noise ratio is high enough to allow acceptable surface-wave arrival time estimation. 2) The illumination is sufficiently uniform (taking into account that a lack of causal illumination can be compensated by the acausal part of the cross-correlated signal). 3) The estimated time picks are consistent. 4) Smooth velocity variations between the tomographic frequency dependent results while inverted independently, and 5) We are choosing the most appropriate resolution.

BSW from cross-correlations
The instrument response is removed by complex deconvolution, after which we apply a spectral domain normalisation (i.e. whitening) (Bensen et al., 2007). We extract coherent EGFs between 0.1 and 0.5 Hz after a spectrogram examination of the ambient noise. This bandwidth is dominated by the microseisms. In the raw spectrogram (Fig. 3 a) we identify a higher time-dependent (seasonal) power spectral density (PSD) band around 0.2 Hz. In the same figure, the sharp lines covering the entire frequency spectrum with large PSD mark the occurrence of earthquakes. The seasonal effect with increased PSD from mid-August to April (top Fig. 3) occurs due to the higher number of ocean storms and larger waves in the autumn and winter (Ardhuin et al., 2011). From the spectrograms we see that there is sufficient energy up to 0.8 Hz.
We compute the cross-correlations per hour and stack the computed cross-correlations using approximately one year and five months of recorded seismic data. Fig. 4 shows the resulting EGF's obtained between the 435 unique station pairs, using data between 0.1 and 0.5 Hz. The extracted EGF's are highly coherent and show a non-symmetrical (in amplitude) 'V' shape of the BSW arrivals indicating, as expected, a non-isotropic azimuthal distribution of ambient noise sources (Froment et al., 2010).
Strictly speaking, the retrieved BSWs only coincide with the surface wave part of the Green's function and its time-reversed version under the condition that (i) the receiver pairs are illuminated uniformly from all angles (Wapenaar and Fokkema, 2006), (ii) a single surface wave mode dominates the recorded ambient vibrations (Halliday and Curtis, 2008), and (iii) the medium is lossless. In practice, and therefore also in our case, these conditions are not fulfilled, leading to deviations of the extracted surface wave velocities from the true surface wave velocities, e.g. (Tsai, 2009;Froment et al., 2010). If the illumination pattern is sufficiently uniform over an angle-range of at least 180 degrees, it suffices to use only the causal or acausal BSW. It has been shown that this still allows us to obtain meaningful tomographic (surface wave) images (e.g. Shapiro et al., 2005;De Ridder et al., 2015). Fig. 5 a), b) and c) show three of the filtered EGFs around 0.18 Hz, 0.28 Hz and 0.38 Hz, respectively, and cross-correlations filtered with a small frequency window around specified frequency values (+/-0.01 Hz). The coherent 'V' shaped wave pattern of the BSW indicates that it is well possible to pick arrival times. From Fig. 5, we can recognise that lower frequencies travel faster. Moreover, it can be seen that at a short interstation distances there is interference between the BSW at causal and acausal times.
We estimate the azimuthal directions of the illumination, using a beamforming analysis applied to the cross-correlation results (Ruigrok et al., 2017).

Surface-wave phase velocity picking
The high coherence of the BSW between adjacent inter-station distances (Fig. 4) allows us to estimate well the timing at local maxima of the correlation waveform (considered to be the correct phase cycle). To obtain individual (i.e., per station pair) phase-velocity estimates, we first extract an average phase velocity dispersion curve for the fundamental mode Rayleigh wave (c(f), where f is frequency) in the frequency domain using the MASW (multichannel analysis of surface waves) algorithm (Park et al., 1998. The average dispersion curve as a function of frequency (Fig. 6) serves as further guidance to avoid phase cycle jumps in the phase picking of individual frequencies (f i ) with a short interval around f i (where f i ∈ [0.1, 0.5] Hz). Realistic velocities are only estimated for frequencies higher than 0.14 Hz. Fig. 6 shows both MASW picking and the average phase velocity dispersion curve (left and right panel, respectively). On top of the average phase velocity, we plot the resulting phase velocity if different cycles were to be picked (green, red and black asterisks Fig. 6 right side). Selecting a correct cycle is straightforward with    the MASW velocity reference. For one cycle velocities are obtained that are close to the MASW velocities. Picking at the maxima of other cycles yields unrealistic velocities (Fig. 6 right). We impose a threshold to withdraw station pairs with distances where too much destructive interference (or low SNR) occurs, with consequent loss of waveform coherence. Higher coherent cross-correlations allow an accurate time-picking arrival.
the same frequency-dependent thresholds as Martins et al. (2019) R i,min = 2 3 i,0 and R i,max = 2.8 i,0 to define the minimum and the maximum admissible inter-station distances, respectively. The "0" subscript refers to the reference (average) phase velocity obtained from MASW. For each station combination, there is (potentially) a causal and an acausal BSW retrieval. We use the BSW that has the highest amplitude. Fig. 7 shows the picks in the time domain. As an additional evaluation, the figure also shows the selected outliers, the values outside the 2 deviation from the mean.
With the estimated time-picks for every f i between 0.1 Hz and 0.5 Hz in steps of 0.02 Hz, we calculate the azimuthal variation of the phase velocity. Azimuthal variations per frequency can be seen as a quality check to detect non-smoothness of time-picks and in cases of success, together with sufficient illumination can also indicate velocity anisotropy. We observe higher azimuthal velocities between 40 • and 80 • , approximately in a northeast-southwest direction, and lower velocities between 100 • and 200 • (Fig. 8). The higher velocities may be correlated with the orientation of the main faults and ridge that cut the WRP (Fig. 1), but further investigation would be required to verify this observation. Azimuthal variations are consistent for all frequencies and have a smooth differential from lower to higher frequencies. These also indicate reliable time  Xia et al. (1999) and Haney and Tsai (2015) theoretical relations and the maximum penetration depth also as function of frequency.
picks (the picks are made independently for each f i ). As a reference for frequency-depth correspondence, we add to Fig. 8, the depth of maximum sensitivity of the Rayleigh waves, as well as the maximum penetration depth, both as a function of frequency.

Tomography
Seismic tomography is an inversion problem that aims to find a slowness field from travel-time observations in three dimensions. The linear forward problem can be written as: where d the data vector (or observations) containing the picked frequency-dependent travel times; m the model vector describing the (unknown) frequency-dependent slowness values of each grid cell; G is the operator matrix containing the ray path lengths of each ray path (first dimension) through each grid cell (second dimension) and; e is a term expressing the noise. We use the same methodology as described in Martins et al. (2019), a first-order Tikhonov regularisation (Tikhonov, 1963) to regularise this inversion problem, which minimises the double objective function: The term ||m|| 2 denotes the norm of the model, multiplied by the so-called regularisation parameter i . The added objective of minimising the norm of the estimated parameters avoids overfitting of the data by enforcing smoothing: where I is the identity matrix and where the superscript T denotes the transpose of a matrix. We use the cross-validation methodology (Golub et al., 1979) to choose a regularisation parameter ( ) per frequency value (f i ). For each frequency, we start by removing one observation tt j where tt represents the travel-time between two stations. Then we use Eq.
(3) to fit a solution without the removed observation (m j ). We predict the neglected travel time using the resolved m j and compare it with the dropped observation itself: where G i is lacking the row associated with station couple j. We perform this for all observations and sum the squares of the residuals to check the robustness of i . We repeat this for multiple and we select the i that minimises P = 1 where n is the number of observations, i.e., the number of station couples for which the phase velocity was estimated.
The tomographic linear inversion is repeated for different frequencies to produce phase-velocity maps. Frequency can be approximated to depth using the theoretical relations in Xia et al. (1999) and Haney and Tsai (2015) (Fig. 8 right). However, for a more accurate frequency-depth conversion it is advised to perform a second inversion using the Rayleigh wave sensitivity kernels (Zhou et al., 2004). This method is described in Section 2.5. As the ray paths per frequency are inverted independently, we also estimate a different regularisation parameter per frequency value, i and f i ∈ [0.18, 0.44] with 0.02 Hz steps.

Checkerboard resolution
An adequate model resolution can help to identify subsurface anomalies' geometries, which is relevant for subsurface characterisation and geothermal purposes. We use checkerboard forward modelling to test the ability of the seismic network geometric coverage to reproduce a simulated checker for each of the tested resolutions versus anomaly size. We test combinations of spatial resolutions (ranging from 1 to 4 km), and size of simulated perturbations (with perturbation sizes ranging from 2 to 6 km). As the number and spatial distribution of the ray paths changes with frequency, we reproduce a checkerboard for each frequency f i . Similarly to the procedure for the field-data inversion, the checkerboard inversion is done independently for each f i using the estimated Tikhonov regularisation parameter of the field-data inversion i . In this section, we only discuss the lateral resolution. The depth resolution is briefly discussed in Section 2.5. Fig. 9 shows some of the most relevant combinations. A feature we extract from these figures is the capability to reproduce the simulated checker inside the area limited by the black polygon. We define the polygon by selecting tested the areas with lower root-mean square error (RMSE), where the vertices are seismic station locations. The RMSE is calculated by first estimating the residuals (the difference between the measured and predicted travel-times) and then by averaging the squared residuals and taking the square root. In equivalent matrix form, the RMSE is . The area inside the polygon is transected by a large number of ray paths with varying orientations for all the analysed frequencies. Even though in some areas outside the polygon (especially along the southwest-northeast direction), Fig. 9. Checkerboard test matrix for combinations of 1 km, 2 km and 3 km grid cells and anomalies sizes between 2 and 6 km. The black polygon identifies the area in which the checker is recovered best. the number of ray paths suffices to estimate slowness values, the lack of multidirectional sampling prevents accurate retrieval of the velocity structure. The bad performance in reproducing the checker outside the polygon occurs because of the broad-band seismic network configuration. Based on these results we drop grid cells outside the polygon for further tomography comparison between resolutions. For the frequency to depth inversion we keep the grid cells outside the polygon if there are enough frequency values to perform the inversion. In Fig. 9 results are shown for 0.18, 0.28 and 0.38 Hz as sampled examples from the used bandwidth. However, we invert all frequencies independently from 0.18 Hz to 0.44 Hz with a spacing of 0.02 Hz and use all these frequencies for further inversion to depth.
The regularisation parameter ( i ) can also be a measure of how well the regularisation is fitting the data given the noise. A small i can indicate that the noise contaminates the solution as it is overfitting d (meaning nearly no regularisation is applied), and a high i can increase solution smoothness and indicate a possible (not reasonable) estimation of m. Fig. 10 shows the variation of the regularisation parameter as function of resolution and frequency. As higher resolutions result in lower regularisation parameters consistently for all frequencies, we divide the regularisation parameter (dx) by the corresponding spatial resolution to provide a fair comparison ( i /dx with 1 ≤ dx ≤ 4 km).
There are no large discrepancies between the estimated regularisation parameters at the same frequencies for different resolutions (when normalised by the spatial resolution of each parameter). Nonetheless, between 0.2 Hz and 0.34 Hz, the preferred regularisation parameter is approximately the double of the value derived from frequencies between 0.16 and 0.18 Hz, and 0.36 and 0.44 Hz. The estimated regularisation parameter seems to be higher for the frequency interval with better quality SNR and a larger number of ray path coverage.
We invert to frequency-dependent Rayleigh wave velocity dispersion curves for different resolutions (1, 2, 3 and 4 km) which we depict in Fig. 10 b. Each line represents the tomographic result for each grid cell with estimated velocities for different grid cell sizes (1 to 4 km). The dispersion curves are smooth for almost all frequencies, with an exception for few grid cells for 1 and 2 km resolution. In theory, and if the ray-path coverage would be the same for each frequency, we would expect the regularization parameters to be the same for all frequencies, indicating that the suitability of the data for use in this inversion problem is equivalent between frequencies. The regularization parameter with the lower standard deviation is the one of 1 km resolution. From Fig. 10 we observe that the depth-dependent velocity estimation using the 3 km resolution is the higher resolution with monotonically decreasing dispersion curves, a requirement for phase velocities (Liu et al., 1976). The depth-dependent velocities for 1 km resolution also seem less noisy, with only few dispersion curves not monotonically decreasing.

Depth velocity estimation
We estimate the depth-dependent S-wave velocity tomography using the methodology of Wathelet (2008), an improved implementation of the neighbourhood algorithm (NA) described by Sambridge (1999). This estimation is another ill-posed problem, as it relies on the inversion of smooth dispersion curves into abrupt depth-dependent velocity changes, imposed by a given depth parametrisation. For each estimated dispersion curve (each grid cell) we run ∼30,000 inversions from which we extract the best model with minimum misfit. We run the models to calculate S-wave velocity at five fixed horizontal layers, ranging from 1000 to 6000 m depth to achieve 1 km resolution depth. We assume a Poisson ratio linearly varying with depth between 0.24 and 0.28 and a fixed density of 2600 kg/m 3 .

Results
We use the fundamental mode Rayleigh-wave arrival time produced by ANSI in a tomographic inversion scheme, to obtain frequency-dependent velocities. Fig. 11 depicts the results of the tomography inversion for the four grid cell tested resolutions (1 to 4 km). The positive and negative anomalies are estimated from an inverted average velocity V 0 per frequency derived from the tomographic inversion. To facilitate the comparison of resolution results in Fig. 11 we show two out of the 14 inverted phase-velocity maps (from 0.18 Hz to 0.44 Hz with a 0.02 Hz spacing). For more details on other frequencies see Section A.
The retrieved velocity anomalies have variations with maxima around 15% from an estimated average velocity V 0 . In Fig. 11 these variations are plotted between -10% and 10% to facilitate visualisation. All the variations above 5% are highlighted and compared between resolutions. We observe that the location of the main anomalies does not differ much between resolutions within the same frequencies (see red and blue circles in Fig. 11), even though higher resolutions (1 km) detect smaller anomalies which lower resolution grid cells fail to recognise (3 and 4 km) (see Section 2.4 and Section 4.1 for more details). The root mean square error (RMSE), measures of imperfections between the fit of the estimator and the data (d), which is lower in the tomographic results with 1 km resolution indicating a better fit to the data. The similarity of anomaly locations between different resolutions is especially interesting considering that each of the sub-figures is inverted independently (with a single design matrix G and regularisation parameter for each frequency value f i as described in section 2.3).
A comparison between both 1 km and 3 km resolution (the resolutions with smoother dispersion curves) allows us to do an additional independent check on the inversion to depth performance while trying to retrieve a higher horizontal spatial resolution. In Fig. 12 we show the depth-dependent S-wave velocity in a 3D field estimated through the procedure described in section 2.5. The velocity variations are displayed with respect to the average dispersion of the inverted results at each depth. We identify low and high velocity anomalies with matching locations between the 1 km and 3 km resolutions. The red and blue circles show the locations with good match for both resolutions and the red and blue arrows in Fig. 12 identify the locations where the match between the anomalies using 1 and 3 km resolution is not good. Most of the seismicity indicated in Fig. 11 happened to the east of the investigated area.

Resolution
The spatial resolution of the tomography is determined by the location, the number of seismic stations (grid spacing and seismic network aperture as function of azimuth) and the frequency content of the data. These parameters define the number of possible ray paths crossing each grid cell, and how well the ray paths are sampling each cell area. Additionally, the frequency dependent inter-station distance filtering (see Section 2.2 for details), adds a relation between the highest achievable resolution and frequency bandwidth. Higher resolutions will require shorter inter-station distances for higher frequencies while larger inter-station distances for lower frequencies. The effect on the resolution will be the loss of ray coverage at the edge of the seismic network for high frequencies and at the middle of the seismic network for low frequencies. From the checkerboard tests, it seems reasonable to use either 1, 2 or 3 km resolution as the checkers are well reproduced, at least inside the defined polygon. Note that the estimated RMSE shown in Fig. 9 is maximized as it also takes into account the area outside the polygon and this bias is the same for all the reproduced checkers. We observe that preferable resolutions depend on the size of the simulated checker anomalies. The larger the simulated anomalies, the better the checkerboard test recreates the checker. Based on these findings, we estimate the depth-dependent velocities using the 3 km resolution (which is the higher resolution without noisy dispersion curves), and for 1 km resolution, while dropping out the dispersion curves which are not monotonically decreasing. The choice of inverting for two resolutions allows an additional independent check on the consistency of the inversion from frequency to depth between resolutions.
The 1 km resolution samples a larger area around the Reykjanes geothermal field, not possible to detect with coarser resolutions, which is a region of interest given the location of most of the geothermal wells in the peninsula. However, while trying to achieve a higher resolution, the direct wave approximation is inherently violated and for resolutions below 3-4 km the quality of the velocity variation estimations may be reduced. This implies that the direct wave assumption only holds if the wavelength of the waves is larger than the scale of the medium heterogeneities. While the observed mismatch between both resolutions (arrows in Fig. 12) could be due to the anomaly size (resolution should be higher than twice the anomaly size), this would also only hold if the direct wave approximation would not play a role. Implications on the deterioration of the tomographic results for resolutions below 3 km would need further research, which is out of the scope of this paper. Therefore, we are careful interpreting the 1km resolution results and focus our interpretation using the 3km resolution tomographic image.

Insights on the geothermal fields
The retrieved S-wave velocity variations can occur as a result of wave propagation through different geological media, specific tectonic features or rock state (solid/melting/partial melting). Conditions such as crack alignment (isotropic or anisotropic fault swarms), composition (e.g. relation between minerals, shale content, fluid content), saturation (porosity, permeability, water content), pressure and temperature determine the speed of seismic waves (Biot, 1956;Gassmann, 1951).
The main contribution due to effective pressure can be observed in Fig. 12, the deeper structures show higher velocities (e.g. mean velocity v 0 at higher depths), as velocity usually increases with effective pressure. However, that is not always the case. As the effective pressure is defined by the difference between the confining pressure and the pore pressure, the pore pressure determines the effective pressure at the same depths. This is under the assumption of the same confining pressure for the same depths, which is not the case for the WRP tectonic setting (see description in Section 1), as expected in a tectonically and magmatically active sites. In the study area there are three zones of fault swarms (Clifton and Kattenhorn, 2006): 1. within the Reykjanes (R) geothermal field, 2. within the Eldvörp (E) and Svartsengi (S) geothermal fields 3. North of Reykjanes (R) and Eldvörp (E).
On top of the aforementioned contributions, in geothermal areas it is well known that alteration in the hydrothermal systems, salinity, shale content, and water saturation (O'Connell and Budiansky, 1974) interfere with the wave speed, as well as with changes in the rock state (solid/melting/partial melting). Variations in fluid content and temperature are likely to be the most determining parameters for velocity anomalies in geo-and hydro-thermally active areas (Nakajima et al., 2001b), and while partial melting areas might not exist, inclusions filled with H 2 O (Nakajima et al., 2001a) seem reasonable to exist in Reykjanes Peninsula. Therefore, the present study should be interpreted together with constraints from other measurements.
The derived velocity anomalies using 3km resolution do not cover the Reykjanes geothermal field preventing us from any possible interpretation. With the 1 km resolution, the observed local low-velocity anomalies at the Reykjanes geothermal field match the location of the intensively exploited part of the geothermal reservoir. This could potentially indicate that the observed lowvelocity anomalies might be related with a heat source, water inclusion or both (O'Connell and Budiansky, 1974;Mavko, 1980). The same location is described by Friðleifsson et al. (2018) as the area of up-flow targeted by IDDP-2 interpreted as hotter and more permeable. And a resistivity model based on 3D inversion of MT data (Karlsdóttir et al., 2018) indicate that the IDDP-02 well was drilled into a low resistivity anomaly, which coincides with the observed low S-wave velocities in Fig. 12 (red arrow in a to d at the Reykjanes geothermal filed location). However, we would need to Fig. 11. Phase-velocity maps after tomographic inversion using two frequencies (0.28 on the left and 0.38 Hz on the right), for the four tested spatial resolutions with edited colorbar between −10 and 10 percentage deviation from the average velocity (V0). The results are interpolated to the same grid cell size for figure comparison. Green triangles identify the seismic stations and the squares show the location of the geothermal fields identified in 1. The black polygon outlines the area where there is enough ray-path coverage to better recreate the simulated checkers of Fig. 9. Red ellipses identify similar low-velocity anomalies at different resolutions in the same frequency band. Blues circles identify similar high-velocity anomalies at different resolutions in the same frequency band. Fig. 12. 3D S-wave percentage velocity variations from the average velocity per depth (V0) between 2 and 6 km depth. Left column shows the results from 1 km and right column 3 km spatial resolution tomographic inversion for all the identified depths. We use the convention of red for low and blue for high velocity anomalies. Black dots indicate the locations of the earthquakes detected and provided by IMO (Icelandic Meterological Office) from 1993 to 2015. Blue lines locate the fault system within the defined polygon. Green triangles identify the seismic stations and the squares show the location of the geothermal fields identified in Fig. 1. Red circles and arrows correspond to areas where anomalies match between resolutions, while blue circles and arrows show where no match is observed. Fig. 13. Tomographic results of all frequencies for 3 km resolution. The polygon outlines the area where there is enough ray-path coverage 9 . In this figure the colorscale shows variations from −15% to 15% from the average velocity per frequency. verify the influence of the direct wave assumption for such high resolution to be able to confirm such interpretation.
We identify low-velocity anomalies which in volcanic areas are usually associated with partially molten magma pockets (or reservoirs depending on their size) or highly fractured conduit systems with possible molten rock (e.g. Benz et al., 1996;Sudo and Kong, 2001;Villasenor et al., 1998) as interpreted by Nakajima and Hasegawa (2003). These anomalies are mostly present in the high density fractured area north of Eldvörp field (between 2 and 5 km depth) and at the location of Svartsengi geothermal field (between 3 and 6 km depth) Fig. 12 from e) to g). The geology in this area is composed by post glacial lavas, with late Pleistocene hyaloclastites at the location of the Svartsengi geothermal field. At the same location of high-temperature thermal areas and high-positive magnetic anomalies (Jakobsson et al., 1978). Given the depth of the Svartsengi anomalies, these might be related with molten or partially molten pockets of magma, which may be the heat source of the surface geothermal manifestations (e.g. the Blue Lagoon, located within red square of Svartsengi geothermal field in Fig. 1). At the Eldvörp geothermal field we do not identify any low-velocity anomaly within the studied depths, but north of Eldvörp field we detect low speed anomalies between 3 and 5 km depth in the area of high fracture density. Franzson (1987) suggests that the dense fissure swarm around Eldvörp provides downflow channels of cold groundwater, while within the central part of the "swarm" an upflow of the hydrothermal system occurs along the same fractures. The derived 3D S-wave could eventually help to support this hypothesis with a detailed comparison with the existing borehole data from the area.
In general, there is an increase of low velocities from west to east indicating that the temperature may be higher towards the interior of the peninsula. Between 4 and 6 km depth a high-velocity anomaly seems to separate the Reykjanes geothermal field from the remaining fields of the peninsula. The high S-wave velocity structure that might indicate the presence of solidified rocks that have been cooled down possibly due to the proximity of the ocean.

Future applications
This study shows the potential of ambient noise tomography as a complementary reservoir characterization tool for field operations. From a seismic network design point of view, the performed resolution tests can also elucidate on the optimization of seismic station locations by using the methodology of Toledo et al. (2018) extended for ambient noise. In a similar manner, our results also show potential to complement the interpretation of deformation studies over the same area by (e.g.) 1. interpreting how the estimated horizontal displacements of Hreinsdóttir et al. (2001) or 3D surface motion of Gudmundsson et al. (2002) can be observed by spatial changes in S-wave velocity field as a result of shear strain; 2. constraining the solutions on the local sources of man-derived subsidence derived by Keiding et al. (2010) and Parks et al. (2018).

Conclusions
We have retrieved surface waves from ambient noise over approximately one year and five months using seismic interferometry and applied tomographic imaging of the subsurface covering the southwestern area of Reykjanes Peninsula. We test tomographic inversion with four spatial horizontal resolutions. We further invert to a depth dependent S-wave velocity model the tomographic result with the highest resolution and the resolution with more stable dispersion curves (1 km and 3 km, respectively). We detect high S-wave velocity structures that might indicate the presence of old intrusions that are now solidified and stiffer or cooling down effects due to ocean proximity. The low S-wave velocity anomalies can be interpreted as pockets of possible partial melt, for which the Vp/Vs ratio could give further insight. We detect S-wave velocity anomalies with variations between −15% and 15% with reference to an estimated average velocity between 2 and 6 km depth with 3 km of lateral resolutions and 1 km of vertical resolution. However, we refrain from interpreting the 1 km resolution tomography as real effects as the assumption inherent in the straight ray approximation might cause deterioration in quality of the results for higher resolutions. Although the used seismic network was not designed for ANT only, it accommodates different seismic imaging techniques and the results underline the capabilities of ANSI in Iceland.
In Fig. 13 we show the tomographic results for each frequency (steps of 0.02 Hz) used for the inversion from frequency to depth. As it can be observed, the anomalies identified for each frequency are smooth when compared with the adjacent frequencies. Considering that each frequency is inverted independently, the observed smoothness reinforces the reliability of the tomographic results.