Interseismic strain accumulation across the central North Anatolian Fault from iteratively ...

20 Pages·2016·5.96 MB·English

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