Constraints from GPS measurements on the dynamics of deformation in Anatolia and the Aegean

29 Pages·2017·6.14 MB·English

Journal of Geophysical Research: Solid Earth
RESEARCHARTICLE Constraints from GPS measurements on the dynamics
10.1002/2016JB013382 of deformation in Anatolia and the Aegean
SpecialSection:
Philip England1,Gregory Houseman2,andJean-Mathieu Nocquet3
StressatActivePlate
Boundaries—Measurement
andAnalysis,andImplications 1DepartmentofEarthSciences,UniversityofOxford,Oxford,UK,2InstituteofGeophysicsandTectonics,UniversityofLeeds,
forSeismicHazard Leeds,UK,3ObservatoiredelaCôted’Azur,UniversitéCôted’Azur,IRD,CNRS,Valbonne,France
KeyPoints: Abstract WeestimatethestrengthofthelithosphereinAnatoliaandtheAegean,andtheboundary
•Physicalmodelthattreatslithosphere
ofAnatoliaandtheAegeanasafluid forcesactinguponit,usingadynamicalmodelthattreatsthelithosphereasathinfluidsheetdeforming
sheetfitsGPSvelocitieswithRMS inresponsetovariationsingravitationalpotentialenergy.Thismodelhasonefreematerialparameter,the
misfitof5mm/yr
powerlawexponent,n,oftheverticallyaveragedrheologyofthelithosphere,andtwoparametersthat
•Distributionsofstrain,earthquakes,
andactivefaultingdeterminedby specifytheforcesperunitlengthappliedtoitsedges.Solutionstothismodelthatbestfitthevelocities
lateralvariationsingravitational of346reliableGPSsitesrequireaneffectiveviscosityofthelithosphereof1022to1021Pasatstrainrates
potentialenergyoflithosphere
of10to100nanostrainperyear.Thebest-fittingforceattheArabia-Anatoliaboundaryisconsistentwith
•Variationinstrengthoflithosphere
notrequiredtoexplainlarge-scale thelithostaticpressureduetothehightopographythere,andtheforceattheNubia-Aegeanboundary
deformationthoughmoderate isconsistentwiththecontrastinlithostaticpressureacrossthatboundary.Noadditionalforce,
weakeninglikelynearNAnatoliafault
from“slabrollback”orbasaltractionsduetoconvectioninthemantle,isrequiredtoexplainthe
observations.Theseresultsaresupportedbyscalingrelationsderivedfromapproximateanalyticalsolutions.
SupportingInformation: Theinverserelationshipbetweentheviscosityofthelithosphereanddeviatoricstressproducesstrong
•SupportingInformationS1 slowlydeformingregionsintheSouthernAegeanandCentralAnatoliawhosemotionsresemblethose
•TableS1
ofmicroplates.Thedistributionofgeodeticstrainrateswithintheregion,andthepartitioningbetween
normalandstrike-slipfaulting,areexplainedbytheinterplaybetweenboundaryconditions,internal
Correspondenceto:
variationsingravitationalpotentialenergy,andthepowerlawrheologyofthelithosphere.
P.England,
[email protected]
1.Introduction
Citation:
England,P.,G.Houseman,and WeinvestigatetheforcesresponsibleforthedeformationofAnatoliaandtheAegean,thecontinentalregion
J.-M.Nocquet(2016),Constraints thatliesbetweeneasternTurkeyandthecoastsoftheEasternMediterraneanandIonianSeas(Figure1).
fromGPSmeasurementsonthe
TheregionisboundedtothesouthbytheNubianandArabianplatesandtothenorthbystableEurasia.
dynamicsofdeformationinAna-
toliaandtheAegean,J.Geophys. Hereinafter,weusethetermEurasiatorefertothepartoftheEurasiancontinentanditssurroundingocean
Res.SolidEarth,121,8888–8916, floorthatexhibitsnegligibleinternaldeformationand,unlessotherwisespecified,velocitiesarerelativeto
doi:10.1002/2016JB013382.
Eurasia.NubiaandArabiamoveroughlynorthwardwhereasAnatoliamovesnorthwest-to-westward,and
westernTurkeyandGreecemovewest-to-southwestward(Figure2a)[Aktug˘ etal.,2009;Floydetal.,2010;
Received25JUL2016
McCluskyetal.,2000;Reilingeretal.,2006].KinematicdescriptionsofthemotionofAnatoliaandtheAegean
Accepted9NOV2016
havelabeleditas“extrusion”or“escape”[e.g.,Armijoetal.,1999,2003;S¸engöretal.,1985]orhaveascribedit
Acceptedarticleonline15NOV2016
Publishedonline3DEC2016 tothe“rollback”oftheNubianslabbeneaththesouthernAegean[e.g.,LePichon,1982;LePichonandKreemer,
2010;Royden,1993].Dynamicalexplanationsforthismotionattributeittodifferencesingravitationalpoten-
tialenergybetweeneasternTurkeyandthesouthernAegean(Figure3),ortotractionsappliedtothebaseof
thelithospherebyflowofthemantle[e.g.,Faccennaetal.,2014;LePichon,1982;LePichonandAngelier,1979;
LePichonandKreemer,2010;McKenzie,1972,1978].
Althoughthelarge-scalemotionsofAnatoliaandtheAegeanrelativetoEurasiacanbedescribedapproxi-
matelybyarigid-bodyrotationwithapoleneartheNiledelta[LePichonandKreemer,2010;Reilingeretal.,
2006],activefaultspervadetheregionandGPSvelocitiesincreasesystematicallyby20mm/yrfromeast
towest[e.g.,Aktug˘ etal.,2009,2013;Floydetal.,2010;GoldsworthyandJackson,2000,2001;Goldsworthy
etal.,2002;Jackson,1994;JacksonandMcKenzie,1988;Reilingeretal.,2006;S¸arogluetal.,1992;Shawand
Jackson,2010;Taymazetal.,1991a](seeFigures1,2,and4b).Thisdeformationresultsinwidespreadlarge
shallowearthquakes(Figure1b)whosefocalmechanismsvaryfromthrustandreversefaultingalongthe
westernandsouthernmarginsoftheregiontonormalfaultingthroughoutmainlandGreeceandwestern
Turkey,tonormal-plus-strike-slipfaultinginthenorthernAegean,andtostrike-slipfaultingonthenorthand
©2016.AmericanGeophysicalUnion.
AllRightsReserved. southeasternmarginsofAnatolia.
ENGLANDETAL. DYNAMICSOFANATOLIAANDTHEAEGEAN 8888
Journal of Geophysical Research: Solid Earth
10.1002/2016JB013382
Figure1.Tectonicsetting.(a)Locationmapshowingplacesreferredtointhetext.KFZ:KefaloniaFaultZone;PST:PlinyandStraboTrenches;NAF:North
AnatolianFault;EAF:EastAnatolianFault.Activefaultsoftheregionareshowninblue(strikeslip)orred(normal)[GoldsworthyandJackson,2000,2001;
Goldsworthyetal.,2002;S¸arogluetal.,1992].YellowtrianglesindicatethelocationsoftheactivevolcanoesoftheAegeanarc.(b)Seismicactivityoftheregion.
Dotsshowtheepicentersofearthquakesbetween1960and2007withhypocentersshallowerthan33km[InternationalSeismologicalCentre,2011].Focal
mechanismsfromtheGCMTcatalog[Dziewonskietal.,1981;Ekströmetal.,2012]andfromJacksonandMcKenzie[1988];ShawandJackson[2010].
ENGLANDETAL. DYNAMICSOFANATOLIAANDTHEAEGEAN 8889
Journal of Geophysical Research: Solid Earth
10.1002/2016JB013382
Figure2.Velocities,relativetoEurasia,ofthe346GPSsitesintheAegeanandAnatoliausedasconstraintsonthephysicalmodel(section2).Doublelineshows
regioninwhichthecalculationsofthispaperarecarriedout(Figure3).
Thevelocitydifferencesof∼20mm/yrwithinAnatoliaandtheAegeanhavebeendescribedbykinematic
modelsthatcontainbetweentwoandsixelasticblocks,withpermanentdeformationbeingconcentrated
aroundtheiredges[e.g.,LePichonetal.,1995;McCluskyetal.,2000;NystandThatcher,2004;Reilingeretal.,
2006].Alternatively,theyhavebeeninterpretedasthestrainofaductilelithospherewithinwhichstronger
enclavesareembedded[e.g.,Billirisetal.,1991;Clarkeetal.,1998;Daviesetal.,1997;Hollensteinetal.,2008]
that deforms in response to boundary tractions and internal gradients in gravitational potential energy
[e.g.,Hatzfeldetal.,1997;LePichon,1982;LePichonandAngelier,1979;LePichonandKreemer,2010;Martinod
etal.,2000;McKenzie,1972,1978;ÖzerenandHolt,2010].
ThedifferencesamongthesedescriptionsofthemotionsofAnatoliaandtheAegean,andoftheirinter-
naldeformation,arerootedinquestionsabouttherelativemagnitudesofforcesactingonthelithosphere,
andabouttheimportanceoflateralvariationsinitsstrength.Weinvestigatetheseissuesusingaphysical
modelconstrainedbyaregional-scaleGPSvelocityfield.First,weshowthatthestressesassociatedwithlat-
eraldifferencesingravitationalpotentialenergyarelargeenoughtocausetheobserveddeformationofthe
lithosphere.Second,weshowthatthedistributionofthisdeformationdoesnotrequirestronglateralvaria-
tionsintheintrinsicphysicalpropertiesofthelithosphere.Furthermore,weshowthattheboundarytractions
requiredtoexplaintheobservedfieldofdeformationareconsistentwiththegravitationalpotentialenergy
(GPE)contrastsacrosstheAegeanandAnatolianboundariesandthattractionsactingonthebaseofthe
lithospherearenotrequiredtoexplainthedeformation.
2.GPSObservations
OuranalysisisconstrainedbymeasurementsofcrustalvelocitiesmadeusingtheGlobalPositioningSystem
(GPS). Nocquet [2012] combined 24 GPS studies in a common reference frame, checking for consistency
amongsolutionsandreevaluatinguncertaintiesinvelocity.WeaddedmeasurementsbyAktug˘etal.[2013,
2015],Ozeneretal.[2010];Tataretal.[2012],andYava¸sog˘luetal.[2010]asindividualsolutionsandmadea
fullcombinationusingtheproceduresofNocquet[2012].Wethenselected,usinginformationfromtheorig-
inalpublications,onlythosestationsthatwereoccupiedonatleastthreeoccasions,overatimespanof
longerthan4years.Allsiteshavingaformalerrorgreaterthan2mm/yrintheirvelocitywereeliminated,
leaving346measurementsofhorizontalvelocitywithintheregionshowninFigure2.Theseareconsiderably
fewerthanthe872velocitiesusedbyÖzerenandHolt[2010],butthatdatasetcontainedmanysitesreported
bytwoormorepublications,andmanysitesthatwehaveexcludedbythecriteriaofnumberanddura-
tionofoccupations.Wethereforeexpectthatthedatasetweuse,inadditiontohavingmoreinformationin
Anatolia[Aktug˘etal.,2013,2015;Ozeneretal.,2010;Tataretal.,2012;Yava¸sog˘luetal.,2010],isappreciablyless
ENGLANDETAL. DYNAMICSOFANATOLIAANDTHEAEGEAN 8890
Journal of Geophysical Research: Solid Earth
10.1002/2016JB013382
Figure3.TopographyandGPEwithinthecalculationdomain(delineatedbylabeledboundarysegments).(a)Surfaceheightwithinthecalculationdomain
smoothedwithaGaussianfilterofwidth100km.(b)GPEcalculatedfromsurfaceheight(section3.1,equation(5),andparametersgiveninTable1).Lineslabeled
AtoFindicatethelocationsatwhichdifferentboundaryconditionsareappliedinthecalculations.ConditionsofzerovelocityareappliedonsegmentsA,B,and
D.ForcesperunitlengthΓ andΓ areappliedtosegmentsCandE,asdiscussedinthetext(section3.2).Zerotangentialtractionandzeronormalvelocityare
W E
appliedtothenorth-southboundarysegmentF.
ENGLANDETAL. DYNAMICSOFANATOLIAANDTHEAEGEAN 8891
Journal of Geophysical Research: Solid Earth
10.1002/2016JB013382
Figure4.Illustrationofanalyticalsolutionforthevelocityinarectangularthinviscoussheet(AppendixA).(a)The
geographicalregionofinterest(boundedbygraylines)isapproximatedbyaboxofeast-westlengthlandnorth-south
dimension2w(doublelines);theslightcurvatureisneglected.DotsshowlocationsofGPSsiteswhosevelocities
(magnitudes)areplottedinFigure4b);siteswithin150kmoftheprofileareused,exceptthoseontheEurasianor
Arabianplates.(b)ProjectionofmagnitudesofGPSvelocitiesalongthelinex=0.Gray-shadedregionshowstherange
ofvelocitiescalculatedforthebest-fittinganalyticalsolution(equation(A6),andseesection4.)Upperboundaryofgray
shadecorrespondstothesolutionaty=0,lowerboundtothesolutionaty=±150km.(c)Comparisonbetween
analyticalsolutionandobservations(intheanalyticalsolutionpositiveu iseastward).
x
ENGLANDETAL. DYNAMICSOFANATOLIAANDTHEAEGEAN 8892
Journal of Geophysical Research: Solid Earth
10.1002/2016JB013382
influencedbymeasurementerror.Themeanspacingbetweensitesisabout50km,thoughsiteswithinthe
AegeanSeaareseparatedbyupto100km.
3.PhysicalModel
3.1.MathematicalFormulation
Wetreatthecontinentallithosphereasasheetoffluidwhosedeformationisdeterminedbythevertical
averagesofitsphysicalpropertiesandoftheforcesactinguponit—thethinsheetapproximation[Birdand
Piper,1980;EnglandandMcKenzie,1983;HousemanandEngland,1986](seeSchmalholzetal.[2014],fora
recent discussion that corrects some misapprehensions that have arisen in the literature concerning this
approximation).Geodeticmeasurementsofthedeformationofthelithospherearetakenatthetopofthe
uppercrust,whichingeneraldeformselastically.Wetaketheviewthat,exceptduringthepostseismicinter-
valsofnearbyearthquakes,thestrainofthisthinelasticlayerfollowscloselytheductilestrainofthelayer
beneath[e.g.,BirdandPiper,1980;Bourneetal.,1998;HousemanandEngland,1986].Thisviewcontrastswith
interpretationsofvelocityfieldsthattreatthelithosphereasanelasticmediumwhosesurfacedeformation
iscontrolledbytheslipofdislocationsburiedwithinit[e.g.,Hammondetal.,2011;Meadeetal.,2002;Meade
andHager,2005;Reilingeretal.,2006;Wallaceetal.,2007].
GarthwaiteandHouseman[2011]compared3-Dcalculationswiththeresultsofthe2-Dthin-sheetmodeland
showedthatintheabsenceofstrongverticalstratificationofphysicalproperties,differencesareminor,pro-
videdthatthehorizontalscaleofthedeformingzoneandthelengthscaleoverwhichsignificantchanges
inboundaryconditionoccur,aregreaterthanthethicknessofthelithosphere.Thefirstoftheseconditions
ismet:thehorizontalscalesoftheregionofinterestareoftheorderofathousandkilometerswhereasthe
lithosphericthicknessisoftheorderofahundredkilometers.Thesecondconditionmaybeviolatedinsmall
regions,suchaseitherendoftheHellenicplateboundary,wherethereareabruptchangesinboundary
condition.
Thedepth-averagedrheologyofthethinsheetisdescribedbyapowerlawrelationbetweendeviatoricstress
andstrainrate:
( )
𝜏 =BĖ n1−1 𝜀̇ (1)
ij ij
where𝜀̇ istheijthcomponentofthestrainratetensor,assumedconstantwithdepththroughthesheet,𝜏
ij ij
istheijthcomponentofthedeviatoricstress,withtheoverbardenotingaveragesoverthethickness,L,ofthe
lithosphere,nisthepowerlawexponentfortherheology,and
( )
Ė = 𝜀̇ 𝜀̇ 12 (2)
kl kl
is the second invariant of the strain rate tensor and the convention of summation over repeated indices
applies.Thefluidisassumedtobeincompressible(𝜀̇ =0).
kk
SonderandEngland[1986]showedthatequation(1)closelyapproximatestheverticallyaveragedrheologyof
alithospheredeformingbyacombinationofcreepmechanisms,whicharelargelytemperature-dependent,
andfrictionalmechanisms,whichdependprincipallyonpressure.Furthermore,ifthebrittleuppercrustis
cutbynumerousfaults,thenitmaybeappropriatetodescribeitslong-termdeformationfieldbyapower
lawrheologyoftheformofequation(1)[e.g.,NanjoandTurcotte,2005].Thesingleparameter,n,allowsthe
depth-averagedrheologyforthelithospheretorangefromNewtonianbehavior(whenn = 1)toalmost
plasticbehaviorforlargen.
Thequantity
( )
𝜂 = 1BĖ n1−1 = 1BnT(1−n), (3)
eff 2 2
( )
1
whereT isthesecondinvariantoftheverticallyaverageddeviatoricstresstensor, 𝜏 𝜏 2,iscommonly
kl kl
referredtoastheeffectiveviscosity,byanalogywiththeviscosityofaNewtonianfluidinthecasethatn=1.
Notethatforn ≠ 1theeffectiveviscositydependsonthedeviatoricstress,orequivalentlyuponstrainrate,
sothatlithosphereofspatiallyinvariantviscositycoefficient,B,willingeneralexhibitlateralvariationinits
effectiveviscosity.
ENGLANDETAL. DYNAMICSOFANATOLIAANDTHEAEGEAN 8893
Journal of Geophysical Research: Solid Earth
10.1002/2016JB013382
Table1.Notation
Parameter Definition Value/Reference
B Viscositycoefficient equation(1)
g Accelerationduetogravity 9.8ms−2
h Surfaceheight Figures1and3
h Referencesurfaceheight 100m
0
L Thicknessofthelithosphere 100km
S Crustalthickness
S Referencecrustalthickness 35km
0
u ithcomponentofvelocity
i
U velocityscale equation(8)
0
𝜀̇ ijthcomponentofstrainrate
ij
̇
E Secondinvariantofstrainratetensor equation(2)
𝜂 Effectiveviscosity equation(3)
eff
Γ Gravitationalpotentialenergy Nm−1
perunitsurfacearea orJm−2
𝜌 Densityofcrust 2800kgm−3
c
𝜌 Densityofmantle 3200kgm−3
m
𝜏 ijthcomponentofdeviatoricstress
ij
Theequationsgoverningthedeformationofthethinsheet(equation(A14))expressthebalancebetween
horizontalgradientsofthedeviatoricstressrequiredtodeformthesheetandgradientsinitsgravitational
potentialenergyperunitarea(GPE).WeestimatetheGPEoflithosphericcolumnsusingthedistributionof
surfaceheightshowninFigure3a,havingappliedaGaussianfilterofwidth100kmtothesurfaceheight
withinthecalculationdomain.Wedefineareferencecolumnoflithosphere,withcrustalthicknessS ,surface
0
heighth ,anddensitiesasstatedinTable1.Surfaceheightabovesealevel,h,isrelatedtocrustalthickness,S,
0
undertheassumptionofAiryisostaticbalance
𝜌 ( )
S−S = m h−h , (4)
0 𝜌 −𝜌 0
m c
where𝜌 and𝜌 aretheaveragedensitiesofcrustandmantle.Foruniformcrustoverlyinguniformmantle,
c m
thedifferenceinGPEbetweenthereferencecolumn,ofGPEΓ ,andacolumnwithcrustofthicknessSis
0
( ( ))
( ) 𝜌 h−h
Γ−Γ =g𝜌 h−h S + m( 0) (5)
0 c 0 0 2 𝜌 −𝜌
m c
wheregistheaccelerationduetogravity;forcolumnswithh<0theseexpressionsareadjustedforthe
contributionfromthewatercolumn[HaxbyandTurcotte,1978].
WithinthecontinentallithosphereoftheregiontheGPE(Γ−Γ )variesfromahighof2TNm−1intheeasttoa
0
lowof−1TNm−1inthesouthernmostAegean(Figure3b).TheassumptionofAiryisostasyisnotcritical;ifiso-
staticcompensationismaintainedinpartbydensityvariationsinthemantle,thentheGPEforagivensurface
heightwouldbeslightlyhigherorlowerthanthatgivenbyequation(5),resultinginminorlocaladjustments
tothevelocityfield.Thesmallsizeofisostaticgravityanomaliesintheregion[Balminoetal.,2012]rulesout
significantnonisostaticcontributiontoGPEfromliftingofthelithospherebyconvectivetractionsappliedto
itsbase,sometimesreferredtoas“dynamictopography”[e.g.,Molnaretal.,2015,Figure5].
3.2.CalculationDomainandBoundaryConditions
Wedividethebordertothecalculationdomainintosixsectionswhoselengthsandlocationsaredictatedby
theirgeologicalsetting(Figure3b).Theadoptionofasmallnumberoflongboundarysegmentsisjustifiedby
afundamentalaspectofthegoverningequations:aconditionofstressorvelocity,whenap√pliedoveralength
wofboundary,influencesstrainratesoveradistancefromthatboundarythatis∼w∕(𝛼 n),wherenisthe
exponentinthepowerlawrheology(equation(1))and𝛼dependsonthenatureoftheboundarycondition
[Englandetal.,1985;Sonderetal.,1986].Whentheboundaryconditionisoneofnormalvelocityortraction,
ENGLANDETAL. DYNAMICSOFANATOLIAANDTHEAEGEAN 8894
Journal of Geophysical Research: Solid Earth
10.1002/2016JB013382
𝛼∼1;whentheconditionistangential,thelengthscaleis4timesshorter(𝛼∼4).Itfollowsthatanyvariation
inboundaryconditionoverascalethatissmallcomparedwiththesizeoftheregionofinteresthasonlyalocal
influenceanddoesnotsignificantlyaffectthedistributionofstressesandvelocitiesfarfromthatboundary.
Thisassertionissupportedbyinvestigationofalternative,morecomplex,boundaryconditionsinsection5.4.
ThereferenceframeisfixedtoEurasiaandwesetvelocitytozeroalongaboundary(segmentAinFigure3b)
thatfollowsthesouthernedgeoftheBlackSea,whichwepresume,fromthelackofdeformationwithinandto
thenorthofit,tobepartofEurasia.Atitswesternend,segmentArunstothenorthoftheregionsofextension
andconvergenceinGreece.Althoughthereissomeactivedeformationtothenorthofthisline,velocities
therearegenerallylowerthan2–3mm/yrrelativetoEurasia[Burchfieletal.,2006;Kotzevetal.,2006].
WeplacetheboundarybetweentheAegeanlithosphereandtheNubianplatealongalinerunningtothe
landwardsideoftheHellenic,Pliny,andStraboTrenches,whichrepresenttheoutcropsofmajorfaultsatthe
southernborderoftheAegean[e.g.,LePichonandAngelier,1979;MascleandChaumillon,1998].Segment
CinFigure3bcorrespondsapproximatelytothedowndiplimitofshallowlydippingthrust-faultingearth-
quakesontheplateinterfaceandrepresentstheboundarybetweentheaccretionaryprismandtheoverriding
Aegeancontinentalcrust[e.g.,ShawandJackson,2010].Itisprobablethat80%to90%oftheconvergence
acrossthisboundaryoccursbysteadyaseismicsliponthesubductioninterface[JacksonandMcKenzie,1988;
Reilingeretal.,2009;Vernantetal.,2014](thoughseeGanasandParsons[2009]foranalternativeviewpoint);
itisthereforedifficulttoconstrainthevelocityconditionforthisboundary.Instead,weapplyanormalforce
perunitlength,Γ ,andallowthistobeafreeparameterinthecalculations.
W
Similarly,itisappropriatetospecifyanormalforceperunitlengthattheboundarybetweenArabiaand
Anatolia(segmentEinFigure3b);againweallowthisforce,Γ ,tobeafreeparameterinthecalculations.We
E
∘
continuethecalculationdomainto41 E(segmentFinFigure3b)forsimplicitysettingtheeast-westvelocity
andthetangentialtractiononthisboundarytobezero.Thissegmentisshortenoughthatitslocationand
conditionsdonotappreciablyinfluencevelocitiesthroughoutmostofthecalculationdomain.
TherelativevelocitiesacrosstheboundariesbetweenTurkeyandNubia(segmentDinFigure3b)andbetween
theIonianSeaandthenorthwesternpartoftheregion(segmentB)arepoorlyconstrained.SegmentDlies
alongadiffusezoneofdeformationbetweenRhodesandtheGulfofIskenderun[e.g.,Aksuetal.,2009;Hall
etal.,2009;tenVeenetal.,2004;Woodsideetal.,2002],anditisunclearwhatfractionoftherelativemotion
betweenNubiaandAnatoliaisabsorbedsouthofthecoastofTurkey.Equally,therateofconvergenceacross
theboundarywiththeIonianSeaisafewmillimetersayear[e.g.,D’Agostinoetal.,2008]anditisunclearwhat
fractionofthatrateisabsorbedoffshore.Inordertomaintainthesimplicityofourcalculationsweinitiallyset
thevelocityoneachsegmenttozero,investigatingmorecomplexboundaryconditionslater(section5.4).
4.PreliminaryAnalysis
Wefirstinvestigatetheobserveddistributionofvelocitiesusingananalyticalsolutionfortheflowofathin
sheetofNewtonian fluid(equation(A6)insection A1). Inthisapproximation,the geometricallycomplex
boundaryconditionsonAnatoliaandtheAegeanarereplacedbyaboxboundedbyplanesatx=0,l,andat
y = ±w(Figure4a).Theplanesaty = ±wrepresentthenorthernandsouthernboundariesofAnatoliaand
theAegean.Theplanex =0approximatestheHellenicplateboundary,andx =lapproximatestheeastern
boundarytoAnatolia.TheGPEinthisboxvarieslinearlyfromΓ atx =0toΓ atx =l,givingagradientof
W E
GPEthatis
ΓE−ΓW = ΔΓ. (6)
l l
Theplaney=wistakentohavezerovelocity,i.e.,tobefixedtoEurasia.Theplaney=−wrunsclosetoseg-
mentsE,D,andpartofC(Figure3b).Asdiscussedinsection3.2,thesetofconditionsalongthisboundary
cannotbecapturedbyasingleparameter;tomaintainthesimplicityoftheapproximation,thevelocityon
thisboundaryisalsosettozero.Theboundaryconditionsonx=0,lareexpressedasforcesperunitlength,
ΔΓ ,ΔΓ ,whichrepresentcontrastsinGPEacrosstheboundaries(positivevaluesofthesequantitiescor-
W E
respondtolowerGPEoutsidetheboundary,therefore,toanetextensionalforceactingonthatboundary:
equations(A4)and(A5)).
ENGLANDETAL. DYNAMICSOFANATOLIAANDTHEAEGEAN 8895
Journal of Geophysical Research: Solid Earth
10.1002/2016JB013382
∘ ∘
The east-west dimension, l, of the geographical region of interest is about 1800 km (from 20 E to 41 E,
Figure4a).Thenorth-southdimensionvariesfromabout400kmintheeastto800kminthewest(Figure4a);
weadoptavalueof600kmforthewidthofthebox(w=300km).Wedonotexpecttheapproximatesolution
toholdnearthecoastlinesofTurkey,wheretheconfigurationoftheboundariesiscomplex.InFigure4bwe
comparethissolutionwiththevariationinGPSvelocitiesalongaswathfromeasternAnatoliatotheHellenic
plateboundarythatis300kmwide.
Theanalyticalsolution(equation(A6))hasthreevelocityscales:ΔΓw2∕(2𝜂Ll)and(ΔΓ orΔΓ )×4w∕(𝜋2𝜂L),
W E
where𝜂istheviscosityofthelithosphereandLitsthickness.Wefound,bylinearleastsquaresinversion,the
valuesoftheseparametersforwhichtheanalyticalsolutionbestfitstheobservations.Thebest-fitvelocity
scaleassociatedwiththeGPEdifference,ΔΓ,is−29mm/yr.WithΔΓ∼3×1012Nm−1(Figure3b),l∼1800km,
w∼300km,andL∼100km,thisscaleyieldsanestimateof∼1021Pasfortheviscosityofthelithosphere.The
best-fittingvelocityscalesassociatedwiththeGPEdifferenceatx =0andl,ΔΓ ,andΔΓ are,respectively,
W E
7mm/yrand22mm/yr;with𝜂 ∼ 1021Pas,thesescalesgiveforcesperunitlengthofΔΓ ∼0.2TNm−1on
W
x=0andΔΓ ∼0.6TNm−1(bothinrelativetension)onx=l.
E
The approximate solution (Figures 4b and 4c) displays the principal features of the variation in velocity
betweeneasternAnatolia,despitetheconsiderabledifferencesbetweentheidealizeddistributionofGPE
(equation(6))andthatshowninFigure3b.Thisanalysisprovidesevidencethatthelarge-scaledeformation
oftheregionisdescribedbythebalancebetweengradientsofGPEanddeviatoricstressinthelithosphere.
5.NumericalExperiments
5.1.MeansofSolution
We solve the governing equations (equation (A14)) subject to the boundary conditions described in
section3.2,andwiththedistributionofGPEderivedfromtopography(section3.1,Figure3b),usingafinite
elementapproximationtotheequationsfordeformationonthesurfaceofasphere(seesectionA2).Thefinite
elementmeshcontains5400nodalpoints,346ofwhichareplacedatthelocationsoftheGPSsites.Themean
spacingbetweennodalpointsisapproximately15km,smallerthanthetypicalspacingsbetweenGPSsites
(Figure2b).Wecheckedtheresolutionofthecalculationsusingmeshesof13,000pointsandfoundthatthe
RMSmisfittotheobservationsdifferedbylessthan1%betweenthetwomeshes.
AsdiscussedinsectionA2,theequationsarereducedtonondimensionalformsthatdependonasingle
parameter,theArgandnumber
g𝜌 L(1−𝜌 ∕𝜌 )
Ar= c( )c m , (7)
1∕n
B U ∕L
0
whose numerator represents a stress scale associated with density variations within the lithosphere and
whosedenominatorrepresentsthestressrequiredtodeformthelithosphereatareferencestrainrateU ∕L
0
(seeequation(1))[EnglandandMcKenzie,1982].
IneachcalculationwechoseavalueforArwhichensuredthatthemaximumdimensionlessvelocitywasof
theorder1anddeterminedU byminimizingthemisfitfunction
0
[ ]
1
1 ∑N | |2 2
M= |uGPS−U u′| (8)
N | i 0 i|
i=1
whereuGPSisthevelocityoftheithGPSsite,u′isthedimensionlessvelocityofthesamesiteinthecalculation,
i i
andNisthenumberofGPSsites.ThisvalueofU thenimpliesavalueoftheviscositycoefficient
0
g𝜌 L(1−𝜌 ∕𝜌 )
B= c ( )c m . (9)
1∕n
Ar U ∕L
0
Foranygivencalculation,thevalueofBisindependentofthechosenvalueofAr(seesectionA2).Withthe
configurationofboundaryconditionsfixed,andwiththedistributionofGPEspecifiedfromthetopogra-
phy(Figure3b),therearethereforeonlythreefreeparameters:theexponentninthepowerlawrheology
(equation(1))andthemagnitudesofΓ andΓ ,theforcesperunitlengthappliedatthewesternandeastern
W E
boundaries(Figure3b).
ENGLANDETAL. DYNAMICSOFANATOLIAANDTHEAEGEAN 8896
Journal of Geophysical Research: Solid Earth
10.1002/2016JB013382
Table2.ValuesofParametersfortheBestFits(Figure5)totheObservedVelocities(Figure2),WithnandtheBoundary
ForcesΓ andΓ asFreeParametersa
W E
Γ Γ B M(equation(8))
W E
TNm−1 TNm−1 Pas1∕n mmyr−1
n=1 0.625 −1.5 2.4×1021 5.1
n=3 1.0 −1.75 6.0×1011 4.7
n=5 1.25 −2.0 8.1×109 4.7
n=9 1.625 −2.25 5.7×108 4.8
aForaNewtonianfluid,theviscosity,𝜂,isequaltohalftheviscositycoefficient,B(equation(3)).
5.2.EstimationofBoundaryForcesandLithosphericRheology
WecalculatedvelocityfieldsvaryingΓ and−Γ betweenzeroand10TNm−1instepsof0.125TNm−1,and
W E
withn = 1,3,5,and9.Solutionsforn = 5exhibitnofeaturethatcannotbeillustratedbythesolutionsfor
n=3and9sotheyarenotshown,butthebest-fittingparametersarelistedinTable2.
ThemisfitsdefineanelongatevalleyontheΓ −Γ plane(Figure5)withminimanearΓ =−(Γ +1TNm−1).
W E E W
ValuesofΓ yieldingthebestfitstoobservationslieintherange0.75to2TNm−1 andthoseforΓ liein
W E
Figure5.Root-mean-squaremisfitsbetweencalculatedandobservedvelocityfieldsforathinviscoussheetwith(a)n=1,(b)3,and(c)9subjectedtothe
boundaryconditionsandtheinternaldistributionofgravitationalpotentialenergyshowninFigure3b.ColorsandblackcontoursshowthevariationofRMS
misfits(equation(8))withΓ andΓ .Best-fittingcombinationsofΓ andΓ areshownbywhitedots(seeTable2).Doublelinesshowcontoursofthelogarithm
W E W E
tobase10oftheviscositycoefficient,B,forthesolutions(unitsofPas1∕n).(d)DistributionsofmisfitsforthecaseinwhichtherearenointernalgradientsofGPE.
ENGLANDETAL. DYNAMICSOFANATOLIAANDTHEAEGEAN 8897