Journal of Geophysical Research: Solid Earth RESEARCHARTICLE Interseismic strain accumulation across the central North 10.1002/2016JB013108 Anatolian Fault from iteratively unwrapped InSAR measurements KeyPoints: •Newiterativephase-unwrapping procedureimprovescoverageof Ekbal Hussain1,Andrew Hooper1,TimJ. Wright1,RichardJ. Walters2,andDavidP.S. Bekaert3,4 InSARmeasurements •Surfacevelocitiesatthecentral 1COMET,SchoolofEarthandEnvironment,UniversityofLeeds,Leeds,UK,2COMET,DepartmentofEarthSciences, NorthAnatolianFaultshoweastward decreaseinsliprate DurhamUniversity,Durham,UK,3JetPropulsionLaboratory,CaliforniaInstituteofTechnology,Pasadena,California,USA, •FaultcreepnearIsmetpasaisreleasing 4FormerlyatCOMET,UniversityofLeeds,Leeds,UK only30-40%oflong-termstraininthe shallowcrust Abstract TheNorthAnatolianFault(NAF)isamajortectonicfeatureintheMiddleEastandisthemost SupportingInformation: activefaultinTurkey.ThecentralportionoftheNAFisaregionofGlobalNavigationSatelliteSystems •SupportingInformationS1 (GNSS)scarcity.Previousstudiesofinterseismicdeformationhavefocusedontheaseismiccreepnearthe townofIsmetpasausingradardataacquiredinasingleline-of-sightdirection,requiringseveralmodeling Correspondenceto: assumptions.WehavemeasuredinterseismicdeformationacrosstheNAFusingbothascendingand E.Hussain, descendingdatafromtheEnvisatsatellitemissionacquiredbetween2003and2010.Ratherthanrejecting [email protected] incorrectlyunwrappedareasintheinterferograms,wedevelopanewiterativeunwrappingprocedure forsmallbaselineinterferometricsyntheticapertureradar(InSAR)processingthatexpandsthespatial Citation: coverage.Ourmethodcorrectsunwrappingerrorsiterativelyandincreasestherobustnessofthe Hussain,E.,A.Hooper,T.J.Wright, unwrappingprocedure.WeremovelongwavelengthtrendsfromtheInSARdatausingGNSSobservations R.J.Walters,andD.P.S.Bekaert(2016), Interseismicstrainaccumulationacross anddeconvolvetheInSARvelocitiesintofault-parallelmotion.Profilesoffault-parallelvelocityreveala thecentralNorthAnatolianFault systematiceastwarddecreaseinfaultslipratefrom30mm/yr(25–34,95%confidenceinterval(CI))to fromiterativelyunwrappedInSAR 21mm/yr(14–27,95%CI)overadistanceof∼200km.Directoffsetmeasurementsacrossthefaultreveal measurements,J.Geophys.Res.Solid Earth,121,doi:10.1002/2016JB013108. faultcreepalonga∼130kmsectionofthecentralNAF,withanaveragecreeprateof8±2mm/yranda maximumcreeprateof14±2mm/yrlocated∼30kmeastofIsmetpasa.Asfaultcreepisreleasingonly 30–40%ofthelong-termstrainintheshallowcrust,thefaultisstillcapableofproducinglarge,damaging Received20APR2016 Accepted27NOV2016 earthquakesinthisregion. Acceptedarticleonline7DEC2016 1.Introduction TheNorthAnatolianFault(NAF)isamajorcontinentalright-lateraltransformfaultlocatedinnorthernTurkey. TogetherwiththeEastAnatolianFault,itfacilitatesthewestwardmotionofAnatolia,caughtintheconver- gencezoneoftheEurasianplatewiththeArabianplate[McKenzie,1972].Sincethe1939M 7.9Erzincan w earthquakeineasternTurkey,theNAFhasrupturedinasequenceoflarge(M >6.7)earthquakeswithadom- w inantwestwardprogressioninseismicity[Barka,1996;Steinetal.,1997].Steinetal.[1997]andHubert-Ferrari etal.[2000]haveinterpretedthissequencetoresultfromstresstransferalongstrike,whereoneearthquake bringstheadjacentsegmentclosertofailure. InordertounderstandtherolethattheNAFplaysinregionaltectonicsandseismichazard,therehavebeen numerousestimatesofthefaultslipratefortheNAFusingpresent-daydeformationmeasuredwithGlobal NavigationSatelliteSystems(GNSS)[e.g.,Straubetal.,1997;Reilingeretal.,2006;Ergintavetal.,2009]oroffset geologicalfeatures[e.g.,Hubert-Ferrarietal.,2002;Puccietal.,2008;Kozacietal.,2009].Therehavealsobeen severalinterferometricsyntheticapertureradar(InSAR)-derivedestimatesofthefaultsliprate,whichhave focusedonthewesternoreasternregionsoftheNAFwheretheInSARcoherenceisbetter[e.g.,Wrightetal., 2001a;Cakiretal.,2005;Waltersetal.,2011;Kanekoetal.,2013;Cakiretal.,2014;Cetinetal.,2014;Waltersetal., ©2016.TheAuthors. Thisisanopenaccessarticleunderthe 2014;CavaliéandJónsson,2014;Hussainetal.,2016]. termsoftheCreativeCommons AttributionLicense,whichpermitsuse, However,sliprateestimatesforthecentralNAFarerelativelypoorlyconstrained,withsparseGNSSdatanorth distributionandreproductioninany ofthisportionofthefault(Figure1)andwideranginggeologicalandgeodeticestimates.Geologicalfault medium,providedtheoriginalworkis properlycited. slipraterangesfromaslowas5mm/yrtoashighas44mm/yr[e.g.,BarkaandHancock,1984;Barka,1992; HUSSAINETAL. INTERSEISMICCENTRALNAF 1 Journal of Geophysical Research: Solid Earth 10.1002/2016JB013108 Figure1.(a)ThecentralsectionoftheNorthAnatolianFault.TheredarrowsarepublishedGNSSvelocitiesfromthe GlobalStrainRateModelproject[Kreemeretal.,2014].Thecoloredsectionsindicatepreviousrupturesalongthissection ofthefault.(b)TheEnvisatsatellitedatatracksusedinthisstudy.Descendingtracksarecoloredinredandascending tracksinblue. Hubert-Ferrarietal.,2002;Kozacietal.,2007;Kozacietal.,2009],whileGNSSstudiesestimatetheslipratefor theregiontoarangeof17–34mm/yr[e.g.,Oraletal.,1993;Noomenetal.,1996;Ayhanetal.,2002;Reilinger etal.,2006]. Shallowaseismicsliponthefaultplane,i.e.,faultcreep,onthecentralportionoftheNAFwasfirstdocumented byAmbraseys[1970],whoobservedincreasingdisplacementsofawallthatwasbuiltacrossthefaultnear thetownofIsmetpasa,overmultipleyears.Ambraseys[1970]estimatedafaultcreeprateof∼20mm/yrfor thetimeperiod1955–1969.Sincethisoriginalinvestigation,thefaultcreephasbeenthefocusofnumerous geodeticstudies[e.g.,Cakiretal.,2005;Kutogluetal.,2010;Karabacaketal.,2011;Ozeneretal.,2013;Cetinetal., 2014].Cetinetal.[2014]suggestedthatthefaultcreepratehasbeendecayingsincethefirstmeasurementsin 1970toacurrentsteadystatevalueof∼6–8mm/yr.MostpreviousInSARstudiesinthisregionhaveonlyused satellitedatafromasingle-lookdirection,e.g.,theuseofdescendingEnvisatdatabyCakiretal.[2005]and Cetinetal.[2014].Kanekoetal.[2013]usedacombinationofascendingtracksfromtheALOSsatelliteandone descendingframefromEnvisattrack207,limitingtheirobservationalperiodto2007–2011.Theysuggested thataseismiccreepatarateof∼9mm/yrislimitedtotheupper5.5–7kmofthecrust,whichexhibitsvelocity strengtheningfrictionalbehavior. Recently,Roussetetal.[2016]usedhigh-resolutionCOSMO-SkyMedsatellitedataspanningthetimewindow betweenJuly2013andMay2014toshowevidenceofperiodsofelevatedfaultcreepspanningamonth withtotalslipof20mm,indicatingthatepisodiccreepeventsmaybeanimportantmechanismproducing aseismicslip. InthisstudyweuseamorecompletedatasetcoveringtheentirecentralNAFinbothascendinganddescend- inggeometriesandspanningthe∼8yeartimewindowbetween2003and2010.Weremovelongwavelength trendsfromtheInSARdatausingpublishedGNSSvelocities[Kreemeretal.,2014]anddeconvolvetheInSAR line-of-sightvelocitiesintofault-parallelandverticalmotion. Weusesimpleelasticdislocationmodelstoestimategeodeticfaultslipratesandlockingdepthsandinves- tigatethespatialvariationoffaultcreepalongthecentralNAF.Wealsodevelopandapplyanewiterative unwrappingalgorithmthatminimizesunwrappingerrorsduringtheInSARprocessing. HUSSAINETAL. INTERSEISMICCENTRALNAF 2 Journal of Geophysical Research: Solid Earth 10.1002/2016JB013108 Table1.DataCoverageforEachEnvisatTrackUsedinThisStudy Track Geometry TimeSpana No.ofImages TotalIntsCreated IntsUsed 250 Descending 2003/12/12–2010/7/23 38 115 59 479 Descending 2003/12/28–2010/7/4 30 90 50 207 Descending 2004/1/13–2010/9/28 40 88 53 436 Descending 2003/7/3–2010/3/18 36 96 65 28 Ascending 2004/7/28–2010/7/7 14 30 21 71 Ascending 2004/1/3–2009/8/29 19 48 29 343 Ascending 2004/6/10–2010/4/15 14 27 20 aDatesareformattedasyear/month/day. 2.InSARProcessing Ourdatasetconsistsof191Envisatimagesfromfourdescendingtracks(250,479,207,and436)andthree ∘ ∘ ascendingtracks(28,71,and343)(Figure1b).TogetherthesecoverthecentralNAFbetween31.5 Eand35 E, andspanthetimeinterval2003–2010.DetailsoftheprocesseddataforeachtrackaregiveninTable1. WefocustheEnvisatimagesusingROI_PAC[Rosenetal.,2004]andusetheDORISsoftware[Kampesetal.,2003] toconstruct494interferograms.Foreachtrackweproducearedundantconnectednetworkofinterferograms whileminimizingthetemporalseparationbetweenacquisitionsandthespatialseparationofthesatellite(the perpendicularbaseline)(FigureS1inthesupportinginformation).Wecorrecttopographiccontributionsto theradarphaseusingthe90mSRTMDigitalElevationModel[Farretal.,2007]andaccountfortheknown oscillatordriftforEnvisataccordingtoMarinkovicandLarsen[2013].Weunwraptheinterferometricphase usinganewiterativeunwrappingprocessdescribedinsection3. WeapplytheStaMPS(StanfordMethodforPersistentScatterers)smallbaselinetimeseriestechnique[Hooper, 2008;Hooperetal.,2012]toremoveincoherentpixelsandreducethenoisecontributiontothedeformation signal,byselectingonlythosepixelsthathavelowphasenoiseonaverageinthesmallbaselineinterferograms usedintheanalysis. Theatmosphericcontributionisoftenthelargestsourceoferrorinradarinterferograms[e.g.,Doinetal.,2009; Waltersetal.,2013;Jolivetetal.,2014;Bekaertetal.,2015a].Tomitigatethisweestimatedatropospherecorrec- tionusingauxiliarydatafromtheERA-Interimglobalatmosphericmodelreanalysisproduct[Deeetal.,2011]. WeusetheTRAIN(ToolboxforReducingAtmosphericInSARNoise)softwarepackage[Bekaertetal.,2015c]to correcteachindividualinterferogramfortroposphericnoise.Afterremovingaplanarphaserampfromeach interferogram,theERA-Icorrectionreducesthestandarddeviationofourtracksby8%onaverage.Theaver- agereductioninstandarddeviationissmallaftercorrection,implyingthatsomeresidualatmosphericsignals remainintheinterferogramsaftertheERA-Icorrection.Theaveragereductioninstandarddeviationforeach trackare10%fortrack207,1%fortrack250,2%fortrack436,12%fortrack479,10%fortrack28,16%for track71,and6%fortrack343(FiguresS2andS3). Ourfinalredundantsmallbaselinenetworksconsistofatotalof297interferogramsovertheseventracks (FigureS1).Weusethesenetworkstocalculatetheaverageline-of-sight(LOS)velocitymapforeachtrack. Anynontectoniclongwavelengthsignals(>100km),includingthoseduetoorbitalerrors,areeffectively removedfromeachtrackwhentheInSARline-of-sight(LOS)velocitiesaretransformedintoaEurasia-fixed GNSSreferenceframe(detailsinsection4).Theuncertaintiesonthefinalvelocityforeachpixelarecalcu- latedusingbootstrapresampling[EfronandTibshirani,1986]andarepresentedatthe1sigmalevelinthe followingwork. WecalculatetheLOSvariance-covariancematrixofthenoiseforeachInSARtrackbycomputingtheaverage radialcovarianceversusdistance(autocorrelation)usingthevelocitiesina50kmby50kmregion∼250kmto thesouthofthefault.Thisregionisassumedtohavenotectonicdeformationandcontainsonlyatmospheric noise.Wefitanexponentialcovariancefunction[e.g.,LohmanandSimons,2005;Parsonsetal.,2006],C(r),as C(r)=𝜎2e−𝜆r, (1) HUSSAINETAL. INTERSEISMICCENTRALNAF 3 Journal of Geophysical Research: Solid Earth 10.1002/2016JB013108 Table2.TheCenterofthe50kmby50kmRegionUsedtoEstimatetheNoise CovarianceFunctionParameters Track Center(Lon,Lat) Variance,𝜎2(mm/yr)2 CharacteristicLength,𝜆(km) 207 33∘E,39.5∘N 8.91 53 250 31.75∘E,39.5∘N 4.95 27 436 34∘E,39.5∘N 3.91 22 479 32.5∘E,39.5∘N 2.88 10 28 34.5∘E,39.5∘N 6.12 25 71 33.2∘5E,39.5∘N 4.00 19 343 32.5∘E,39.5∘N 1.00 4 whereweestimatethevariance(𝜎2)andthecharacteristiclength(𝜆),whichgivethespatialcorrelationof noiseasafunctionofdistancebetweenpixels(r).Ourvaluesforeachtrackandthecenteroftheregion usedtocalculatethecovariancefunctionareshowninTable2.Thesecovariancesareusedinsection5when modelingthehorizontalvelocitiesandfaultcreeprates. 3.IterativePhaseUnwrapping 3.1.MethodDescription Phase unwrapping is the process of recovering continuous phase values from phase data that are mea- suredmodulo2𝜋radians(wrappeddata)[GhigliaandPritt,1998].Original2-Dphase-unwrappingalgorithms unwrappedthephaseofeachindividualinterferogramindependently[e.g.,Goldsteinetal.,1988;Costantini, 1998; ZebkerandLu, 1998]. However, a time series of selected interferogram pixels can be considered a 3-D data set, the third dimension being that of time.HooperandZebker [2007] showed that treating the unwrappingproblemasone3-Dproblemasopposedtoaseriesof2-Dproblemsleadstoanimprovement intheaccuracyofthesolutioninasimilarwaytowhich2-Dunwrappingprovidesanimprovementover one-dimensionalspatialmethods. Fully 3-D phase-unwrapping algorithms commonly assume that the phase difference between neighbor- ingpixelsisgenerallylessthanhalfaphasecycle(2𝜋radians)inalldimensions[HooperandZebker,2007]. However,duetoatmosphericdelays,InSARsignalsareeffectivelyuncorrelatedintime,violatingthisassump- tion.Otherunwrappingalgorithmsrequiretheassumptionofatemporalparametricfunction,suchasalinear phaseevolutionintime[Ferrettietal.,2001],tounwrapthephasesignals. Thestandardunwrapping algorithmused intheStanford Methodfor PersistentScatterers(StaMPS) soft- ware[Hooper,2010]usestheactualphaseevolutionintimetoguideunwrappinginthespatialdimension without assuming a particular temporal evolution model. The phase difference between nearby pixels (double-differencephase)isfilteredintimetogiveanestimateoftheunwrappeddisplacementphaseforeach satelliteacquisitionandanestimateofthephasenoise.Thisisusedtoconstructprobabilitydensityfunctions foreachunwrappeddouble-differencephaseineveryinterferogram.Anefficientalgorithm(SNAPHU)[Chen andZebker,2000,2001]thensearchesforthesolutioninspacethatmaximizesthetotaljointprobability,i.e., minimizesthetotal‘cost’. Foraconnectednetworkofsmallbaselineinterferograms,thephaseunwrappingofindividualinterferograms canbecheckedfornetworkconsistencybysummingthephasearoundclosedinterferometricloops[e.g.,Pepe andLanari,2006;Biggsetal.,2007;Cavaliéetal.,2007;Jolivetetal.,2011](Figure2).Inthestandardunwrapping approachusedinStaMPS,anyinterferogramsidentifiedtohavelargeunwrappingerrorsareremovedfrom thesmallbaselinenetwork,whichcanresultinlossofinformationand/orreductioninnetworkredundancy. NotethatsomeotherInSARpractitioners[e.g.,Biggsetal.,2007;Wangetal.,2009;Waltersetal.,2011]generally donotdropbadlyunwrappedinterferograms,butattempttocorrectunwrappingerrorsbymanuallyadding integermultiplesof2𝜋tobadlyunwrappedregionsofpixels.However,thisisatime-consumingprocess. In our method, we iterate the standard StaMPS unwrapping procedure while calculating the sum of the unwrappedphasearoundclosedloopsforeverypixelineveryinterferogram,usingthefollowingequation: ∑n−1 UW{𝜙 −𝜙}+𝜖=0, (2) (i+1)modn i i=0 HUSSAINETAL. INTERSEISMICCENTRALNAF 4 Journal of Geophysical Research: Solid Earth 10.1002/2016JB013108 Figure2.Asimpleinterferometricloopconsistingofthreeacquisitions(redpoints)withphase𝜙 .Theinterferograms 0∶2 aredenotedbythebluelinesandarethedifferenceinphasefortwoacquisitions.UWistheStaMPSunwrapping operator,seetextfordetails.Foreverypixelunwrappedcorrectlyineachinterferogramthephasesumaroundtheloop isequaltozero,i.e.,UW(𝜙 −𝜙 )+UW(𝜙 −𝜙 )+UW(𝜙 −𝜙 )=0. 1 0 2 1 0 2 whereUWistheStaMPSunwrappingoperator,nisthenumberofinterferogramsonthepatharoundaninter- ferometricloop,(𝜙 −𝜙)aretheinterferometricphasevaluesofapixelintheinterferogramscreatedbycal- i+1 i culatingthephasedifferencebetweenimagei+1andirelativetoareferencepoint,and𝜖istheerrorterm.The referencepointischosentobenorthofthefaultforalltracks.Anypixelssatisfyingtherequirementof|𝜖|<1 radaredefinedas“error-freepixels”andareassumedtobecorrectlyunwrapped.Anerrortermisneeded becausetheinterferogramsaremultilookedbeforeunwrapping,andsowedonotexpecttohaveperfectclo- surearoundeachinterferometricloop.Using𝜖=1isreasonableasitiswellbelowthe2𝜋radiansrequiredto produceunwrappingerrorsandallowsforasmallamountofclosureerrorintroducedbythenonlinearnature ofmultilooking.Inourtestssetting𝜖to0.5madenosignificantimpactontheacceptancerates. In each iteration, we keep all unwrapping parameters fixed (such as the number of interferograms and filtering)butassumethatpixelsidentifiedaserror-freeinthepreviousiterationarelikelyunwrappedcorrectly andapplyahighcosttochangingthephasedifferencebetweenthesepixelsinthenextiteration.TheStaMPS unwrappingalgorithmusesthedouble-differencephaseevolutionintimetocalculateaprobabilitydensity functionofunwrappedphaseforeachpixelpairineachinterferogram.Forinterferogramswherebothpixels inapairareidentifiedasunwrappedcorrectly,wesettheweightingto100timesthoseoftheotherinterfer- ograms,toeffectivelyensurethattheevolutionintimeisfixed.Inthisway,theiterativeunwrappingmethod usestheerror-freepixelsasaguidetounwrappingtheregionsthatcontainedunwrappingerrorsinprevious iterations. López-Quirozetal.[2009]describeaprocesswhereunwrappingisiteratedontheresidualinterferogramafter theremovalofanestimateofthedeformationsignalwhileourtechniqueiteratestheStaMPSunwrapping procedureontheactualinterferometricphase. 3.2.TestingtheIterativeUnwrappingProcedure We tested the new algorithm on data from Envisat descending track 207, which covers a region roughly 100kmby400kmincentralTurkey(Figure1b).Eachiterationconsistsofthefollowingsteps:runningthe StaMPSunwrappingalgorithm,determiningthepixelsunwrappedcorrectlyineachinterferogramusingthe methoddescribedaboveandinAppendixA,applyingahighcosttounwrappingacrossthesepixels,and rerunningtheunwrappingalgorithmagain.Weiteratethisprocedure30times.Theresultsfromstandard unwrappingdoesnotchangeasnomodificationsaremadetoitsinputsandisrepresentedbythestraight lineindicatingnochangeinthenumberoferror-freepixelsperiteration.Figure3showsthatthepercentage oferror-freepixelsintheentiresmallbaselinenetworkincreasessharplywiththefirsteightiterationsfrom 70%to83%,reachingamaximumof84%after30iterations;meaningthattherearesomeunwrappingerrors themethodisunabletofix.Thisisalsoevidentfromtheindividualinterferograms(Figure4),whichshowthis samerapidincreaseinthepercentageoferror-freepixelsfollowedbyaplateau.Itisclearthattherearesome unwrappingerrorsthatcannotbecorrected(bluecolorsinFigure5)usingtheiterativemethod.However, theiterativeproceduregreatlyreducesthetotalnumberofunwrappingerrorsandthusincreasestheInSAR coveragewhileminimizingerrors. Aftereightiterationsthepercentageoferror-freepixelsincreasedfrom90%to94%fortrack250,from65% to80%fortrack436,from92%to95%fortrack479,from83%to87%fortrack343,from71%to77%fortrack 28,andfrom91%to93%fortrack71. HUSSAINETAL. INTERSEISMICCENTRALNAF 5 Journal of Geophysical Research: Solid Earth 10.1002/2016JB013108 Figure3.Totalpercentageofpixelsinthesmallbaselinenetworkfordescendingtrack207thatwereidentifiedas closed,i.e.,correctlyunwrapped,usingouriterativeunwrappingprocedure(blue)andthestandardunwrapping(red) algorithm.Thereisarapidincreaseinthenumberoferror-freepixelsforthefirsteightiterationsafterwhichitreachesa plateau.Asnomodificationismadetotheinputoftheunwrappingalgorithm,thereisnochangeforeachiterationof thestandardunwrappingalgorithm. 4.InterseismicVelocityFieldAcrosstheCentralNAF Toinvestigatethepatternofinterseismicstrainaccumulationalongthefault,wedecomposeourfullInSAR velocityfieldintothefault-parallelandfault-perpendicularcomponentsofmotion.Followingthemethod describedinHussainetal.[2016],wedothisfirstbyresamplingourInSARLOSvelocities(Figure6)ontoa 1kmby1kmgridencompassingthespatialextentofallourtracks.Weuseanearestneighborresampling techniqueincludingonlythosepersistentscattererpixelswithanearestneighborwithin1kmofthecenter ofeachgridpoint.WereferenceeachtracktoaEurasia-fixedGNSSreferenceframebyfirstaveragingthe InSARvelocitiesthatfallina1kmradiusaroundeveryGNSSstationwithintheboundariesofeachInSAR track.WeprojecttheGNSSvelocitiesintothelocalsatellitelineofsightandcalculatethedifferencefromthe Figure4.Changesinthepercentageoferror-freepixels(correctlyunwrappedpixels)periterationshownforselected interferograms.Inbluearethechangesfortheiterativeunwrappingalgorithm,whileredindicatesthestandard unwrapping. HUSSAINETAL. INTERSEISMICCENTRALNAF 6 Journal of Geophysical Research: Solid Earth 10.1002/2016JB013108 Figure5.Evolutionofthenumberoferror-freepixels(correctlyunwrappedpixels)periterationshownforinterferogram29.Error-freepixelsareidentifiedinred, whilepixelsthatdidnotclose,i.e.,haveunwrappingerrors,areinblue.TheunwrappedphaseforeachiterationisshowninFigureS7inthesupporting information. InSARvelocities.TheverticalcomponentoftheGNSSvelocitiesarenotavailableontheGlobalStrainRate Modelwebsite.Ergintavetal.[2009]showedthattheverticalGNSScomponentissmallandverynoisyover westernTurkey;therefore,weonlyusethehorizontalvelocitiesinouranalysis.Wedeterminethebestfitplane throughtheresidualvelocitiesandremovethisfromtheInSARvelocitiestotransformtheLOSvelocitiesinto aEurasia-fixedGNSSreferenceframe.Thisprocedureisdoneseparatelyforeachtrack. Toestimatetheuncertaintiesinthedata,wecalculatetheRMSresidualinhorizontalvelocitiesintheover- lappingareasbetweenneighboringtracksassumingnegligibleverticalmotion(FigureS4).Theresidualsare approximatelyGaussianwithmeanvaluesclosetozero.TheaverageRMSmisfitis5mm/yr,whichgivesan empiricaluncertaintyof∼4mm/yrfortheindividualtracks. Foreverypixelwhereinformationfrombothascendinganddescendinggeometriesareavailable,weuse equation(3)toinvertfortheeast-westandverticalcomponentsofmotionfollowingthemethoddescribed byWrightetal.[2004]andHussainetal.[2016]whiletakingintoaccountthelocalincidenceangles: ⎡D ⎤ E D =[sin(𝜃)cos(𝛼) −sin(𝜃)sin(𝛼) −cos(𝜃)]⎢D ⎥, (3) LOS ⎢ N⎥ ⎣D ⎦ U whereD istheLOSvelocity,𝜃isthelocalradarincidenceangle,𝛼istheazimuthofthesatelliteheading LOS vector,and[D ,D ,D ]Tisavectorwiththeeast,north,andverticalcomponentsofmotion,respectively. E N U HUSSAINETAL. INTERSEISMICCENTRALNAF 7 Journal of Geophysical Research: Solid Earth 10.1002/2016JB013108 Figure6.Descendingandascendingline-of-sightvelocitieswitheachtrackreferencedtoaEurasiafixedGNSSreference frame.Redcolorsindicatemotionawayfromthesatellite,whilebluecolorsindicatemotiontowardthesatellite. Equation (3) contains three unknowns (D ,D and D ), but we only have two input velocities with large E N U differences in satellite look angle in the inversion (the ascending and descending InSAR LOS velocities). Therefore,itisimpossibletocalculatethefull3-Dvelocityfieldwithoutapriorassumption.Thecommon assumptionmadeinpreviousstudiesisthatthereisnoverticalmotionacrosstheregionofinterest[e.g., Waltersetal.,2014;Hussainetal.,2016].Inourcasewenotethatboththeascendinganddescendingtracksare equallyinsensitivetomotioninthenorth-southdirection.Wethereforeusethesmoothinterpolatednorth componentoftheGNSSvelocities(FigureS5)toconstrainthenorth-southcomponent(D )intheinversion, N andsolvefortheeast-westandverticalcomponentsofmotionusingtheInSARLOSvelocities.Wecalculate thefault-parallelcomponentofthehorizontalvelocitybyassumingthatmotionoccursonastrike-slipfault ∘ trendingatN81 E. Ourfault-parallelvelocities(Figure7a)showtheexpectedright-lateralinterseismicmotionacrosstheNAF, with red colors representing motion to the northeast and blue to the southwest. Our estimated vertical componentshowsthatthereislittleverticalmotionacrosstheNAFinthisregion(Figure7b). Thereisarelativelysharpchangeinfault-parallelvelocitysouthoftheNAF(Figure7)thatcoincideswiththe B-B′profileline.Webelievethatthisisduetoacombinationofpostseismicdeformationfromthe2000Orta HUSSAINETAL. INTERSEISMICCENTRALNAF 8 Journal of Geophysical Research: Solid Earth 10.1002/2016JB013108 Figure7.LOSInSARvelocitiesdecomposedintothe(a)fault-paralleland(b)verticalcomponentsofmotion,wherethe north-southcomponentisconstrainedbytheGNSSnorthcomponent(FigureS5),seetextfordescription.Negative fault-parallelvelocitiesindicatemotiontowardthewest,andnegativefault-perpendicularvelocitiesindicatemotionto thesouth.UncertaintymapsforthesecomponentsareinFigureS6.ThelineslabeledA-A′,B-B′,andC-C′areprofiles throughthefault-parallelvelocityshowninFigure8.EarthquakemomenttensorsarefromtheGlobalCentroidMoment Tensorcatalogforalleventsgreaterthanmagnitude4between1976and2016.The2000M 6Ortaearthquake w locationisshowninFigure7a. earthquake(M 6)[Taymazetal.,2007],residualatmosphereintroducedmainlyfromascendingtrack71,and w postseismicdeformationfromthe1999IzmitandDüzceearthquakes. 5.ModelingProfileVelocities Weanalyzethreeprofilesacrossthefaultwherevelocitiesfromwithin20kmareprojectedontotheprofiles showninFigure7a.Waltersetal.[2014]notedthatthereisavariationinthefault-parallelvelocityawayfrom thefaultthatisnotduetointerseismicloadingbutduetotheproximitytotheEulerpoleofrotation.For example,GNSSvelocitiespresentedbyNocquet[2012]showfault-parallelvelocityvectorswithmagnitude ∼25mm/yrclosetotheNAFbut∼8mm/yrinCyprusroughly800kmawayfromthefault.Thisvariationis mostlyduetotheproximityoftheCyprusGNSSstationstothepoleofrotationofAnatoliawithrespectto Eurasia.WeusethepoleofrotationcalculatedforAnatoliawithrespecttoEurasiabyReilingeretal.[2006], ∘ ∘ ∘ whoestimatedarotationrateof1.23 /Myraboutapolelocatedat32.1 E,30.8 NneartheNiledelta.Ina Eurasia-fixedreferenceframethisrotationeffectonlyappliestotheregionsouthoftheNAFandcorresponds toavalueof𝜃 =0.0215mm/yr/kmor2.15mm/yratadistanceof100kmfromthefault. rot Assumingthefault-parallelvelocitiesfartosouthofthefault(>200km)aremostlyduetoatmosphericnoise and contain no tectonic deformation, we calculate the variance-covariance matrix of the noise using the ∘ ∘ methoddescribedinsection2,usingvelocitiesfroma50kmby50kmregioncenteredon32.5 E,39 N.The estimatedvariance(𝜎2)andcharacteristiclength(𝜆)forthecovariancefunction(equation(1))is6.35(mm/yr)2 and35.8km,respectively. ProfilesA-A′ andC-C′ donotcrossthecreepingsectionofthefault.Fortheseprofileswefita1-Dmodel [Savage and Burford, 1973] through the profiles where the fault-parallel velocity, v , at a fault normal par HUSSAINETAL. INTERSEISMICCENTRALNAF 9 Journal of Geophysical Research: Solid Earth 10.1002/2016JB013108 Figure8.Profilesthroughthefault-parallelvelocitiesalongthreelinesshowninFigure7.Theredpointsare fault-parallelvelocitiesprojectedfromwithin±25kmdistanceontotheprofile.Thebluepointsarethefault-parallel componentoftheGNSSvelocities.Theboldblackdashedlineisthebestfit,maximumaposterioriprobability(MAP) solution,whilethelightgrey-shadedregionisthe95%modelconfidencerange.Thebestfitmodelparametersare showninthetextwiththe95%confidencerangeinbrackets. distancex,isafunctionofthefaultsliprate,S,andthelockingdepth,d .Includingtherotationeffectdiscussed 1 above,our1-Dmodelis ( ) { S x 0.0215, ifx>0 v (x)= arctan +x𝜃 +a, where𝜃 = , (4) par 𝜋 d rot rot 0, ifx≤0 1 whereaisastaticoffset. However,profileB-B′ crossesthecreepingsectionofthefault.Forthisprofilewemodelthefault-parallel velocityasacombinationoftwosignals:alongwavelengthsignalthatrepresentsinterseismicloadingatrate Sandlockingdepthd ,andashortwavelengthsignalthatrepresentsthefaultcreepatarateCfromthe 1 surfacedowntodepthd [e.g.,Wrightetal.,2001a;Elliottetal.,2008;Hussainetal.,2016]. 2 ( ) [ ( ) ] { S x 1 x 0.0215, if x>0 v (x)= arctan +C arctan −(x) +x𝜃 +a, where𝜃 = , (5) par 𝜋 d 𝜋 d rot rot 0, if x≤0 1 2 where(x)istheHeavisidefunction. Wefindbestfitvaluesforeachmodelparameter(S,d ,C,andd )andanoffseta,usingaBayesianapproach, 1 2 implementingtheGoodmanandWeare[2010]affine-invariantensembleMarkovchainMonteCarlo(MCMC) samplerwhileaccountingforthecovariance.FordetailsseeHussainetal.[2016]. OurMCMCsampleruses600walkerstoexploretheparameterspaceconstrainedby0 < S(mm/yr)< 60, 0 < d (km),< 60,0 < C(mm/yr),< 30,0 < d (km),< 40,−40 < a(mm/yr)< 40,assumingauniform 1 2 priorprobabilitydistributionovereachrange.Animportantconstraintweimposeisthatthemaximumcreep depthcannotbegreaterthanthelockingdepth,i.e.,d ≤d .OurMCMCmodelrunsover300,000iterations 2 1 andproduces48,000randomsamplesfromwhichweestimateboththemaximumaposterioriprobability (MAP)solutionandcorrespondingparameteruncertainties. TheresultsofouranalysisareshowninFigure8,withtheobservedprofilevelocityinredandtheMAPsolution inbolddashedline.Thesampledmarginalprobabilitydistributionsforthefaultsliprate,thelockingdepth, HUSSAINETAL. INTERSEISMICCENTRALNAF 10