Upper mantle structure beneath the Galápagos Archipelago from surface wave tomography

25 Pages·2007·2.16 MB·English

JOURNAL OF GEOPHYSICAL RESEARCH, VOL.112, B07303,doi:10.1029/2006JB004672, 2007
Click
Here
for
Full
Article
Upper mantle structure beneath the Gala´pagos Archipelago from
surface wave tomography
DarwinR.Villago´mez,1DouglasR.Toomey,1EmilieE.E.Hooft,1andSeanC.Solomon2
Received31July2006;revised2February2007;accepted7March2007;published4July2007.
[1] We present a Rayleigh wave tomographic study of the upper mantle beneath the
Gala´pagos Archipelago. We analyze waves in 12 separate frequency bands (8–50 mHz)
sensitive to shear wave velocity (V ) structure in the upper 150 km. Average phase
S
velocities are up to 2 and 8% lower than for 0- to 4-My-old and 4- to 20-My-old Pacific
seafloor, respectively. Laterally averaged V is 0.05–0.2 km/s lower between 75- and
S
150-km depth than for normal Pacific mantle of comparable age, corresponding to an
excess temperature of 30 to 150(cid:1)C and (cid:1)0.5% melt. A continuous low-velocity volume
that tilts in a northerly direction as it shoals extends from the bottom of our model to the
base of a high-velocity lid, which is located at depths varying from 40 to 70 km. We
interpret this low-velocity volume as an upwelling thermal plume that flattens against the
base of the high-velocity lid. The high-velocity lid is (cid:1)30 km thicker than estimated
lithospheric thickness beneath the southwestern archipelago, above the main region of
plume upwelling. We attribute the thicker-than-normal high-velocity lid to residuum
fromhotspot melting. The thicknessofthelid appearstocontrol thefinal depthof
melting andthevariability ofbasalt composition inthearchipelago. Atdepths lessthan
100–120km,plumematerialspreadsindirectionsbothtowardandagainsteastwardplate
motion, indicating that plumebuoyancyforcesdominate overplatedragforces and
suggesting arelatively highplume buoyancyflux(B(cid:2)2000kg/s).
Citation: Villago´mez,D.R.,D.R.Toomey,E.E.E.Hooft,andS.C.Solomon(2007),UppermantlestructurebeneaththeGala´pagos
Archipelago fromsurfacewavetomography,J.Geophys.Res., 112, B07303,doi:10.1029/2006JB004672.
1. Introduction be driven by pressure gradients associated with plate crea-
tion at ridges [Yale and Phipps Morgan, 1998; Toomey et
[2] Hot spot volcanism is widely thought to be the result al., 2002b].
of upwelling and melting of hot, buoyant mantle [Morgan,
[3] Although regional seismic tomography has provided
1971]. Gravity and topography observations of hot spot
compelling evidence for plume-like upwelling in the upper
swells and results from modeling suggest that these mantle
mantle[e.g.,Granetetal.,1995;Wolfeetal.,1997;Ritteret
upwellings or plumes rise to the base of the lithosphere
al., 2001; Allen et al., 2002; Li and Detrick, 2006], the
where they spread laterally [e.g., Ribe and Christensen,
resolution of images of off-axis hot spots in the uppermost
1994; Feighner and Richards, 1995; Sleep, 1996]. Some
mantle has not been adequate to provide clear tests of
plumes interact with nearby spreading centers and produce
models of plume spreading and plume-ridge interaction.
physical and chemical anomalies along some 15–20% of
To address these issues we present a surface wave tomo-
the global mid-ocean ridge system, although the precise
graphicstudyoftheGala´pagosArchipelago.Themaingoal
mechanismanddepthoftransportofplumematerialtomid-
ofourstudyistocharacterizetheupwelling,spreading,and
oceanridgesarestillmattersofdebate[Itoetal.,2003,and
melting of the shallow mantle beneath the Gala´pagos hot
references therein]. Lubrication theory and models predict
spot. The Gala´pagos Archipelago is an excellent setting to
thatforoff-axisplumesthemainfactorsthatcontrolplume-
study the dynamics of the interaction among hot spots, the
ridgeinteractionaregravitationalspreadingofaplumelayer
lithosphere, and a mid-ocean ridge, because of the proxim-
that pancakes beneath the sloping base of the lithosphere
ityofthehotspottotheGala´pagosSpreadingCenter(GSC)
and drag by the overriding plate [Ribe, 1996; Ito et al.,
and because the direction of plate drag over the hot spot
1997; Ribe and Delattre, 1998]. Alternatively, for a suffi-
(eastward) is approximately perpendicular to the relative
ciently low-viscosity asthenosphere plume flow could also
spreading direction of the GSC (north-south).
1Department of Geological Sciences, University of Oregon, Eugene, 2. Regional Setting
Oregon,USA.
2Department of Terrestrial Magnetism, Carnegie Institution of [4] TheGala´pagoshotspot,locatedintheequatorialeast
Washington, Washington, District of Columbia, USA. Pacific, consists of 10 major volcanic islands and 21
emergentvolcanoes(Figure1).Theislandssitonashallow
Copyright2007bytheAmericanGeophysicalUnion. and broad submarine platform (the Gala´pagos swell) that is
0148-0227/07/2006JB004672$09.00
B07303 1 of 25
B07303 VILLAGO´MEZ ETAL.: GALAPAGOS SURFACE WAVE TOMOGRAPHY B07303
Figure 1. (a) Map of the Gala´pagos Islands and seismic network. Triangles indicate seismic stations.
The black square and solid circle with a 100-km radius indicate the approximate center and area of a
region of anomalously thin mantle transition zone [Hooft et al., 2003]. The black arrow indicates the
direction of motion of the Nazca plate in a hotspot reference frame [Gripp and Gordon, 2002].
Bathymetry is from W. Chadwick (http://newport.pmel.noaa.gov/~chadwick/galapagos.html) 1000-m
contour interval. (b) Vector velocities of plates and other features in the Gala´pagos region.
elevated more than 2000 m above the surrounding ocean greater than 410 km [Hooft et al., 2003]. Above the
floor. An oceanic fracture zone crosses the northern part of downward-deflected 410-km discontinuity, body wave
theGala´pagosnear91(cid:1)Wandcreatesa(cid:1)5-Mylithospheric tomography resolves low seismic velocities at depths of
age offset, with thinner lithosphere beneath the eastern part 50–250 km, consistent with upwelling of anomalously hot
of the archipelago. The hot spot sits on the Nazca plate, mantle [Toomey et al., 2002a].
which moves eastward with respect to the hot spot, in a [6] Several uniquefeaturesmaketheGala´pagosArchipel-
direction approximately perpendicular to the north-south ago different from more conventional hot spots such as
(N-S) spreading of the GSC. Hawaii.First,whilevolcanism shows ageneral ageprogres-
[5] Evidence of a plume origin for the Gala´pagos hot sion in the direction of plate motion, the progression is not
spot includes basalts enriched in incompatible elements monotonic;almostalloftheGala´pagosIslandshaveeruptedin
(for example, higher 3He/4He and 87Sr/86Sr and lower theHolocene[SimkinandSiebert,1994].Second,geochem-
143Nd/144Nd) [White and Hofmann, 1978; Geist et al., icalsignaturesofbasaltsshowanunusualspatialdistribution:
1988; White et al., 1993; Kurz and Geist, 1999; Harpp depletedbasaltsappearnearthecenterandnortheasternpartof
and White, 2001], a general progression of the age of archipelago,while enrichedlavasappearprimarily alongthe
volcanism away from the hot spot in the direction of plate westernandsouthernparts[WhiteandHofmann,1978;Geist
motion[McBirneyandWilliams,1969;Sintonetal.,1996], etal.,1988;Whiteetal.,1993;KurzandGeist,1999;Harpp
and seismic images of anomalous upper mantle structure. andWhite,2001].Moreover,theGala´pagoshotspotinfluen-
Hooft et al. [2003] used receiver functions to show that the cesmagmatismandtectonicsontheGSC,aninferenceevident
410-km mantle discontinuity is deflected downward within in correlated variations in geophysical, geochemical, and
an area approximately 100 km in radius centered beneath volcanologicalmanifestationsalongtheridge[Detricketal.,
the southwestern corner of the archipelago (Figure 1). This 2002].ThesevariationsalongtheGSCaxisaremoreorless
anomaly reflects higher temperatures (130 ± 60 K) across symmetricabout91.5(cid:1)W[Schillingetal.,2003].
thatphasetransition,consistentwithupwellingfromdepths
2 of25
B07303 VILLAGO´MEZ ETAL.: GALAPAGOS SURFACE WAVE TOMOGRAPHY B07303
Figure2. Azimuthaldistributionofthe189events(M >5.9)forwhichRayleighwaveswereanalyzed
S
in this study. Epicentral distance varies from 40(cid:1) to 140(cid:1). Solid lines correspond to great circle paths.
Azimuthal equidistant projection centered at 0(cid:1)N, 90(cid:1)W.
[7] Two general models of mantle upwelling for the were deployed on nine islands of the archipelago between
Gala´pagos hot spot have been proposed. The first model September 1999 and March 2003 (Figure 1). The network
is based on geochemical data from the archipelago and consisted of 10 portable broadband stations and the Global
accounts for the Nazca plate moving eastward with respect Seismographic Network station PAYG. The station spacing
to the hot spot. Richards and Griffiths [1989] and White et was between 50 and 70 km. Three-component Streckeisen
al.[1993]suggestedthattheparticularspatialdistributionof STS-2sensorswereusedatallportablestations;twoGuralp
incompatibleelementsinthearchipelagocouldbetheresult CMG-3ESP instruments were initially deployed but were
ofthermalentrainmentofdepleteduppermantle,asaresult replaced after the first year. Data loggers were PASSCAL-
oflocalconvectiveoverturnwithinthecenterofadeflected equivalent Reftek units recording continuously at 20 sam-
mantle upwelling or plume. In this view, the plume is ples per second. The seismic network spanned an area
deflectedtotheeastintheshallowmantleinresponsetoplate approximately 200 km in diameter. For this study we use
drag. However, this model does not take into account the three-component recordings of Rayleigh waves generated
observed geochemical and geophysical variations along the by 186 teleseismic events with M > 5.9 at epicentral
S
axisoftheGSC.Thesecondmodelisbasedongeodynamical distances ranging between 40(cid:1) and 140(cid:1) (Figure 2).
modeling of hot spot-ridge interaction and accounts for the
effectofN-Sseafloorspreading[Itoetal.,1997].Themodel
3. Method
includes northward migration of the GSC relative to the hot
3.1. Imaging of Phase Velocity
spotbutdoesnotaccountforeastwardNazcaplatemotion.Hot
spot-derivedmaterialistransportedtotheridgesymmetrically [10] Theseismicdataarefirstusedtoderiveone-andtwo-
dimensional images of Rayleigh wave phase velocity. We
totheeastandwest,accountingforthesymmetricalalong-axis
obtain phase and amplitude information from the vertical-
geochemical variations, but the model does not consider the
component seismograms. After correcting for instrument
asymmetricalgeochemicalpatternsobservedwithinthearchi-
response, the data are windowed and filtered into 12 diffe-
pelago.Todate,nomodelcanaccountforallofthegeophys-
rent frequency bands using a 10-mHz-wide, fourth-order,
ical and geochemical observations at both the GSC and the
zero-phase Butterworth filter. The center frequencies are
Gala´pagosArchipelago.
between 8 and 50 mHz, or 20- to 125-s period (Figure 3),
[8] TostudytheuppermantlebeneaththeGala´pagosand
test models of plume-ridge-lithosphere interaction we con- corresponding to seismic wavelengths of (cid:1)80–500 km.
Fundamental Rayleigh waves are sensitive to structure to
ducted a broadband seismic experiment. Seismic stations
a depth of approximately one wavelength, with peak sensi-
3 of25
B07303 VILLAGO´MEZ ETAL.: GALAPAGOS SURFACE WAVE TOMOGRAPHY B07303
particular event is represented as the sum of two plane
waves of the form
U ðwÞ ¼ A ðwÞexp½(cid:7)iðk (cid:8)x(cid:7)wtÞ(cid:9)
z 1 1
þ A ðwÞexp½(cid:7)iðk (cid:8)x(cid:7)wtÞ(cid:9); ð1Þ
2 2
where U is vertical displacement, A is the amplitude of
z i
each incoming plane wave, k is the horizontal wave
i
number vector, x is the position vector, and t is time. Li et
al. [2003] showed that when this method was used in
Rayleigh wave tomography in eastern North America it
provided 30–40% variance reduction compared with the
standard one-plane wave method.
[12] To characterize the heterogeneous structure within
thenetwork,thetargetvolumeisparameterizedusingagrid
of nodes. The phase velocity is defined at each of these
nodes by
Vðw;qÞ ¼ B ðwÞ þ B ðwÞcosð2qÞ þ B ðwÞsinð2qÞ; ð2Þ
0 1 2
where B is the azimuthally averaged phase velocity, B are
0 i
the anisotropic phase velocity coefficients, q is the azimuth
ofpropagation, andwis frequency.Weassume that higher-
order azimuthal terms (4q terms) are small for Rayleigh
waves [Smith and Dahlen, 1973]. The direction of fast
propagation is 1 arctan (B / B ), and the peak-to-peak
2 2 1
amplitude or degree of anisotropy is 2(B2 + B2)1/2 / B .
1 2 0
We invert the frequency-dependent phase and amplitude
data separately for each period band.
[13] Because of their finite frequency, surface waves
are sensitive to two-dimensional structure near the propa-
gation path. To account for these effects we calculate two-
dimensional sensitivity kernels for fundamental Rayleigh
wavesbymeansofasingle-scattering(Born)approximation
[Zhou et al., 2004; Yang and Forsyth, 2006]. For each
frequency band, the phase and amplitude sensitivity kernels
Figure3. (a)VerticalseismogramsforstationPAYGforan arecalculatedforphasevelocityperturbationsandareincor-
eventthatoccurredintheVanuatuIslandson26November poratedintotheisotropicphasevelocityinversions(Figure4).
1999(MS7.3,epicentraldistance100(cid:1),backazimuth253(cid:1)). [14] The solution of the nonlinear inverse problem for
Unfiltered seismogram on top and bandpass-filtered Ray- planewaveandphasevelocityparametersisperformedasa
leighwavesforperiods20to125sbelow.(b)Rayleighwaves two-stage iteration process [Forsyth and Li, 2005]. In the
filtered at 50-s period, on vertical seismograms for all first stage of each iteration, velocity is held fixed and the
recordingstations.(c)SensitivitykernelsforRayleighwaves best fitting parameters for the two-wave approximation are
asfunctionsofdepthforperiodsof20,29,50,80,and125s. found for each using the downhill simplex method of
simulated annealing [Press et al., 1992]. In the second
stage, corrections to the velocity model and wave parame-
tersaredeterminedusingthelinearizedinversion technique
tivity in depth at about 1/3 of the wavelength (Figure 3c).
ofTarantolaandValette[1982].Theobserveddata,thereal
We measure the amplitude and phase of each bandpassed-
and imaginary components at a single frequency, are ini-
filteredseismogramusingthediscreteFouriertransform.To
tially assigned equal variance. Experience shows that a
ensure data quality we select Rayleigh waves having
typical misfit to the normalized real and imaginary terms
amplitudes at least 2.5 times greater than that of the
isontheorderof0.1,whichwechooseasaninitialapriori
preceding body waves. Furthermore, we use seismograms
estimate of standard deviation.
only from events for which the waveforms from station to
[15] We also assume that the solution (velocity parame-
station are similar, i.e., for which the average normalized
ters) is not too far from an initial estimate, so we penalize
cross-correlation coefficient is greater than 0.9.
changes from this starting model. This penalty is achieved
[11] Propagation effects outside the network, as well as
byintroducingnonzerotermsinthediagonalsoftheapriori
heterogeneous structure within the network, can affect
model covariance matrix. The amount of penalization is
Rayleigh waves. In order to account for wave-propagation
controlledbytheparameters ,whichistheapriorivalueof
effects outside the network, such as multipathing, we use a o
the standard deviation for the velocity terms in the inver-
two-plane wave approximation technique [Forsyth and Li,
sion. This parameter is an estimation of the allowed varia-
2005]. At each frequency w, the incoming wavefield of a
4 of25
B07303 VILLAGO´MEZ ETAL.: GALAPAGOS SURFACE WAVE TOMOGRAPHY B07303
Changes in the phase velocity of a Rayleigh wave are
mainly sensitive to perturbations in shear wave velocity
and less to perturbations in compressional wave velocity
(V )ordensity.Weperforminversionsforone-dimensional
P
V structure at each grid node by finding the best fit
S
between the observed phase velocities and those predicted
by the code DISPER80 [Saito, 1988], which calculates
normal modes for laterally homogeneous media. This tech-
nique yields predicted isotropic phase velocities from a
given shear wave velocity model, as well as sensitivity
kernelsforV ,V ,anddensity.Thesesensitivitykernelsare
S P
used in the inversion for V perturbations from a starting
S
model in an iterative process using the linearized inversion
technique of Tarantola and Valette [1982]. The inversion
resultsarevaluesofV asafunctionofdepthandestimates
S
of standard deviation.
[18] Because the inversion for shear wave velocity is
underdeterminedwemustassumesomeaprioriinformation
about the model parameters. We use an a priori model
covariance matrix of the form [e.g., Tarantola and Valette,
1982]:
Cm ¼ s2 exph(cid:7)(cid:1)D (cid:7)D(cid:2)2=(cid:1)2D2(cid:2)i; ð3Þ
ij i i j
where D is depth, D is the characteristic length of
smoothing, and s is the a priori estimate of the standard
i
deviation of the ith velocity term in the inversion. We
assume that the resulting shear wave velocities are not too
far from an initial estimate, so we penalize changes from
this starting model by introducing nonzero terms in the
diagonals of Cm. The amount of penalization is controlled
Figure 4. Normalized two-dimensional phase sensitivity bytheparameters.Alowervalueofs representsahigher
i i
kernelstolocalphasevelocityperturbationshowingthefirst penaltyandgreaterdampingofthesolution.Usingdifferent
negative and positive sensitivity regions for the event values of s for different layers allows us to constrain
i
depictedinFigure3andstationPAYGat(a)50-speriodand selectivelythedifferentpartsofthemodel.Additionally,we
(b) 80-s period. White triangles represent seismic stations, introducetheassumptionthattheresultingvelocitymodelis
the larger triangledenotes stationPAYG,andthe blackline smooth. We impose smoothness on the model for the ith
indicates the great circle path. parameterbypenalizingdifferencesinvelocitywithrespect
to neighboring points, through the introduction of non-zero
tions with respect to the starting model. In addition, we terms to the off-diagonals of the a priori model covariance
positionasetofnodessurroundingtheregionofinterest,for matrix weighted using the characteristic distance D.
which we allow more variation; this outer ring of nodes [19] A three-dimensional VS model is constructed by
absorbs additional traveltime variations not accounted for mergingalltheone-dimensionalVSresultsobtainedateach
by structure inside the target volume. node. This three-dimensional model is generally smoother
[16] To remove the influence of events that are not well in the vertical direction than laterally, and so we apply
described by the two-plane wave approximation, each horizontal smoothing within each depth layer using a two-
inversion is performed in two sets of iterations. In the first dimensional moving average of neighboring points. When
set, all the observations are assigned equal variance or smoothing laterally, we allow changes only up to a small
weight, as expressed above. Then, in the second set of fraction of the standard deviation, usually 10–20%. The
iterations, the observations are assigned variances based on resulting model is thus smooth in both the vertical and
the resulting standard deviations found after the first set lateral directions, and in our experience the maximum
of iterations. This sequence ensures that poor wavefield magnitude of a typical velocity anomaly is somewhat
models are given less weight and donot bias the inversion. decreased but its spatial extent is preserved.
To describe the quality of fit to the data we use the root-
mean square (RMS) misfit of the phase in seconds, which 4. Results
represents the misfit that is most directly related to travel-
times and the velocity structure. [20] We first present results concerning the validity of the
two-planewaveapproximationmethodandacomparisonwith
3.2. Inversion for V estimates of the direction of propagation obtained indepen-
s
[17] The estimates of isotropic phase velocity are in turn dently from polarization analysis. Second, we show results
used to constrain the shear wave velocity (V ) structure. of inversions for one- and two-dimensional phase velocity,
S
5 of25
B07303 VILLAGO´MEZ ETAL.: GALAPAGOS SURFACE WAVE TOMOGRAPHY B07303
for cases with and without anisotropy. Third, we present
images of three-dimensional V structure derived from the
S
two-dimensionalisotropicphasevelocityinversions.
4.1. Validity of the Two-Plane-Wave Approximation
[21] Results of the two-plane wave approximation show
thatingeneraltheprimarywaveismuchlargerinamplitude
than the secondary wave. The average ratio of the primary
waveamplitudetothesecondarywaveamplitudedecreases
with frequency from 7.9 at 8 mHz (125 s period) to 2.8 at
50 mHz (20 s period). This decrease is expected because
higher frequency waves are more strongly affected by
focusing and multipathing. Furthermore, deviations from
great circle path are less than 30(cid:1) for the primary waves.
[22] To test the validity of the two-plane wave approxi-
mation we compare the two-plane wave results with those
obtained independently from polarization analysis [Vidale,
1986]. To ensure measurement quality we use the cutoff
parameters of Larson and Ekstro¨m [2002] for Rayleigh
waves. Measurements of the direction of propagation of
Rayleigh waves using polarization analysis confirm that
deviations from great circle paths are small (less than 30(cid:1)).
Moreover, polarization analysis and the two-plane wave
approximation are generally in good agreement on the
primary direction of propagation (Figure 5).
[23] To investigate if local topography affects the propa-
gation of Rayleigh waves, we measure the scattering of the
arrivalanglesforalleventsateachstationfrompolarization
analysis. If there is a local topographic effect, there should
be noticeable scattering of the individual arrival angles.
Moreover, if topographic effects are important, scattering
should be frequency dependent, because higher-frequency
wavesaremoresensitivetotopographythanlower-frequency
waves. The amount of scattering is quantified using the
standard deviation of the individual measurements. We
found that the scattering of measurements is relatively
small (averaging 8.4 ± 5.1(cid:1) for all events), and that there
is no frequency dependence at the 95% confidence level.
We conclude that the effect of local topography on the
Figure5. Resultsfromthetwo-planewaveapproximation
propagation of Rayleigh waves is not significant and that
and polarization analysis. (a) Polarization direction of
the incoming wavefield can be accurately described by the
Rayleigh waves for the event shown in Figure 3. Filled
two-plane wave approximation.
white bars indicate the direction from polarization analysis
4.2. Phase Velocity Inversion and its uncertainty. Dashed lines show great circle paths,
[24] We present results from four sets of inversions for and solid lines indicate the direction of the primary wave
frequency-dependent phase velocity. In all inversions we fromthetwo-waveapproximation.(b)Deviationfromgreat
usearegulargridofnodesseparatedby0.2(cid:1)inlatitudeand circle path of the first plane wave from the two-wave
longitude (Figure 6a); this grid is encompassed by a set of approximation compared with the mean deviation from
nodes with larger prior uncertainties. Figure 6b shows the greatcirclepathfrompolarizationanalysisfora29-speriod
typicalpathcoverageusedforthephasevelocityinversions. and events with an amplituderatio of primary to secondary
Table 1 shows the number of events and observations and waves of greater than 4.
the resulting RMS phase misfits for all inversions.
[25] Inthefirstsetofinversionswesolvefortheisotropic consistentlylowerthanvaluesforPacificlithosphereofcom-
component of phase velocity, B , which is kept constant at
0 parable age [Nishimura and Forsyth, 1989], although this
all grid nodes. These results provide a uniform isotropic
differenceislesspronouncedatlongerperiods.Atperiodsof
phasevelocity for theentire regionforeachperiod.Weuse
20to67sthephasevelocityis2to2.5%lowerthanforPacific
aninitialvalueofphasevelocityof3.8km/sands =0.1as
o lithosphere0–4Myold(0–4NF89)and4.5to8%lowerthan
the a priori value of the standard deviation for the phase
for Pacific lithosphere 4–20 My old (4–20NF89). At longer
velocity.
periods,80to125s,phasevelocitiesare0to2%lowerthanfor
[26] Results of the inversion show that phase velocity 0–4NF89 and 2.5 to 4% lower than for 4–20NF89. Phase
increaseswithperiodfrom3.625±0.005to4.05±0.02km/s
velocities are 1–2% higher than values for Iceland [Li and
for periods from 20 to 125 s (Figure 7). These values are
Detrick, 2006], except at periods 25 to 50 s where they are
6 of25
B07303 VILLAGO´MEZ ETAL.: GALAPAGOS SURFACE WAVE TOMOGRAPHY B07303
provideanimportantconstraintonmantleflow.Finitestrain
induceslattice-preferredorientation(LPO)ofminerals,such
as the alignment of the olivine a axis [e.g., Christensen,
1984]. Because olivine, the most abundant mineral in the
upper mantle, is seismically anisotropic, the alignment of
crystallographic a axes caused by mantle flow can produce
measurable anisotropy.
[28] The peak-to-peak amplitude of azimuthal anisotropy
varies between 0.2 and 1% (0.3–0.5% standard deviation)
for periods 20–50 s, and between 1.2 and 3% (0.6–1.4%
standard deviation) for periods 67–125 s (Figure 8a). The
values at 100-s period have been omitted because of the
large uncertainties associated with the resulting parameters,
andbecausewedidnotachieveamisfitreductionrelativeto
the isotropic inversion (see Table 1). However, because of
the relatively high uncertainties we cannot reject the null
hypothesis of isotropy at the 95% confidence level, espe-
cially for periods 20–50 s. This result suggests that at
shallower depths the magnitude of the regional azimuthal
anisotropy is small or variable in direction so that effective
anisotropyislow.Significantseismicanisotropyisobserved
at periods longer than 50 s, indicating that its source is
likelylocatedatdepthsgreaterthan(cid:1)100km(seeFigure3c).
Forperiodsgreaterthan50s,adegreeofanisotropyofabout
1–3% agrees with regional estimates of Nishimura and
Forsyth[1988]thatshowazimuthalvariationsof1–2%.
[29] The direction of fast Rayleigh wave propagation is
generallyclosetoeast-west(73–101(cid:1)),comparablewiththe
easterly direction of Nazca plate motion in the hot spot
reference frame (90.1(cid:1) azimuth at 0(cid:1)N, 91(cid:1)W, for HS3-
NUVEL1A) [Gripp and Gordon, 2002]. At 25- and 29-s
period the fast direction of propagation changes to almost
N-S (14 ± 9(cid:1) and 9 ± 51(cid:1), respectively), close to the
direction of Nazca-Cocos spreading (7.15(cid:1) at 1(cid:1)N, 91(cid:1)W,
for NUVEL-1A) [DeMets et al., 1994]. However, because
only two period bands show this anomalous direction, and
becauseofthehighuncertaintyofthemeasuredazimuthfor
the 29-s band and the low degree of anisotropy at lesser
periods, we consider that the predominant direction of
azimuthal anisotropy in the region is east-west (E-W). We
could not resolve lateral variations of anisotropy, and thus
ourresultsareaverageestimatesofazimuthalanisotropyfor
the entire region, which includes the Gala´pagos platform
and its surroundings (Figure 6). However, SKS splitting
indicates that anisotropy within the Gala´pagos platform
Figure6. (a)Gridnodeparameterizationusedinthephase
varies laterally, with isotropy in the center of the archipel-
velocity inversions. (b) Path coverage for 50-s phase
ago and anisotropy with nearly E-W fast directions along
velocity inversion. White triangles denote seismic stations.
the western edge (81(cid:1)–109(cid:1) at seismic stations G05, G06,
G07, and G10) [Fontaine et al., 2005]. Regional observa-
tions of Rayleigh wave 2q azimuthal anisotropy across the
similar,andarecomparabletoyoungPacificlithospherenear
eastern Pacific also indicate an E-W fast direction of
theEastPacificRise(EPR)[Forsythetal.,1998](Figure7).
anisotropy [Nishimura and Forsyth, 1988; Montagner and
[27] In a second set of inversions, we add uniform
Tanimoto, 1990]. We suggest that the observed Rayleigh
azimuthal anisotropy and solve for the phase velocity
wave azimuthal anisotropy represents an average between
parameters (B , B , and B ), which are kept constant at all
0 1 2 an E-W direction of regional mantle flow and isotropy
gridnodes.Weuses =0.1asanapriorivalueofstandard
o beneath the center of the archipelago.
deviation for the velocity and anisotropy terms. Results
from the inversion show that the coefficient B changes by [30] In a third set of inversions, we obtain lateral varia-
0 tions in phase velocity, but we do not allow for azimuthal
less than 0.3% from the previous isotropic inversion. From
anisotropy.Wesolveforisotropicphasevelocity,B ,ateach
the coefficients B and B we obtain average values of the 0
1 2 grid node, while including two-dimensional sensitivity
fastdirectionofpropagationandamplitudeof2qanisotropy
kernels. We use the value of B from the uniform velocity
for the entire region. Measurements of seismic anisotropy 0
inversions as the initial value in the two-dimensional
7 of25
B07303 VILLAGO´MEZ ETAL.: GALAPAGOS SURFACE WAVE TOMOGRAPHY B07303
Table1. Comparison ofPhase Velocity Inversions
RMSPhaseMisfit(s)
Numberof Numberof UniformVelocity, UniformVelocity, Two-DimensionalVelocity,
Period(s) Events Observations NoAnisotropy UniformAnisotropy NoAnisotropy
20 94 1330 0.64 0.62 0.51
22 110 1498 0.58 0.53 0.51
25 123 1636 0.71 0.62 0.59
29 119 1584 0.89 0.78 0.69
33 118 1660 0.68 0.63 0.62
40 120 1564 0.79 0.72 0.68
50 120 1540 0.85 0.80 0.77
67 110 1394 0.93 0.87 0.86
80 88 1114 0.83 0.81 0.83
100 76 912 0.82 0.82 0.78
111 69 828 0.92 0.88 0.86
125 63 734 1.01 0.97 0.92
inversions. We also tested the use of a two-dimensional velocities vary laterally by up to ±1.5% with respect to the
perturbational model resulting from adopting the inversion uniform phase velocity model (Figure 9). However, phase
solutionatoneperiodastheinitialmodelfornearbyperiods velocities are consistently lower than values for Pacific
(e.g.,Weeraratneetal.,Rayleighwavetomographybeneath lithosphere of comparable ages: 0 to 3% lower than 0–
intraplate volcanic ridges in the South Pacific, submitted to 4NF89, and 2 to 9% lower than 4–20NF89. By examining
Journal of Geophysical Research, 2007),and weconfirmed the a posteriori model covariance matrix and the values of
that our results are independent of the starting model. The uncertainties inmodel parameters we define anareaof best
results of two-dimensional inversions show significant path coverage and resolution, which is used to plot the
improvement (up to 40% variance reduction) when com- phase velocity maps in Figure 9.
pared with the uniform isotropic and anisotropic phase [31] There aretwomainregionsofanomalouslylowphase
velocity inversions, suggesting that lateral variations of velocity. The first is near the southwestern corner of the
phasevelocityarerequiredbythedata.Theresultingphase archipelago,beneaththevolcanoesofFernandinaandsouthern
Figure7. AveragephasevelocityasafunctionofperiodfortheGala´pagosArchipelago(whitesquares
andboldline).DashedanddottedlinesindicateresultsfromthestudyofNishimuraandForsyth[1989]
forPacificOceanlithosphereofage0–4Myand4–20Myold,respectively.Circlesindicateresultsfrom
thestudyofLiandDetrick[2006]forIceland.TrianglesindicateresultsfromthestudyofForsythetal.
[1998] for the East Pacific Rise. All error bars represent one standard deviation.
8 of25
B07303 VILLAGO´MEZ ETAL.: GALAPAGOS SURFACE WAVE TOMOGRAPHY B07303
Figure8. Resultsfrominversionswithazimuthalanisotropy.(a)Amplitudeofanisotropyasafunction
of period and 1-s error bars. (b) Azimuth of fast direction of propagation and 1-s error bars. The solid
horizontallineindicatesthedirectionofplatemotioninahotspotreferenceframe(89.4(cid:1)at0(cid:1)N,91(cid:1)W,
for HS3-NUVEL1A) [Gripp and Gordon, 2002]. The dashed horizontal line indicates the direction of
Nazca-Cocos spreading (7.15(cid:1) at 1(cid:1)N, 91(cid:1)W, for NUVEL-1A) [DeMets et al., 1994].
Isabela. The anomaly is more evident at shorter periods, the procedure described in section 3.2. We apply lateral
especially20–25s.Thesecondlow-velocityregioniscentered smoothing while allowing for changes of less than 10% of
near 0.5(cid:1)S, 90.5(cid:1)W, beneathSantiago and Santa Cruz and is the standard deviation.
evidentinthephasevelocitymapsfrom29-to80-speriod.At [34] The smooth form of the phase velocity kernels
100-to125-speriod,thesecondanomalydecreasesinintensity (Figure 3c) shows that surface waves cannot resolve sharp
andmovesslightlysouthward. vertical velocity changes, including the expected variations
[32] In a fourth set of inversions, we assume uniform at the crust-mantle interface. Crustal velocity structure and
anisotropy but allow for lateral variation in phase velocity. thickness estimates are, however, available for the Gala´pa-
Thetwo-dimensionalphasevelocities(B )varybylessthat gos platform [Feighner and Richards, 1994; Toomey et al.,
0
0.3% compared with the isotropic case. The amplitude of 2001]. We tested the results of the inversion with three
anisotropy and the direction of fast propagation are also different assumptions about the crustal structure. Under the
verysimilartothoseobtainedintheinversionswithuniform firstassumptionweuseaconstantcrustalthicknessof15km
phase velocity: the direction of fast propagation varies andanaveragecrustalvelocityprofile[Toomeyetal.,2001]
between71and93(cid:1),exceptat25-sand29-speriods,where as our initial model. Changes with respect to these initial
the direction of anisotropy is 19 ± 7(cid:1) and 11 ± 74(cid:1), crustal velocities are penalized more than changes of
respectively, and the peak-to-peak amplitude of anisotropy mantle velocities in the inversion (s = 0.01 versus s =
i i
variesbetween0.15and2.9%.However,theresultsofthese 0.1,respectively).Becauseweexpectbathymetricdepthtobe
inversions do not provide significant variance reduction negatively correlated with crustal thickness, under the
with respect to the isotropic case. second assumption we assign different crustal thicknesses
to different grid nodes as a function of bathymetry: crustal
4.3. Shear Wave Velocity Inversion
thicknessistakentobe5kmifbathymetricdepthisgreater
[33] We usethetwo-dimensionalisotropicphasevelocities than 2000 m, 10 km if bathymetric depth is between 2000
(thirdsetofinversions)toconstructthree-dimensionalimages and 1000 m, and 15 km otherwise. Again, changes with
of shear wave velocity structure. Phase velocity data at each respect to the initial crustal velocities are penalized more
grid node are first inverted for one-dimensional VS. We then than changes of mantle velocities in the inversion (si =
parameterize a one-dimensional model (0- to 410-km depth) 0.01 versus s = 0.1, respectively). The third assumption is
i
in layers of 5-km thickness and use si = 0.1, D = 10 km. that crustal thickness and velocities are nowhere con-
Last, we merge all the resulting one-dimensional VS models strained (si = 0.1 everywhere). We find that changing
to obtain the three-dimensional velocity structure following the assumption about crustal structure has no significant
9 of25
B07303 VILLAGO´MEZ ETAL.: GALAPAGOS SURFACE WAVE TOMOGRAPHY B07303
Figure9. Results of inversion for two-dimensional isotropicphase velocity forallperiodbands. Units
are percent variation with respect to the frequency-dependent value of the isotropic uniform phase
velocity(fromFigure7).Contoursshownare(cid:7)1,(cid:7)0.5,0,and0.5%.ThicklinesoutlinetheGalapagos
Islands (0-m isobath), and white triangles indicate the locations of seismic stations.
10of 25