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

Document Version

Final published version

Published in

Journal of Volcanology and Geothermal Research

Citation (APA)

Martins, J. E., Weemstra, K., Ruigrok, E., Verdel, A., Jousset, P., & Hersir, G. (2019). 3D S-wave velocity

imaging of Reykjanes Peninsula high-enthalpy geothermal fields with ambient-noise tomography. Journal of

Volcanology and Geothermal Research, 391, [106685]. https://doi.org/10.1016/j.jvolgeores.2019.106685

Important note

To cite this publication, please use the final published version (if applicable).

Please check the document version above.

Copyright

Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy

Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.

This work is downloaded from Delft University of Technology.

ContentslistsavailableatScienceDirect

Journal

of

Volcanology

and

Geothermal

Research

jou rn al h om ep a g e :w w w . e l s e v i e r . c o m / l o c a t e / j v o l g e o r e s

Invited

Research

Article

3D

S-wave

velocity

imaging

of

Reykjanes

Peninsula

high-enthalpy

geothermal

fields

with

ambient-noise

tomography

J.E.

Martins

_,

_C.

_Weemstra

_,

_E.

_Ruigrok

_,

_A.

_Verdel

_,

_P.

_Jousset

_,

_G.P.

_Hersir

a_{NetherlandsOrganizationforAppliedScientificResearch,TNO,Utrecht,TheNetherlands} b_{DelftUniversityofTechnology,Delft,TheNetherlands}

c_{RoyalNetherlandsMeteorologicalInstitute,DeBilt,TheNetherlands} d_{UtrechtUniversity,Utrecht,TheNetherlands}

e_{Helmholtz-ZentrumPotsdamGFZ,Potsdam,Germany} f_{IcelandGeoSurveyISOR,Reykjavik,Iceland}

a

r

t

i

c

l

e

i

n

f

o

Received31January2018 Accepted8October2019 Availableonline6November2019 Keywords:

Velocityanomalies Seismicinterferometry Empiricalgreenfunctions Surface-waves

Reservoircharacterization Modelresolution

a

b

s

t

r

a

c

t

Tomographicimagingbasedonambientseismicnoisemeasurementshasshowntobeapowerfultool, especiallyinareaslikeIceland,wherethemicroseismilluminationisexcellent.Inthispaper,weproduce a3DS-wavetomographicimageoverthewesternReykjanesPeninsulahigh-enthalpygeothermalfields andevaluatethereliabilityofthetomographicresultsfordifferentresolutionsthroughsimulatedand realdata.Weuse30broadbandstationsoperatingforapproximatelyone-and-a-halfyearandapply ambientnoiseseismicinterferometryforeachstation-pair.ThisresultsinempiricalGreen’sfunctionsin whichespeciallytheballisticsurfacewaves(BSW)arewellresolved.TheretrievedBSWexhibitahigh signal-to-noiseratiobetween0.1and0.5Hz,andthebeamforminganalysisindicatesanapparent surface-wavevelocityof3km/soverabroadazimuthalrange.Forthetomographicinversion,weinvertthe estimatedphasevelocitiesbetweenallstationpairstofrequency-dependentphasevelocitymapsinfour differentresolutions(1,2,3,and4km)usingaTikhonovregularisation.Withtheestimatedregularisation parameterperfrequencyperresolution,weinvertsimulateddataforcheckerboardsensitivitytestsper frequencyfordifferentcombinationsofvelocityanomalysizesandresolutions.

Finally,aftertheinversiontodepth,wedetectS-wavevelocityanomalieswithvariationsbetween−15% and15%withreferencetoanestimatedaveragevelocityusing1kmand3kmoflateralresolutionsand1 kmofverticalresolution.Thisstudyshowsthepotentialofambientnoisetomographyascomplementary seismologicaltoolforreservoircharacterization.

©2019TheAuthors.PublishedbyElsevierB.V.ThisisanopenaccessarticleundertheCCBYlicense (http://creativecommons.org/licenses/by/4.0/).

1. Introduction

TheWesternReykjanesPeninsula(WRP),locatedatthe

south-westerntipof Iceland (Fig.1)is a transitional zone betweena

tectonicspreadingcentre(theextensionoverIcelandofthe

Mid-AtlanticRidge)andatransformzone(PálmasonandSæmundsson,

1974; Sigmundsson et al., 2018). At this transitional zone, an

extensionalcomponentwithnortheast-southwesttrendingnormal

faulting,andatransformcomponentwithstrike-slipfaulting

ori-entednorth-southcanexplainobserveddisplacementcomponents

(Kleinetal.,1973;Einarsson,2008).Thistectonicsettingfacilitates

theintrusionofmagmaindikesalongtheNE-SWtrendingareaof

denselyspacedfissuresandfaults.Thepresenceofmagmatic

intru-∗ Correspondingauthor.

E-mailaddress:joana.estevesmartins@tno.nl(J.E.Martins).

sionsandfaultswarms,inturn,promotescooler”fresh”ground

waterinflowandthermalupflow(Franzson,1987).Theresultant

geothermalmanifestationshavemadetheWRPsubjectof

geother-malpowerproductionforover30years.

Currently, attempts are made to explore deep geothermal

sources, the supercritical part of hydrothermal systems, which

offersdifferentpotentialforharnessinggeothermalenergy.One

ofthemostnoticeableadvantagesofsupercriticalsystemsisthat

theincreaseofpressure(>221barsor22.06MPa)andtemperature

(>374 ◦C higher for seawater) can potentially induce an

esti-matedpoweroutput∼10timeshigherthantraditionalIcelandic

geothermalwells(Friðleifssonand Albertsson,2000;Albertsson

etal.,2003).InIceland,whilecustomarygeothermalsystems

usu-allyreachdownto3kmdepthwithsteamtemperaturesbetween

290◦Cand320◦C,predictionsbyFriðleifssonetal.(2014a)report

thatdeepwellsreachingdepthsbetween4kmand5kmmayhave

temperaturesbetween400◦Cand600◦C.Theseresultsweretested

https://doi.org/10.1016/j.jvolgeores.2019.106685

Fig.1. MapofWRPseismiccampaign,faults,fracturesandlocationinIceland.a)WRPcoastaloutline.Thegreentrianglesdenotetheonshorebroadbandstations.Thered dotsidentifywellsheadsinReykjanesgeothermalfields,andtheblacklinestheidentifiedfaultsandfracturesinthevicinityofthegeothermalproductionareas.Thesquares locatehightemperatureareas(mappedbytheIcelandicGeoSurveyISOR(Guðnasonetal.,2014))fromwesttoeast,Reykjanes(R),Eldvörp(E),Svartsengi(S),ofwhichthered squaresindicatethelocationofthetwoexistingpowerplants.b)IcelandiccoastalboundarywithReykjanespeninsulawithintheredrectangle.Reddashedlinelocatesthe plateboundaryandblacklineslocatethetransformfaultsystem.Theneovolcaniczonesarefilledinlightshadedgreyalong(EinarssonandSaemundsson,1987;Erlendsson andEinarsson,1996;Einarssonetal.,2002;Einarsson,2008;Sigmundssonetal.,2018).

andpartiallyconfirmedwithIDDP2inReykjanes(drilledbetween

2016and2017)withmeasuredtemperature427◦Candfluid

pres-sureof340barsat∼4.5kmdepth(Friðleifssonetal.,2017a).In

2009,withintheIcelandicDeepDrillingProject(IDDP),thefirst

IDDPwell(IDDP-1)drilledinIcelandintheKraflaarea

encoun-teredrhyoliticmagmaatadepthof2104m(Hólmgeirssonetal.,

2010).Magnetotelluric(MT)electromagneticsurveysperformedat

Kraflawronglyestimatedamagmasourceat4.5kmdeep(Elders

etal., 2011;Gasperikovaet al.,2015), and thedrilling

encoun-teredshallowermagma.EventhoughMTsurveysarethepreferred

geophysicalmeasurement procedureto locate thewater

reser-voirsrequiredfortheclosedloopofgeothermalproduction(e.g.

Newmanetal.,2008;Árnasonetal.,2010;Hersiretal.,2018)and

havebeenutilizedfrequentlyinIceland(e.g.Hersiretal.,1984;

EysteinssonandHermance,1985;Oskooietal.,2005),themagma

pocketthatwasdrilledintoattheKraflaIDDP-1welllocationwas

possiblybelowthelevelofresolutionoftheareasurveyedusingMT

(Árnasonetal.,2007).TheIDDP-2welldrilledinDecember2016at

theWRPsalinegeothermalsysteminsouthwestIceland

success-fullyreachedapproximately4.5km(vertical)depth(Friðleifsson

etal.,2017b;Eldersetal.,2014;Friðleifssonetal.,2018).Forthe

locationofthewell,theoperatorscouldconsiderresultsfrom

resis-tivitymeasurements (Karlsdóttir etal., 2018;Friðleifssonet al.,

2014b;Darnetetal.,2018),aswellastraveltimeseismic

tomogra-phy(Joussetetal.,2017).Theseismictomographywasestimated

fromaseismiccampaigndevelopedundertheEuropean-funded

programIntegratedMethodsforAdvancedGeothermalExploration

(IMAGE).The Joussetet al. (2017) tomographic studyconfirms

thepreviousresultsofBjarnasonetal.(1993)andofTryggvason

etal. (2002)over the sameareawithenhanced details around

welllocations. The authors interpret thelow-ratio anomaly of

compressional-overshear-velocityasbeingduetotheabsenceof

asizeablemagmaticbodyatthetipofReykjanesPeninsula,which

wasconfirmedbyFriðleifssonetal.(2018).

TheunexpecteddrillingintoamagmasourceinKraflain2009

highlighted the need to explore high-resolution imaging

tech-niquesasacomplementtocurrentmeasurementmethodsandto

improvetheexisting ones.Seismictomographictechniquesand

recent advances using ambientnoise-based methodologies can

playaroleinassessingthenecessarydepthandresolution

infor-mationandinconstrainingothergeophysicalestimations.Inthis

regard,ambientnoise seismicinterferometry(ANSI)techniques

canofferadditionaladvantagesbyavoidingthecostofactive

seis-mic methods (therefore, making ANSItechniques economically

moreattractive),andcircumventinglimitationssuchasthe

lim-itednumberofearthquakesandirregularearthquakedistribution.

Thestraightforwarddataacquisitionandthetheoreticalconcepts,

extendedfrom1Dmedia(Claerbout,1968)toarbitrarily

hetero-geneous 3D media by Wapenaar (2004), make ANSI attractive

for tomography applications.The ANSIconcept relies ona

vir-tual sourcethat is generated atthe locationof one of thetwo

receiversbycross-correlationandsummationof(noise)

record-ingsfromsurroundingambientnoisesources. Thetomographic

resultsare subsequentlyderivedfromRayleigh(orLove)waves

retrievedbetweenthevirtualsourcesandthereceivers.The

num-berofapplicationsofambientnoisetomography(ANT)studieshas

increasedinrecentyears,especiallyinIceland(Obermannetal.,

2016;Benediktsdóttiretal.,2017;Jeddietal.,2017;Greenetal.,

2017;Martinset al.,2019).Along withthedirect advantageof

characterisingthesubsurfacewithinthedepthrangeofIcelandic

geothermaloperations,ANTcancontributetoconstrainingother

geophysicalmeasurementsorinterpretationsthathavepreviously

beenacquiredoverthesameareaandviceversa.ANTcanbeusedto

improveasubsurfaceimagewithcomplementaryseismicstudies,

(Verdeletal.,2016;Blancketal.,2019;Joussetetal.,2017).

Inthisstudy,wederivea3DS-wavevelocitytomographicimage

oftheWRP’ssubsurface byapplyingANT totheseismicsurvey

deployedundertheIMAGEprojectframework.Ontopofassessing

thereliabilityoftheretrieveddispersioncurves,wedevote

partic-ularattentiontothemodelresolutiongiventhedeployednetwork

configurationandtoobtainresultsthatallowtoconstrainother

geophysicalmeasurements.

2. Dataandmethodology

Withinthe IMAGEprojectframework, theGerman Research

Center for Geosciences (GFZ Potsdam) and Iceland Geosurvey

Table1

Stationcoordinatesoftheseismicnetworkusedinthisstudy.

Stationcode Latitude[◦N] Longitude[◦N] Sensor

BER 63,818466 −22,562944 TrilliumComp

EIN 63,856934 −22,619266 TrilliumComp

GEV 63,828094 −22,466464 TrilliumComp

HAH 63,928601 −22,638602 TrilliumComp

HAS 63,882949 −22,715219 TrilliumComp

HOS 63,947675 −22,089654 TrilliumComp

KEF 64,016213 −22,627548 TrilliumComp

KUG 64,009625 −22,139352 TrilliumComp

LFE 63,883868 −22,535548 TrilliumComp

ONG 63,818635 −22,727764 TrilliumComp

PAT 63,953996 −22,532067 TrilliumComp

PRE 63,886239 −22,633336 TrilliumComp

RAH 63,852855 −22,567946 TrilliumComp

RAR 63,825801 −22,678328 TrilliumComp

RET 63,806745 −22,700812 TrilliumComp

SDV 63,821768 −22,633443 TrilliumComp

SKG 63,863371 −22,32921 TrilliumComp

SKH 63,904584 −22,414933 TrilliumComp

STA 63,854302 −22,697544 TrilliumComp

SUH 63,852128 −22,502155 TrilliumComp

ARN 63,862844 −22,04607 Marksensor

HOP 63,845073 −22,39242 Marksensor

KHR 63,832367 −22,596432 Marksensor

KRV 63,812744 −22,660527 Marksensor

MER 63,883236 −22,228253 Marksensor

NEW 63,932983 −22,385217 Marksensor

SKF 63,811912 −22,687761 Marksensor

STF 63,913445 −22,547067 Marksensor

STK 63,899571 −22,697866 Marksensor

VSR 63,873686 −22,588137 Marksensor

WRP.Ofthese,24 areoceanbottomstations (OBS),and30 are seismometersplacedonshore(Blancketal.,2019,Joussetetal., 2017).

TheonshoreinstrumentswereoperatingfromMarch2014,and

theOBSwereplacedinAugust2014.Alltheequipmentwas

col-lectedinAugust2015.TheOBSdeployedduringtheIMAGEproject

werenotusedinthisstudyduetoaphaseshiftintheinstruments

(Weemstraetal.,2016).Weusedthe30onshore(Table13

compo-nentsseismometers(20broadbandTrilliumcompactsensorsand

10short-periodMarkSensors)withacornerfrequencyaslowas

0.005Hzandasamplingrateof200Hzforadurationofalmost

oneyearand fivemonths.Weonlyusethevertical-component

displacementsandtoreducecomputationtime,wedown-sampled

therecordsto25samplespersecond(Nyquistof12.5Hz).

Theappliedmethodologyfollowstheprocessingapproachof

Martinsetal.(2019)integrallyandtheprocessingchainisdepicted

inFig.2withadivisionbetweendata-processingand inversion

schemes (identified in Fig. 2 by the orange and green dashed

squares,respectively).We dividetheprocessingchain intotwo

steps: data processing of the retrieved surface-waves and the

inversionprocedure. Thepre-processing includesdeconvolution

oftheinstrumentresponses,spectralwhitening,temporal

aver-aging(Bensenetal.,2007)andfrequencyfiltering.Thenextstep

isEmpiricalGreen’sFunctions(EGF)retrievalwithseismic

inter-ferometry.Thetomographicinversion schememakesuseofthe

estimatedBallistic SurfaceWaves (BSW)arrival-times fromthe

EGFtoestimatefrequencydependentspatialvelocityanomalies.

Finally, we estimate the depth-dependent 3D S-wave velocity

anomalies.

Intherestofthissection,wedescribethedatapre-processing,

EGFretrieval,BSWarrivaltimepicking,tomographyandinversion

toS-wave-velocity.Betweenthesesteps,weperformqualitychecks

toensurethat:1)Thesignal-to-noiseratioishighenoughtoallow

acceptablesurface-wavearrivaltimeestimation.2)The

illumina-tionissufficientlyuniform(takingintoaccountthatalackofcausal

illuminationcanbecompensatedbytheacausalpartofthe

cross-correlatedsignal).3)Theestimatedtimepicksareconsistent.4)

Smooth velocityvariations betweenthe tomographicfrequency

dependentresultswhileinverted independently,and5)Weare

choosingthemostappropriateresolution.

2.1. BSWfromcross-correlations

Theinstrumentresponseisremovedbycomplexdeconvolution,

afterwhichweapplyaspectraldomainnormalisation(i.e.

whiten-ing)(Bensenetal.,2007).WeextractcoherentEGFsbetween0.1

and0.5Hzafteraspectrogramexaminationoftheambientnoise.

Thisbandwidthisdominatedbythemicroseisms.Intheraw

spec-trogram(Fig.3a)weidentifyahighertime-dependent(seasonal)

powerspectraldensity(PSD)bandaround0.2Hz.Inthesame

fig-ure,thesharplinescoveringtheentirefrequencyspectrumwith

largePSDmarktheoccurrenceofearthquakes.Theseasonaleffect

withincreasedPSDfrommid-AugusttoApril(topFig.3)occurs

duetothehighernumberofoceanstormsandlargerwavesinthe

autumnandwinter(Ardhuinetal.,2011).Fromthespectrograms

weseethatthereissufficientenergyupto0.8Hz.

Wecomputethecross-correlationsperhourandstackthe

com-puted cross-correlationsusing approximatelyoneyear andfive

monthsofrecordedseismicdata.Fig.4showstheresultingEGF’s

obtainedbetweenthe435uniquestationpairs,usingdatabetween

0.1and0.5Hz.TheextractedEGF’sarehighlycoherentandshow

a non-symmetrical(inamplitude) ’V’shape oftheBSWarrivals

indicating,asexpected,anon-isotropicazimuthaldistributionof

ambientnoisesources(Fromentetal.,2010).

Strictlyspeaking, theretrievedBSWsonly coincidewiththe

surfacewavepartoftheGreen’sfunctionanditstime-reversed

ver-sionundertheconditionthat(i)thereceiverpairsareilluminated

uniformlyfromallangles(WapenaarandFokkema,2006),(ii)a

singlesurfacewavemodedominatestherecordedambient

vibra-tions(HallidayandCurtis,2008),and(iii)themediumislossless.

Inpractice,andthereforealsoinourcase,theseconditionsarenot

fulfilled,leadingtodeviationsoftheextractedsurfacewave

veloci-tiesfromthetruesurfacewavevelocities,e.g.(Tsai,2009;Froment

etal.,2010).Iftheilluminationpatternissufficientlyuniformover

anangle-rangeofatleast180degrees,itsufficestouseonlythe

causaloracausalBSW.Ithasbeenshownthatthisstillallowsusto

obtainmeaningfultomographic(surfacewave)images(e.g.Shapiro

etal.,2005;DeRidderetal.,2015).

Fig.5a),b)andc)showthreeofthefilteredEGFsaround0.18

Hz,0.28Hzand0.38Hz,respectively,andcross-correlationsfiltered

withasmallfrequencywindowaroundspecifiedfrequencyvalues

(+/-0.01Hz).Thecoherent’V’shapedwavepatternoftheBSW

indicatesthatitiswellpossibletopickarrivaltimes.FromFig.5,we

canrecognisethatlowerfrequenciestravelfaster.Moreover,itcan

beseenthatatashortinterstationdistancesthereisinterference

betweentheBSWatcausalandacausaltimes.

Weestimatetheazimuthaldirectionsoftheillumination,using

a beamforminganalysisapplied tothecross-correlation results

(Ruigroketal.,2017).Fig.5showsthebeamformingresults(d,eand

f)forthreefrequencybandscoveringthehighestSNRoftheBSW

bandwidth[0.16,0.22],[0.22,0.32]and[0.32,0.44],respectively.

PersistentBSWarrivewithinmostdirectionsofthreeazimuthal

quadrants,between90◦ and360◦.Besidesthebackazimuth,the

beamformingalsoyieldsthehorizontalrayparameter(inverseof

velocityforsurfacewaves)oftheincomingwaves.The

beamform-ingresultsindicatethatnostationpairshouldbedroppedbecause

of the lack of illumination. The quadrant lacking surface wave

arrivals(between0◦and90◦)canbecompensatedbytheopposite

quadrant(between180◦and270◦).Theseresultsareinagreement

withsimilaranalysesinIcelandforthethreeanalysedfrequency

Fig.2. Illustrativeprocessing-chainflowchart.Eachstepisidentifiedwiththecorrespondingreferenceorfigureofthisstudy.

Fig.3.Powerspectraldensity(PSD)oftherecordings(April2014-August2015)bythestationEINlocatedinthecentreoftheseismicnetwork.Powerspectraareaveraged overfourhoursand1.38×10−3_{Hz.a)showsthe(non-normalized)PSDandb)thepowerspectraindividuallynormalized(i.e.,withrespecttothemaximumpowerina}

2.2. Surface-wavephasevelocitypicking

ThehighcoherenceoftheBSWbetweenadjacentinter-station

distances(Fig.4)allowsustoestimatewellthetimingatlocal

max-imaofthecorrelationwaveform(consideredtobethecorrectphase

cycle).Toobtainindividual(i.e.,perstationpair)phase-velocity

estimates,we firstextractan averagephase velocity dispersion

curveforthefundamentalmodeRayleigh wave(c(f),wherefis

frequency)inthefrequencydomainusingtheMASW

(multichan-nelanalysisofsurfacewaves)algorithm(Parketal.,1998,1999).

Theaveragedispersioncurveasafunctionoffrequency(Fig.6)

servesasfurtherguidancetoavoidphasecyclejumpsinthephase

pickingofindividualfrequencies(fi)withashortintervalaroundfi

(wherefi∈[0.1,0.5]Hz).Realisticvelocitiesareonlyestimatedfor

frequencieshigherthan0.14Hz.

Fig.6showsbothMASWpickingandtheaveragephase

veloc-itydispersioncurve(leftandrightpanel,respectively).Ontopof

theaveragephasevelocity,weplottheresultingphasevelocityif

differentcyclesweretobepicked(green,redandblackasterisks

Fig.6rightside).Selectingacorrectcycleisstraightforwardwith

Fig.4. EGFwithambientnoiseseismicinterferometryappliedtovertical-componentdata.Eachseismictracecorrespondstoastation-paircombinationorderedby inter-stationdistance.Forminguniquestationpairswith30broad-bandstationsresultsin435seismograms.Approximately1.5yearsofdataFig.3wasusedbetween0.1and0.5 Hz.

Fig.5.EGFfor3frequencybands(top)andcorrespondingillumination(bottom).a),b)andc)representtheEGFafterwhiteningandfrequencyfilteringaround(+/-0.01Hz) 0.18,0.28and0.38Hz,respectively.d),e)andf)representthebeamformingresultsforthreefrequencybandwidthscoveringthehighestSNRoftheBSW,bandwidth[0.16, 0.22],[0.22,0.32]and[0.32,0.44],respectively.

Fig.6. Leftpanel:pickinginfrequency-wavenumberdomainusingMASWalgorithm.Rightpanel:bluefilledlinerepresentsthetimedomainaveragedispersioncurveas resultofMASWpicking.Thegreen,redandblackasterisksshowthetravel-timepicks(reworkedtovelocities)ofthreerandomstationpairsfordifferentphaseunwrapping integers.Theblacklinesconnectingtheasterisksindicatetheselectedphasevelocitycurvec(f).

Fig.7. a)Cross-correlatedresultsandcorrespondingtime-picksforfrequencies(lefttoright)0.18,0.28and0.38Hz.b)Time-picks(reworkedtovelocities)within2sigma velocityvariationfromthemeanforfrequencies(lefttoright)0.18,0.28and0.38Hz.

theMASWvelocityreference.Foronecyclevelocitiesareobtained

thatareclosetotheMASWvelocities.Pickingatthemaximaof

othercyclesyieldsunrealisticvelocities(Fig.6right).Weimposea

thresholdtowithdrawstationpairswithdistanceswheretoomuch

destructiveinterference(orlowSNR)occurs,withconsequentloss

ofwaveformcoherence.Highercoherentcross-correlationsallow

anaccuratetime-pickingarrival. Considering i,0=

v_i,0

fi , weuse

thesamefrequency-dependentthresholdsasMartinsetal.(2019)

Ri,min= 23i,0andRi,max=2.8i,0todefinetheminimumandthe

maximumadmissibleinter-stationdistances,respectively.The“0”

subscriptreferstothereference(average)phasevelocityobtained

fromMASW.Foreachstationcombination,thereis(potentially)a

causalandanacausalBSWretrieval.WeusetheBSWthathasthe

highestamplitude.Fig.7showsthepicksinthetimedomain.Asan

additionalevaluation,thefigurealsoshowstheselectedoutliers,

thevaluesoutsidethe2deviationfromthemean.

Withtheestimatedtime-picksforeveryfibetween0.1Hzand

0.5Hzin stepsof0.02 Hz,wecalculatetheazimuthalvariation

ofthephasevelocity.Azimuthalvariationsperfrequencycanbe

seenasaqualitychecktodetectnon-smoothnessoftime-picksand

incasesofsuccess,togetherwithsufficientilluminationcanalso

indicatevelocityanisotropy.Weobservehigherazimuthal

veloci-tiesbetween40◦and80◦,approximatelyinanortheast-southwest

direction,andlowervelocitiesbetween100◦and200◦(Fig.8).The

highervelocitiesmaybecorrelatedwiththeorientationofthemain

faultsandridgethatcuttheWRP(Fig.1),butfurtherinvestigation

wouldberequiredtoverifythisobservation.Azimuthalvariations

areconsistentforallfrequenciesand havea smoothdifferential

Fig.8. Left:Azimuthalphasevelocityperfrequency.Right:EstimateofthedepthofmaximumsensitivityasfunctionoffrequencybybothXiaetal.(1999)andHaneyand Tsai(2015)theoreticalrelationsandthemaximumpenetrationdepthalsoasfunctionoffrequency.

picks(thepicksaremadeindependentlyforeachfi).Asa

refer-encefor frequency-depthcorrespondence,weadd toFig.8,the

depthofmaximumsensitivityoftheRayleighwaves,aswellas

themaximumpenetrationdepth,bothasafunctionoffrequency.

2.3. Tomography

Seismictomographyisaninversionproblemthataimstofinda

slownessfieldfromtravel-timeobservationsinthreedimensions.

Thelinearforwardproblemcanbewrittenas:

d=G·m+e, (1)

wheredthedatavector(orobservations)containingthepicked

frequency-dependenttraveltimes; mthemodel vector

describ-ingthe(unknown)frequency-dependentslownessvaluesofeach

gridcell;Gistheoperatormatrixcontainingtheraypathlengths

ofeachraypath(firstdimension)througheachgridcell(second

dimension)and;eisatermexpressingthenoise.Weusethesame

methodologyasdescribed in Martinset al.(2019),a first-order

Tikhonovregularisation(Tikhonov,1963)toregularisethis

inver-sionproblem,whichminimisesthedoubleobjectivefunction:

min

mRn{||d−Gm||

2_+||m||2_{} (2)}

Theterm||m||2 _{denotesthenormofthemodel,multiplied}

bytheso-calledregularisationparameteri.Theaddedobjective

ofminimisingthenormoftheestimatedparametersavoids

over-fittingofthedatabyenforcingsmoothing:

ˆ

m=(GTG+I)−1GTd, (3)

whereIistheidentitymatrixandwherethesuperscriptTdenotes

thetransposeofamatrix.Weusethecross-validation

methodol-ogy(Golubetal.,1979)tochoosearegularisationparameter()per

frequencyvalue(fi).Foreachfrequency,westartbyremovingone

observationttjwherettrepresentsthetravel-timebetweentwo

stations.ThenweuseEq.(3)tofitasolutionwithouttheremoved

observation( ˆmj).Wepredicttheneglectedtraveltimeusingthe

resolvedmjandcompareitwiththedroppedobservationitself:

rj=dj−Gj. ˆmj,whereGiislackingtherowassociatedwithstation

couplej.Weperformthisforallobservationsandsumthesquaresof

theresidualstochecktherobustnessofi.Werepeatthisfor

mul-tipleandweselecttheithatminimisesP= 1_n



j=1(Gjmˆj−dj) 2

,

wherenisthenumberofobservations,i.e.,thenumberofstation

couplesforwhichthephasevelocitywasestimated.

The tomographic linear inversion is repeated for different

frequencies to produce phase-velocity maps.Frequency can be

approximatedtodepthusingthetheoreticalrelationsinXiaetal.

(1999)andHaneyandTsai(2015)(Fig.8right).However,foramore

accuratefrequency-depthconversionitisadvisedtoperforma

sec-ondinversionusingtheRayleighwavesensitivitykernels(Zhou

etal.,2004).Thismethodisdescribed inSection2.5.Astheray

pathsperfrequencyareinvertedindependently,wealsoestimate

adifferentregularisationparameterperfrequencyvalue,iand

fi∈[0.18,0.44]with0.02Hzsteps.

2.4. Checkerboardresolution

Anadequatemodelresolutioncanhelptoidentifysubsurface

anomalies’geometries,whichisrelevantforsubsurface

character-isationandgeothermalpurposes.We usecheckerboardforward

modelling to test the ability of the seismic network

geomet-ric coverage to reproducea simulated checker for each of the

testedresolutionsversusanomalysize.Wetestcombinationsof

spatial resolutions (ranging from1 to4 km), and size of

sim-ulatedperturbations (withperturbationsizesrangingfrom2 to

6 km).Asthenumber and spatialdistributionof theray paths

changes withfrequency, wereproducea checkerboardfor each

frequencyfi.Similarly totheprocedurefor thefield-data

inver-sion,thecheckerboardinversionisdoneindependentlyforeach

fiusingtheestimatedTikhonov regularisationparameterof the

field-data inversion i.In this section,weonly discussthe

lat-eralresolution.ThedepthresolutionisbrieflydiscussedinSection

2.5.

Fig.9showssomeofthemostrelevantcombinations.A

fea-tureweextractfromthesefiguresisthecapabilitytoreproduce

thesimulatedcheckerinsidethearealimitedbytheblack

poly-gon. We define thepolygon by selecting tested theareas with

lowerroot-meansquareerror(RMSE),wheretheverticesare

seis-micstationlocations.TheRMSEiscalculatedbyfirstestimating

theresiduals(thedifferencebetweenthemeasuredandpredicted

travel-times) and then by averaging the squared residuals and

taking thesquare root.In equivalentmatrix form,the RMSEis

estimatedbyRMSE= d−G

ˆ

m



d−Gmˆ

.Theareainsidethepolygonis

tran-sectedbyalargenumberofraypathswithvaryingorientations

foralltheanalysedfrequencies.Eventhoughinsomeareasoutside

Fig.9.Checkerboardtestmatrixforcombinationsof1km,2kmand3kmgridcellsandanomaliessizesbetween2and6km.Theblackpolygonidentifiestheareainwhich thecheckerisrecoveredbest.

thenumberofraypathssufficestoestimateslownessvalues,the

lackofmultidirectionalsamplingpreventsaccurateretrievalofthe

velocitystructure.Thebadperformanceinreproducingthechecker

outside the polygon occurs because of the broad-band seismic

networkconfiguration.Basedontheseresultswedropgridcells

outsidethepolygonforfurthertomographycomparisonbetween

resolutions.Forthefrequencytodepthinversionwekeepthegrid

cellsoutsidethepolygonifthereareenoughfrequencyvaluesto

performtheinversion.InFig.9resultsareshownfor0.18,0.28and

0.38Hzassampledexamplesfromtheusedbandwidth.However,

weinvertallfrequenciesindependentlyfrom0.18Hzto0.44Hz

withaspacingof0.02Hzanduseallthesefrequenciesforfurther

inversiontodepth.

Theregularisationparameter(i)canalsobeameasureofhow

welltheregularisationisfittingthedatagiventhenoise.Asmalli

canindicatethatthenoisecontaminatesthesolutionasitis

over-fittingd(meaningnearlynoregularisationisapplied),andahigh

icanincreasesolutionsmoothnessandindicateapossible(not

reasonable)estimationofm.Fig.10showsthevariationofthe

reg-ularisationparameterasfunctionofresolutionandfrequency.As

higherresolutionsresultinlowerregularisationparameters

con-sistentlyforallfrequencies,wedividetheregularisationparameter

(dx)bythecorrespondingspatialresolutiontoprovideafair

com-parison(i/dxwith1≤dx≤4km).

Therearenolargediscrepanciesbetweentheestimated

regular-isationparametersatthesamefrequenciesfordifferentresolutions

(whennormalisedbythespatialresolution ofeach parameter).

Nonetheless,between0.2Hzand0.34Hz,thepreferred

regular-isationparameterisapproximatelythedoubleofthevaluederived

fromfrequenciesbetween0.16and0.18Hz,and0.36and0.44Hz.

Theestimatedregularisationparameterseemstobehigherforthe

frequencyintervalwithbetterqualitySNRandalargernumberof

raypathcoverage.

Weinverttofrequency-dependentRayleighwavevelocity

dis-persioncurvesfordifferentresolutions(1,2,3and4km)which

wedepictinFig.10b.Eachlinerepresentsthetomographicresult

foreachgridcellwithestimatedvelocitiesfordifferentgridcell

sizes(1to4km). Thedispersioncurves aresmooth foralmost

allfrequencies,withanexception forfewgridcellsfor 1and2

kmresolution.Intheory,andiftheray-pathcoveragewouldbe

thesameforeachfrequency,wewouldexpecttheregularization

parameterstobethesameforallfrequencies,indicatingthatthe

suitabilityofthedataforuseinthisinversionproblemis

equiv-alentbetweenfrequencies.Theregularizationparameterwiththe

lowerstandarddeviationistheoneof1kmresolution.FromFig.10

weobservethatthedepth-dependentvelocity estimationusing

the3kmresolutionisthehigherresolutionwithmonotonically

decreasingdispersioncurves,arequirementforphasevelocities

(Liuetal.,1976).Thedepth-dependentvelocitiesfor1km

reso-lutionalsoseemlessnoisy,withonlyfewdispersioncurvesnot

monotonicallydecreasing.

2.5. Depthvelocityestimation

Weestimatethedepth-dependentS-wavevelocity

tomogra-phy using the methodology of Wathelet (2008), an improved

implementationoftheneighbourhoodalgorithm(NA)described

bySambridge(1999).Thisestimationisanotherill-posed

prob-lem,asitreliesontheinversionofsmoothdispersioncurvesinto

abrupt depth-dependent velocity changes, imposed by a given

depthparametrisation.Foreachestimateddispersioncurve(each

gridcell)werun∼30,000inversionsfromwhichweextractthe

bestmodelwithminimummisfit.Werunthemodelstocalculate

S-wavevelocityatfivefixedhorizontallayers,rangingfrom1000

to6000mdepthtoachieve1kmresolutiondepth.Weassumea

Poissonratiolinearlyvaryingwithdepthbetween0.24and0.28

andafixeddensityof2600kg/m3_.

3. Results

We use the fundamental mode Rayleigh-wave arrival time

producedbyANSIinatomographicinversionscheme,toobtain

frequency-dependentvelocities.Fig.11depictstheresultsofthe

tomographyinversionforthefourgridcelltestedresolutions(1to

4km).Thepositiveandnegativeanomaliesareestimatedfroman

invertedaveragevelocityV0perfrequencyderivedfromthe

tomo-graphicinversion.Tofacilitatethecomparisonofresolutionresults

inFig.11weshowtwooutofthe14invertedphase-velocitymaps

(from0.18Hzto0.44Hzwitha0.02Hzspacing).Formoredetails

onotherfrequenciesseeSectionA.

Theretrievedvelocityanomalieshavevariationswithmaxima

around15%fromanestimatedaveragevelocityV0.InFig.11these

variationsareplottedbetween-10%and10%tofacilitate

visuali-sation.Allthevariationsabove5%arehighlightedandcompared

betweenresolutions. We observethatthe locationof themain

anomaliesdoesnotdiffermuch betweenresolutionswithinthe

samefrequencies(seeredandbluecirclesinFig.11),eventhough

higherresolutions(1km)detectsmalleranomalieswhichlower

resolution gridcells failtorecognise(3and4 km)(seeSection

2.4andSection4.1formoredetails).Therootmeansquareerror

(RMSE),measuresofimperfectionsbetweenthefitofthe

estima-torandthedata(d),whichislowerinthetomographicresultswith

1kmresolutionindicatingabetterfittothedata.Thesimilarity

ofanomalylocationsbetweendifferentresolutionsisespecially

interesting considering that each of the sub-figures is inverted

independently(withasingledesignmatrixGandregularisation

parameterforeachfrequencyvaluefiasdescribedinsection2.3).

A comparisonbetweenboth 1 kmand 3 kmresolution (the

resolutions with smoother dispersion curves) allows us to do

anadditionalindependentcheckontheinversiontodepth

per-formance while trying to retrieve a higher horizontal spatial

resolution.InFig.12weshowthedepth-dependentS-wavevelocity

ina3Dfieldestimatedthroughtheproceduredescribedinsection

2.5.Thevelocityvariationsaredisplayedwithrespecttothe

aver-agedispersionoftheinvertedresultsateachdepth.Weidentify

lowandhighvelocityanomalieswithmatchinglocationsbetween

the1kmand3kmresolutions.Theredandbluecirclesshowthe

locationswithgoodmatchforbothresolutionsandtheredandblue

arrowsinFig.12identifythelocationswherethematchbetween

theanomaliesusing1and3kmresolutionisnotgood.Mostofthe

seismicityindicatedinFig.11happenedtotheeastofthe

investi-gatedarea.

4. Discussion

4.1. Resolution

Thespatialresolutionofthetomographyisdeterminedbythe

location,thenumber ofseismicstations(gridspacingand

seis-micnetworkapertureasfunctionofazimuth)andthefrequency

contentofthedata.Theseparametersdefinethenumberof

pos-sibleraypathscrossingeachgridcell,andhowwelltheraypaths

aresamplingeachcellarea.Additionally,thefrequencydependent

inter-stationdistancefiltering(seeSection2.2fordetails),addsa

relationbetweenthehighestachievableresolutionandfrequency

bandwidth. Higherresolutionswillrequireshorter inter-station

distancesforhigherfrequencieswhilelargerinter-stationdistances

forlowerfrequencies.Theeffectontheresolutionwillbethelossof

raycoverageattheedgeoftheseismicnetworkforhighfrequencies

Fig.10. Left:Regularisationparameternormalizedbythecorrespondingresolution,asfunctionoffrequency.Right:dispersioncurvesoftheinvertedsurface-wavetravel timesforthecorresponding1,2,3and4kmresolutionsidentifiedintheleftfigure.

Fromthecheckerboardtests,itseemsreasonabletouseeither

1,2or3kmresolutionasthecheckersarewellreproduced,atleast

insidethedefinedpolygon.NotethattheestimatedRMSEshownin

Fig.9ismaximizedasitalsotakesintoaccounttheareaoutsidethe

polygonandthisbiasisthesameforallthereproducedcheckers.

Weobservethatpreferableresolutionsdependonthesizeofthe

simulatedcheckeranomalies.Thelargerthesimulatedanomalies,

thebetterthecheckerboardtestrecreatesthechecker.Basedon

thesefindings,weestimatethedepth-dependentvelocitiesusing

the3kmresolution(whichisthehigherresolutionwithoutnoisy

dispersioncurves),andfor 1kmresolution,whiledropping out

thedispersioncurveswhicharenotmonotonicallydecreasing.The

choiceofinvertingfortworesolutionsallowsanadditional

inde-pendentcheckontheconsistencyoftheinversionfromfrequency

todepthbetweenresolutions.

The1kmresolutionsamplesalargerareaaroundtheReykjanes

geothermalfield,notpossibletodetectwithcoarserresolutions,

which is a regionof interest given thelocation of mostof the

geothermal wells in the peninsula. However, while trying to

achieveahigherresolution,thedirectwaveapproximationis

inher-entlyviolatedandforresolutionsbelow3–4kmthequalityofthe

velocityvariationestimationsmaybereduced.Thisimpliesthatthe

directwaveassumptiononlyholdsifthewavelengthofthewaves

islargerthanthescaleofthemediumheterogeneities.Whilethe

observedmismatchbetweenbothresolutions(arrowsinFig.12)

couldbeduetotheanomalysize(resolutionshouldbehigherthan

twicetheanomalysize),thiswouldalsoonlyholdifthedirectwave

approximationwouldnotplay a role.Implicationsonthe

dete-riorationofthetomographicresultsforresolutionsbelow3 km

wouldneedfurtherresearch,whichisoutofthescopeofthispaper.

Therefore,wearecarefulinterpretingthe1kmresolutionresults

andfocusourinterpretationusingthe3kmresolutiontomographic

image.

4.2. Insightsonthegeothermalfields

TheretrievedS-wavevelocityvariationscanoccurasaresult

of wave propagation through different geological media,

spe-cifictectonicfeaturesorrockstate(solid/melting/partialmelting).

Conditions such as crack alignment (isotropic or anisotropic

faultswarms),composition(e.g.relationbetweenminerals,shale

content,fluid content),saturation(porosity,permeability, water

content),pressureandtemperaturedeterminethespeedofseismic

waves(Biot,1956;Gassmann,1951).

Themaincontributionduetoeffectivepressurecanbeobserved

inFig.12,thedeeperstructuresshowhighervelocities(e.g.mean

velocity

v

effectivepressure.However,that isnotalwaysthecase. Asthe

effectivepressureisdefinedbythedifferencebetweenthe

confin-ingpressureandtheporepressure,theporepressuredetermines

theeffectivepressureatthesamedepths.Thisisunderthe

assump-tionofthesameconfiningpressureforthesamedepths,whichis

notthecasefortheWRPtectonicsetting(seedescriptionin

Sec-tion1),asexpectedinatectonicallyandmagmaticallyactivesites.

Inthestudyareatherearethreezonesoffaultswarms(Cliftonand

Kattenhorn,2006):1.withintheReykjanes(R)geothermalfield,2.

withintheEldvörp(E)andSvartsengi(S)geothermalfields3.North

ofReykjanes(R)andEldvörp(E).

Ontopoftheaforementionedcontributions,ingeothermalareas

itiswellknownthatalterationinthehydrothermalsystems,

salin-ity,shalecontent,andwatersaturation(O’ConnellandBudiansky,

1974)interferewiththewavespeed,aswellaswithchangesin

therockstate(solid/melting/partialmelting).Variations influid

contentand temperaturearelikely tobethemostdetermining

parameters for velocityanomaliesin geo- and hydro-thermally

activeareas (Nakajima et al.,2001b), andwhile partialmelting

areasmightnotexist,inclusionsfilledwithH2O(Nakajimaetal.,

2001a)seemreasonabletoexistinReykjanesPeninsula.Therefore,

thepresentstudyshouldbeinterpretedtogetherwithconstraints

fromothermeasurements.

Thederived velocityanomaliesusing 3kmresolutiondo not

covertheReykjanesgeothermalfieldpreventingusfromany

pos-sibleinterpretation.Withthe1kmresolution,theobservedlocal

low-velocityanomaliesattheReykjanesgeothermalfieldmatch

thelocationoftheintensivelyexploitedpartof thegeothermal

reservoir.Thiscouldpotentiallyindicatethattheobserved

low-velocity anomaliesmight be related with a heat source, water

inclusionorboth(O’ConnellandBudiansky,1974;Mavko,1980).

ThesamelocationisdescribedbyFriðleifssonetal.(2018)asthe

areaofup-flowtargetedbyIDDP-2interpretedashotterandmore

permeable.Andaresistivitymodelbasedon3DinversionofMT

data(Karlsdóttiretal.,2018)indicatethattheIDDP-02wellwas

drilledintoalow resistivityanomaly, whichcoincides withthe

observedlowS-wavevelocitiesinFig.12(redarrowinatodatthe

Fig.11.Phase-velocitymapsaftertomographicinversionusingtwofrequencies(0.28ontheleftand0.38Hzontheright),forthefourtestedspatialresolutionswithedited colorbarbetween−10and10percentagedeviationfromtheaveragevelocity(V0).Theresultsareinterpolatedtothesamegridcellsizeforfigurecomparison.Greentriangles

identifytheseismicstationsandthesquaresshowthelocationofthegeothermalfieldsidentifiedin1.Theblackpolygonoutlinestheareawherethereisenoughray-path coveragetobetterrecreatethesimulatedcheckersofFig.9.Redellipsesidentifysimilarlow-velocityanomaliesatdifferentresolutionsinthesamefrequencyband.Blues circlesidentifysimilarhigh-velocityanomaliesatdifferentresolutionsinthesamefrequencyband.

Fig.12.3DS-wavepercentagevelocityvariationsfromtheaveragevelocityperdepth(V0)between2and6kmdepth.Leftcolumnshowstheresultsfrom1kmandright

column3kmspatialresolutiontomographicinversionforalltheidentifieddepths.Weusetheconventionofredforlowandblueforhighvelocityanomalies.Blackdots indicatethelocationsoftheearthquakesdetectedandprovidedbyIMO(IcelandicMeterologicalOffice)from1993to2015.Bluelineslocatethefaultsystemwithinthe definedpolygon.GreentrianglesidentifytheseismicstationsandthesquaresshowthelocationofthegeothermalfieldsidentifiedinFig.1.Redcirclesandarrowscorrespond toareaswhereanomaliesmatchbetweenresolutions,whilebluecirclesandarrowsshowwherenomatchisobserved.

Fig.13.Tomographicresultsofallfrequenciesfor3kmresolution.Thepolygonoutlinestheareawherethereisenoughray-pathcoverage9.Inthisfigurethecolorscale showsvariationsfrom−15%to15%fromtheaveragevelocityperfrequency.

geothermalmanifestations (e.g.theBlueLagoon,locatedwithin

redsquareof Svartsengigeothermalfieldin Fig.1).Atthe

Eld-vörpgeothermalfieldwedonotidentifyanylow-velocityanomaly

withinthestudieddepths,butnorthofEldvörpfieldwedetectlow

speedanomaliesbetween3and5kmdepthintheareaofhigh

frac-turedensity.Franzson(1987)suggeststhatthedensefissureswarm

aroundEldvörpprovidesdownflowchannelsofcold

groundwa-ter,whilewithinthecentralpartofthe“swarm”anupflowofthe

hydrothermalsystemoccursalongthesamefractures.Thederived

3DS-wavecouldeventuallyhelptosupportthishypothesiswitha

detailedcomparisonwiththeexistingboreholedatafromthearea.

Ingeneral,thereisanincreaseoflowvelocitiesfromwestto

eastindicatingthatthetemperaturemaybehighertowardsthe

interiorofthepeninsula.Between4and6kmdepthahigh-velocity

anomalyseemstoseparatetheReykjanesgeothermalfieldfromthe

remainingfieldsofthepeninsula.ThehighS-wavevelocity

struc-turethatmightindicatethepresenceofsolidifiedrocksthathave

beencooleddownpossiblyduetotheproximityoftheocean.

4.3. Futureapplications

Thisstudyshowsthepotentialofambientnoisetomography

asacomplementaryreservoircharacterizationtoolforfield

opera-tions.Fromaseismicnetworkdesignpointofview,theperformed

resolutiontestscanalsoelucidateontheoptimizationofseismic

stationlocationsbyusingthemethodologyofToledoetal.(2018)

extendedforambientnoise.Inasimilarmanner,ourresultsalso

showpotentialtocomplementtheinterpretationofdeformation

studiesoverthesameareaby(e.g.)1.interpretinghowthe

esti-matedhorizontaldisplacementsofHreinsdóttiretal.(2001)or3D

surfacemotionofGudmundssonetal.(2002)canbeobservedby

spatialchangesinS-wavevelocityfieldasaresultofshearstrain;

2.constrainingthesolutionsonthelocalsourcesofman-derived

subsidencederivedbyKeidingetal.(2010)andParksetal.(2018).

with3kmoflateralresolutionsand1km ofverticalresolution.

However,werefrainfrominterpretingthe1kmresolution

tomog-raphyasrealeffectsastheassumptioninherentinthestraightray

approximationmightcausedeteriorationinqualityoftheresults

forhigherresolutions.Althoughtheusedseismicnetworkwasnot

designedforANTonly,itaccommodatesdifferentseismic

imag-ingtechniquesandtheresultsunderlinethecapabilitiesofANSIin

Iceland.

Acknowledgments

The research leading to the results in this manuscript has

received funding from the EC Seventh Framework Programme

undergrantagreementNo.608553(ProjectIMAGE).Theauthors

wouldlike tothankISOR,H.S.Orka, theIcelandMeteorological

Office(IMO),theGeophysicalinstrumentalPoolofPotsdamand

theDEPAS(DeutscheGerätePoolfürAmphibischeSeismologie)for

theirworkingatheringandprovidingustheseismicdataforthis

study.Wewouldliketothanktheanonymousreview,Dr.Kasper

vanWijkandProf.Gudmundssonfor theirinsightfulcomments

andconstructivecriticismduringthereviewprocess.Finally,the

authorswouldliketothankPallEinarsonandKristján

Sæmunds-sonfortheirimensegeologicalcontributionsinIceland,forwhich

somearenicelydescribedinVoightetal.(2018).

AppendixA.

InFig.13weshowthetomographicresultsforeachfrequency

(stepsof0.02Hz)usedfortheinversionfromfrequencytodepth.

Asitcanbeobserved,theanomaliesidentifiedforeachfrequency

aresmoothwhencomparedwiththeadjacentfrequencies.

Consid-eringthateachfrequencyisinvertedindependently,theobserved

References

Albertsson,A.,Bjarnason,J.“O.,Gunnarsson,T.,Ballzus,C.,Ingason,K.,etal., 2003.TheIcelandDeepDrillingProject:fluidhandling,evaluation,and utiliza-tion.,pp.000599234.

Ardhuin,F.,Stutzmann,E.,Schimmel,M.,Mangeney,A.,2011.Oceanwavesources ofseismicnoise.J.Geophys.Res.:Oceans116.

Árnason,K.,Eysteinsson,H.,Hersir,G.P.,2010.Joint1DinversionofTEMandMT dataand3DinversionofMTdataintheHengillarea,SWIceland.Geothermics 39,13–34.

Árnason,K.,etal.,2007.AstudyoftheKraflavolcanousinggravity,microearthquake andMTdata.,pp.001045504.

Benediktsdóttir,Á.,Gudmundsson,Ó.,Brandsdóttir,B.,Tryggvason,A.,2017. Ambi-entnoisetomographyofEyjafjallaj”okullvolcano,Iceland.J.Volcanol. GeothermalRes.347,250–263.

Bensen,G.,Ritzwoller,M.,Barmin,M.,Levshin,A.,Lin,F.,Moschetti,M.,Shapiro, N., Yang,Y., 2007. Processingseismic ambient noise data to obtain reli-ablebroad-bandsurfacewavedispersionmeasurements.Geophys.J.Int.169, 1239–1260.

Benz,H.,Chouet,B.,Dawson,P.,Lahr,J.,Page,R.,Hole,J.,1996.Three-dimensionalP andSwavevelocitystructureofRedoubtVolcano,Alaska.J.Geophys.Res.:Solid Earth101,8111–8128.

Biot,M.A.,1956.Theoryofpropagationofelasticwavesinafluid-saturatedporous solid.II.Higherfrequencyrange.J.AcousticalSoc.Am.28,179–191.

Bjarnason,I.T.,Menke,W.,Flóvenz,Ó.G.,Caress,D.,1993.Tomographicimageof themid-AtlanticplateboundaryinsouthwesternIceland.J.Geophys.Res.:Solid Earth98,6607–6622.

Blanck,H.,Jousset,P.,Hersir,G.P.,Ágústsson,K.,Flóvenz,Ó.G.,2019.Analysisof 2014-2015on-andoff-shorepassiveseismicdataontheReykjanesPeninsula, SWIceland.J.Volcanol.GeothermalRes.

Claerbout,J.F.,1968.Synthesisofalayeredmediumfromitsacoustictransmission response.Geophysics33,264–269.

Clifton,A.E.,Kattenhorn,S.A.,2006.Structuralarchitectureofahighlyoblique diver-gentplateboundarysegment.Tectonophysics419,27–40.

Darnet,M.,Wawrzyniak,P.,Coppo,N.,Nielsson,S.,Schill,E.,Friðleifsson,G.Ó.,2018. MonitoringgeothermalreservoirdevelopmentswiththeControlled-Source Electro-Magneticmethod-AcalibrationstudyontheReykjanesgeothermal field.Elsevier.

DeRidder,S.,Biondi,B.,Nichols,D.,2015.Elliptical-anisotropiceikonalphase veloc-itytomography.Geophys.Res.Lett.42,758–764.

Einarsson,P.,2008.Plateboundaries,riftsandtransformsinIceland.J“okull 58,35–58.

Einarsson,P.,Þorbjarnarson,S.,B”ottger,M.,2002.Faultsandfracturesofthe SouthIcelandseismiczonenearÞjórsá.Landsvirkjun.

Einarsson,P.,Saemundsson,K.,1987.EarthquakeEpicenters1982-1985and Vol-canicSystemsinIceland:UpptokJardskjalfta1982-1985OgEldstodvakerfia Islandi.Menningarsjodur.

Elders,W.,Friðleifsson, G.,Albertsson,A., 2014.Drillinginto magmaandthe implicationsoftheIcelandDeepDrillingProject(IDDP)forhigh-temperature geothermalsystemsworldwide.Geothermics49,111–118.

Elders, W.A., Friðleifsson, G.Ó., Zierenberg, R.A., Pope, E.C., Mortensen, A.K., Guðmundsson,Á.,Lowenstern,J.B.,Marks,N.E.,Owens,L.,Bird,D.K.,etal.,2011. OriginofarhyolitethatintrudedageothermalwellwhiledrillingattheKrafla volcano,Iceland.Geology39,231–234.

Erlendsson,P.,Einarsson,P.,1996.TheHvalhnúkurFault,astrike-slipfaultmapped withintheReykjanesPeninsulaobliquerift,Iceland.SeismologyinEurope.XXV ESCGeneralAssembly,ReykjavıkIceland,pp.9–14,September.

Eysteinsson,H.,Hermance,J.F.,1985.Magnetotelluricmeasurementsacrossthe easternneovolcaniczoneinsouthIceland.J.Geophys.Res.:SolidEarth90, 10093–10103.

Franzson,H.,1987.TheEldvorpHigh-TemperatureArea,SW-Iceland.Geothermal GeologyofFirstExplorationWell.Proceedingsofthe9thNewZealand Geother-malWorkshop,179–185.

Friðleifsson,G.,Elders,W.,Albertsson,A.,2014a.TheconceptoftheIcelanddeep drillingproject.Geothermics49,2–8.

Friðleifsson,G.,Elders,W.A.,Zierenberg,R.A.,Stefánsson,A.,Fowler,A.P., Weisen-berger,T.B.,Harðarson,B.S.,Mesfin,K.G.,2017a.Theicelanddeepdrillingproject 4.5kmdeepwell,iddp-2,intheseawater-rechargedreykjanesgeothermalfield inSWIcelandhassuccessfullyreacheditssupercriticaltarget.ScientificDrilling 23,1–12.

Friðleifsson, G.,Albertsson,Ó.A., 2000.Deepgeothermaldrilling atReykjanes Ridge:opportunityforaninternationalcollaboration.ProceedingsoftheWorld GeothermalCongress,F7-5.

Friðleifsson, G.Ó., Elders,W.A., Zierenberg,R., Steafánsson,A., Sigurðsson,Ó., Gíslason,Þ., Weisenberger,T.B.,Harðarson,B.S.,Mesfin,K.G., 2017b.ICDP supportedcoringinIDDP-2atReykjanes-theDEEPEGSdemonstratorin Iceland-Supercriticalconditionsreachedbelow4.6kmdepth.EGUGeneralAssembly ConferenceAbstracts,14147,volume19.

Friðleifsson,G.Ó.,Elders,W.A.,Zierenberg,R.A.,Fowler,A.P.,Weisenberger,T.B., Mesfin,K.G.,Sigurðsson,Ó.,Níelsson,S.,Einarsson,G.,Óskarsson,F.,etal.,2018. TheIcelandDeepDrillingProjectatReykjanes:Drillingintotherootzoneofa blacksmokeranalog.J.Volcanol.GeothermalRes.

Friðleifsson,G.Ó.,Sigurdsson,Ó.,Þorbj“ornsson,D.,Karlsdóttir,R.,Gíslason, Þ.,Albertsson,A.,Elders,W.,2014b.PreparationfordrillingwellIDDP-2at Reyk-janes.Geothermics49,119–126.

Froment,B.,Campillo,M.,Roux,P.,Gouedard,P.,Verdel,A.,Weaver,R.L.,2010. Estimationoftheeffectofnonisotropicallydistributedenergyontheapparent arrivaltimeincorrelations.Geophysics75,SA85–SA93.

Gasperikova,E.,Rosenkjaer,G.K.,Arnason,K.,Newman,G.A.,Lindsey,N.J.,2015. ResistivitycharacterizationoftheKraflaandHengillgeothermalfieldsthrough 3DMTinversemodeling.Geothermics57,246–257.

Gassmann,F.,1951.Elasticwavesthroughapackingofspheres.Geophysics16, 673–685.

Golub,G.H.,Heath,M.,Wahba,G.,1979.Generalizedcross-validationasamethod forchoosingagoodridgeparameter.Technometrics21,215–223.

Green,R.G.,Priestley,K.F.,White,R.S.,2017.Ambientnoisetomographyreveals uppercrustalstructureofIcelandicrifts.EarthPlanet.Sci.Lett.466,20–31. Gudmundsson,S.,Sigmundsson,F.,Carstensen,J.M.,2002.Three-dimensional

sur-facemotionmapsestimatedfromcombinedinterferometricsyntheticaperture radarandGPSdata.J.Geophys.Res.:SolidEarth107,ETG–ETG13.

Guðnason,E.Á.,etal.,2014.Analysisofseismicactivityonthewesternpartofthe ReykjanesPeninsula.SWIceland,December2008-May2009.

Halliday,D.,Curtis,A.,2008.Seismicinterferometry,surfacewavesandsource dis-tribution.Geophys.J.Int.175,1067–1087.

Haney,M.M.,Tsai,V.C.,2015.Nonperturbationalsurface-waveinversion:ADix-type relationforsurfacewaves.Geophysics80,EN167–EN177.

Hersir,G.P.,Árnason, K.,Vilhjálmsson,A.M.,Saemundsson,K.,Ágústsdóttir,Þ., Friðleifsson,G.Ó.,2018.Kr‘ysuvíkhightemperaturegeothermalareainSW Iceland:Geologicalsettingand3Dinversionofmagnetotelluric(MT)resistivity data.J.Volcanol.GeothermalRes.

Hersir,G.P.,Bj”ornsson,A., Pedersen,L.B.,1984.Magnetotelluricsurvey acrosstheactivespreadingzoneinsouthwestIceland.J.Volcanol.Geothermal Res.20,253–265.

Hólmgeirsson,S.,Gudmundsson,A.,Pálsson,B.,Bóasson,H.,Ingasson,K., Thorhalls-son,S.,2010.DrillingoperationsofthefirstIcelanddeepdrillingwell(IDDP). ProceedingsoftheWorldGeothermalCongress2010,10.

Hreinsdóttir,S.,Einarsson,P.,Sigmundsson,F.,2001.Crustaldeformationatthe obliquespreadingReykjanesPeninsula,SWIceland:GPSmeasurementsfrom 1993to1998.J.Geophys.Res.:SolidEarth106,13803–13816.

Jakobsson,S.P.,Jónsson,J.,Shido,F.,1978.PetrologyofthewesternReykjanes penin-sula,Iceland.J.Petrol.19,669–705.

Jeddi,Z.,Gudmundsson,O.,Tryggvason,A.,2017.Ambient-noisetomographyof Katlavolcano,southIceland.J.Volcanol.GeothermalRes.347,264–277. Jousset,P.,Ágústsson,K.,Blanck,H.,Metz,M.,Franke,S.,PàllHersir,G.,Bruhn,D.,

Flovenz,Ó.,Friðleifsson,G.,2017.TraveltimeseismictomographyonReykjanes, SWIceland.EGUGeneralAssemblyConferenceAbstracts,17398,volume19. Karlsdóttir,R.,Vilhjálmsson,A.M.,Guðnason,E.Á.,2018.Threedimensional

inver-sionofmagnetotelluric(MT)resistivitydatafromReykjaneshightemperature fieldinSWIceland.J.Volcanol.GeothermalRes.

Keiding,M.,Árnadóttir,T.,Jonsson,S.,Decriem,J.,Hooper,A.,2010.Plate bound-arydeformationandman-madesubsidencearoundgeothermalfieldsonthe ReykjanesPeninsula,Iceland.J.Volcanol.GeothermalRes.194,139–149. Klein,F.W.,Einarsson,P.,Wyss,M.,1973.MicroearthquakesontheMid-Atlantic

PlateBoundaryontheReykjanesPeninsulainIceland.J.Geophys.Res.78, 5084–5099.

Liu,H.-P.,Anderson,D.L.,Kanamori,H.,1976.Velocitydispersionduetoanelasticity; implicationsforseismologyandmantlecomposition.Geophys.J.Int.47,41–58. Martins,J.E.,Ruigrok,E.,Draganov,D.,Hooper,A.,Hanssen,R.F.,White,R.S.,Soosalu, H.,2019.ImagingTorfajökull’sMagmaticPlumbingSystemWithSeismic Inter-ferometryandPhaseVelocitySurfaceWaveTomography.J.Geophys.Res.:Solid Earth124,2920–2940.

Mavko,G.M.,1980.Velocityandattenuationinpartiallymoltenrocks.J.Geophys. Res.:SolidEarth85,5173–5189.

Nakajima,J.,Hasegawa,A.,2003.Tomographicimagingofseismicvelocitystructure inandaroundtheOnikobevolcanicarea,northeasternJapan:implicationsfor fluiddistribution.J.Volcanol.GeothermalRes.127,1–18.

Nakajima,J.,Matsuzawa,T.,Hasegawa,A.,Zhao,D.,2001a.Seismicimagingofarc magmaandfluidsunderthecentralpartofnortheasternJapan.Tectonophysics 341,1–17.

Nakajima,J.,Matsuzawa,T.,Hasegawa,A.,Zhao,D.,2001b.Three-dimensional struc-tureofVp,Vs,andVp/VsbeneathnortheasternJapan:Implicationsforarc magmatismandfluids.J.Geophys.Res.:SolidEarth106,21843–21857. Newman,G.A.,Gasperikova,E.,Hoversten,G.M.,Wannamaker,P.E.,2008.

Three-dimensionalmagnetotelluriccharacterization oftheCosogeothermalfield. Geothermics37,369–399.

Obermann,A.,Lupi,M.,Mordret,A.,Jakobsdóttir,S.S.,Miller,S.A.,2016.3D-ambient noiseRayleighwavetomographyofSnæfellsj“okullvolcano,Iceland.J. Volcanol.GeothermalRes.317,42–52.

O’Connell,R.J.,Budiansky,B.,1974.Seismicvelocitiesindryandsaturatedcracked solids.J.Geophys.Res.79,5412–5426.

Oskooi,B.,Pedersen,L.B.,Smirnov,M.,Árnason,K.,Eysteinsson,H.,Manzella,A., Group,D.W.,etal.,2005.ThedeepgeothermalstructureoftheMid-Atlantic RidgededucedfromMTdatainSWIceland.Phys.EarthPlanet.Interiors150, 183–195.

Pálmason,G.,Sæmundsson,K.,1974.IcelandinrelationtotheMid-AtlanticRidge. Annu.Rev.EarthPlanet.Sci.2,25–50.

Park,C.B.,Miller,R.D.,Xia,J.,1999.Multichannelanalysisofsurfacewaves. Geo-physics64,800–808.

Park,C.B.,Miller,R.D.,Xia,J.,etal.,1998.Imagingdispersioncurvesofsurfacewaves onmulti-channelrecord.SEGExpandedAbstracts,1377–1380,volume17.

imagingof theP-andS-wave velocitystructure andearthquake locations beneathSouthwestIceland.Geophys.J.Int.151,848–866.

Tsai,V.C.,2009.Onestablishingtheaccuracyofnoisetomographytravel-time mea-surementsinarealisticmedium.Geophys.J.Int.178,1555–1564.

Zhou,Y.,Dahlen,F.,Nolet,G.,2004.Three-dimensionalsensitivitykernelsforsurface waveobservables.Geophys.J.Int.158,142–168.