May 10, 2000 - William C. Hammond and Eugene D. Humphreys. Department of ... the Earth [Forsyth, 1992], the relationship between the physical state of the ...
JOURNAL
OF GEOPHYSICAL
RESEARCH,
VOL. 105, NO. B5, PAGES 10,975-10,986, MAY
10, 2000
Upper mantle seismic wave velocity' Effects of realistic partial melt geometries William C. Hammond and Eugene D. Humphreys Department of GeologicalSciences,University of Oregon, Eugene
Abstract. We investigate seismic wave velocity reduction resulting from the presence of partial melt in the upper mantle. The amount of shear and bulk modulus reduction produced by the presenceof a connectednetwork of realistically shaped and naturally organized melt inclusionsis found using finite element calculations. The geometriesof the inclusionsare taken directly from laboratory experiments of mantle melting, with finite element meshesconstructedto conform to these shapes. The shear and bulk moduli of the compositematerial are found for
both the unrelaxed(isolatedinclusions) and relaxed(pressure equalizedinclusions) casesby assigningappropriate material propertiesto the fluid. Modulus reduction from deformation simulations of a solid containing realistically shaped and ellipseshaped melt inclusions quantify the effect of melt pocket cuspatenessand melt pocket organization on seismicvelocity reduction. The three-dimensionalresponse is estimated from two-dimensionaldistributionsof the melt phase by determining the mode II and mode III componentsof elastic modulus reduction separately and summing their effects. In general, cuspate and naturally organized melt inclusions causegreater velocity reduction. It is shownthat V? and Vs reduction per percent
partial melt are at least 3.6% and 7.9%, respectively.Even highervaluesfor velocity reductionare possibleabove 1% melt fraction if melt existsonly in tubulesbelow 1% melt fraction. The lower, more conservativevalues of velocity reduction are •070% greater for Vp and 84% greaterfor Vs than the analyticallydetermined values for ellipsoidal inclusions. Somewhat greater effectsare possibleif nonrandom organization of melt occurson scalesgreater than our model.
1.
Introduction
contributingeffectsof variousmaterial qualities (e.g., temperature,composition)and to the complexityof the
Measurements of seismicvelocity and attenuation are responseof realistic partially molten material subject to frequently employedto infer the degreeand spatial ex- simple stresses. tent of partial melt in the Earth's uppermantle [e.g., The crux of the problem in understandingthe effects Zhao et al., 1992; Humphreysand Dueker, 1994; Sobolev of partial melt lies in understandinghow the microstrucet al., 1996; Xu and Wiens, 1997; Dunn and Toomey, tural geometry of the pore and conduit shapescontain1997; Toomeyet al., 1998]. While tomographicimaging the melt influencesrock elastic properties. Work to ing affords the highest resolution of any geophysical date has beenfocusedin two directions:(1) laboratory technique that images volumes of the deep interior of
the Earth [Forsyth,1992],the relationshipbetweenthe experimentaldeterminationof materialproperties[e.g., physicalstate of the mantle and the magnitude of seis- $ato et al., 1988; Hirth and Kohlstedt,1995a,b] and (2) analytic determinationof the mechanicalresponse mic velocity anomaliesremains imprecise. Recent studusing simplifiedgeometriessuchas cracks[O'Connell ies illustrate just how unconstrained the partial melt and Budianski, 1974],circularand three-grainjunction content of the upper mantle is (e.g., Sobolevet al. prismatic tubes [Mavko,1980],spheres [Einstein,1906], [1996],who find that for inferredP wavevelocityvariand ellipsoids [Eshelby, 1957; Wu, 1966; Schmeling, nt.icm.qc•f .•%. the nartial melt content can be anywhere 1985; Mainprice, i997]. Grea/. progre• ila• been made between0 and 3%). Improved understandingof the in laboratory determination of the mechanical response relationship between mantle physical state and seismic observationshas remained problematic because of the to mantle rocksboth above[i•.g.,Hirth and Kohlstedt, 1995a,b] and below[Jacksonet al., 1992]the solidus.
Paper number 2000JB900041.
These experiments,however,have not achievedthe temperature and pressure conditions that produce partial melting while also deformingthe sampleon time-scales
0148-0227/00 / 2000JB900041$09.00
that
Copyright 2000by the American Geophysical Union.
10,975
simulate
observed seismic waves.
10,976
HAMMOND
AND HUMPHREYS:
PARTIAL
MELT AND SEISMIC VELOCITY
Analytic results thus far have two limitations. First, times (frequencies)associatedwith the transition bethe shapes of melt inclusionshave been idealized as tweendiscretestatesof stress[O'C'onnelland Budiancracks, spheres,ellipsoids,or simplified cuspate forms ski, 1977]. Our strategyfor determinationof the viscoelastic mawhich only roughly approximate the true geometries of melt in the upper mantle [e.g., Waft and Bulau, terial propertiesof partially molten peridotirehasthree 1979; Cooper and Kohlstedt, 1982; Kohlstedt, 1992; Waft parts. It is based on the principle that the dispersion and Faul, 1992; Faul et al., 1994; Hirth and Kohlst- due to anelasticity can be derived from the stressre-
edt, 1995a]. While it hasbeenarguedthat the actual shapescan,in termsof summarystatistics(suchas aspect ratio and form factor), be adequatelyrepresented by ellipsoids[Schmeling, 1985;Faul et al., 1994],the mechanicalresponseof theseshapessystematicallydiffersfrom simpler geometrieswith the sameaspectratio. A fluid-filled, convexcircular cylinder, for example,is naturally more rigid than a cuspate,concavecylinder
sponseof the viscoelasticmaterial to a step function in strain. We calculate(1) the elasticresponsefor the "unrelaxed" casein which fluid flow between pores has not yet diminished the differential pore pressureexcited by the imposedstrain; (2) the elasticresponse for the "relaxed" case in which pore pressureequilibrium has been attained by fluid flow between pores; and (3) the relaxationtime spectrumby estimatingthe time-scalesof interpore flow and pressureequilibrium. Developmentof the anelasticprocessand the frequency dependentrelaxation, i.e., step 3, is discussedin the
with the samemelt fractionand aspectratio [Mavko, 1980].This leadsto an underestimation of the amount of velocity reduction for a given partial melt fraction paperby HammondandHumphreys [thisiswhenusingrelationshipsbasedon ellipsoidalgeometry. companion Second,determinationof the frequencydependence of sue],who determinethat relaxationoccurssufficiently the material's elastic modulus and attenuation has been quicklysothat it is the relaxedmodulusexcitedby obcomplicatedby the need to isolatespecificrelaxation servedseismicwaves. The first two stepsuse finite ele-
Figure 1. Backscattered electronmicroscopic imageof the melt phasefromlaboratoryexperimentby Faul et al. [1994].Blackrepresents quenched liquidmelt phase;whiteis crystalline solid.
HAMMOND
AND HUMPHREYS:
PARTIAL MELT AND SEISMIC VELOCITY
ment modelsdesignedto representrealisticcrystal-melt configurations, as inferredfrom laboratoryexperiments
of partially molten peridotite,as in Figure 1 [Faul et al., 1994].
10,977
Th ree-dimensio .n. al
b)
a)
The advantage of using finite element calculations in the first two steps is that any simple or complex pore shapecan be representedwith equal easeand elastic interactions
between distinct bodies are accounted
for. A disadvantageis that certain important intuitions are lost when not treating the problem analytically, so further analysis is required for a detailed understand-
ing of the underlyingprocesses.In theory, viscoelastic finite elements can be used to evaluate the details
of frequency dependencein seismicvelocity and attenuation. This is, however, computationally expensive sincelower-frequencyinteractionsrequire unreasonably
Two-dimensional
c)
d)
long run times. A lowerboundon the frequencies that can be investigatedpracticallyis thus imposed. These two problemsare addressedby usingan alternative approach, the network representationof Hammond and Humphreys[thisissue]. Creation of a fully three-dimensionalelastic finite element model of realistically shaped and organized melt inclusionsis difficultfor the followingreasons'(1) laboratorymantle melting experimentalproductsfrom which we derive inclusion shape data are polished fiat
to exposea singletwo-dimensionalcrosssectionfrom
x
which we must infer the three-dimensionalshapeof the
Figure 2. Idealized single melt. inclusionsillustratmelt phaseand (2) the crystalsthat surroundthe melt ing the inference of three-dimensional behavior from are too complexto be meshedby our finite element two-dimensional melt inclusionshapedata. (a) Threemeshingalgorithms. Becauseof these difficultieswe dimensionaloblate ellipsoid;arrows showsenseof shear usetwo finite element analysesthat use simpler geome- in the xy plane. (b) Three-dimensional cuspateinclutries to approximatethe full three-dimensionalelastic sion with same aspect ratio. In Figures 2a and 2b responseto inclusioncuspateness. The first calculation the edges excited by mode II and mode III cracktip determines the effect of cuspatenessand organization stressare shownby II and III, respectively.(c) Twocylinder. (d) Cuson the mode II componentof modulusreduction(see dimensionalellipsecross-sectioned pate cross-sectionedcylinder. In Figures 2c and 2d Figure 2). The seconddeterminesthe effectof mode shear sensearrows show mode II excited by shear in III modulus reduction owing to inclusion cuspateness xy plane and mode III excited by shear in yz and xz and organization. The total shear modulus reduction planes. is the sum of the mode II and mode III
effects.
From
the resulting reduced elastic moduli we calculate the amount of Vp and V$ velocity reduction per percent partial melt and compare it to the velocity reduction calculated from ellipse-shapedmelt inclusionswith the same melt fraction and aspect ratio. The ratio between these two is a seismicvelocity correction factor." The corrected velocities provide new values for the partial derivatives of seismicvelocity with respect to melt fraction. These valuescan be applied in the internrof, nf,icm c•f •qoi•qrnict,omogra,nhv
in the same manner
as
the derivatives of seismicvelocity with respect to temperature, density, and composition. 2.
Elastic
Element
Effects-
Finite
Calculations
Realistic melt geometrieshave been determined in
Kohlstedt,1995a]that simulatemantle conditionsto depths of -•70 km. We use imagesof the melt phase
from Faul et al. [1994](Figure 1) to createfinite element representationsof the compositematerial. Details of the experimental configuration and imaging proce-
duresare givenby Faul et al. [1994]and Waft and Bulau [1982]. The melt distributionresultingfrom these experiments is determined by slicing the run product and by polishingand imaging the glassyinclusionsus-
ing backscatterelectronmicroscopy.Once the inclusion geometryis in suitable digital format, the shapes are represented in two-dimensional finite elements(Figure 3). We determine
the unrelaxed
and relaxed shear mod-
uli with an elastic calculation by using the appropriate elastic constants for elements inside the inclusions.
laboratory experiments[e.g., Waft and Bulau, 1979; See Table 1 for specificvalues. To simulate the unreWaft and Faul, 1992; Faul et al., 1994; Hirth and
laxed case,the shearmodulusof the melt is assigneda
10,978
HAMMOND
AND HUMPHREYS'
PARTIAL MELT AND SEISMIC VELOCITY
(a)
to macroscopicshear strain is zero, we can simulatethis stressstate by assigningshear and bulk elastic moduli to essentiallyzero inside the inclusions,effectivelyfill-
ing the inclusionswith vacuum[Mavko,1980]. This technique is appropriate for modeling shear stressrelaxation from melt squirt but not for bulk relaxation.
It is arguedby Hammondand Humphreys[thisissue], however,that bulk effectsare probably negligiblecompared to shear effects. Thus we model the effects of melt squirt arising only from imposedshear. The
finite
element
calculations
determine
the dis-
placements to all nodal positions that minimize the globalelasticstrain energy. Stressesand strainsare calculated everywherewith element shape functionsand elastic constitutive relations subject to the displacements applied as boundary conditions. Figure 4 is a close up view of the hydrostatic strains near a threegrain junction resultingfrom a macroscopically applied shear displacementboundary condition. For the twodimensional results presented here, isoparametric triangular quadratic elementswith isotropic solid elasticity were used, while for the three-dimensionalcalculations isoparametricwedge-shapedlinear elementswere
used[e.g.,Zienkiewiczand Taylor,1989].Meshgeneration is accomplishedthrough a semiautomaticroutine. All the finite element calculations presentedhere were performedwith the commercialsoftwarepackageMSC
PATRAN/AdvancedFEA, a product of the MacNeal SchwendlerCorporation. 2.1.
Numerical
Accuracy
To verify that our finite element modelingprovides accurate results, two tests were performed. First, the nodal densitywas progressivelyincreaseduntil displacements and stressesdid not change significantlyeven around the complicated cuspate inclusion tips and in areasof inclusioninteraction. Second,a relatively complex case for which an analytic solution exists was
tested. Mavko [1980]solvedfor the amountof elastic modulus reduction due to the presenceof a cuspate cylindrical pore of "triple junction" crosssectionshape.
Figure 3. (a) Partial melt geometriesusedto calcu- For a singledry pore (i.e., an inclusionwith no fluid inlate compositematerial elasticproperties. (b) Finite side) with shapeparameter e - 0 (the most cuspate element mesh built around same geometrieswith elaspossible),length d, radius R, solid shearmodulus tic finite elements. Poisson'sratio rs, bulk modulus Ks, and total volume
V, the bulk modulusof the compositeK• is given(after correctionfor a minor error from Appendix A of Mavko
value that is essentiallyzero, while its bulk modulus is approximately that of a basaltic liquid. The response thus represents times after shear stressesin the fluid Table 1. Elastic Constants Used in Calculation of Mahave relaxed but hydrostatic stressesremain. When the terial Properties. external strain is imposed, the pressureinside each inSolid Melt clusion,which is a function of its shapeand orientation, is calculatedby the finite elementalgorithm. Unrelaxed Relaxed For the relaxed case we assume that all inclusions are
interconnectedand that fluid flow has occurred,equalizing the inclusionpressures.Sincethe differencesin pressure between inclusionsare identically zero and since the averagehydrostaticstressthroughoutthe soliddue
In GPa.
124
40
0
64
0
0
HAMMONDAND HUMPHREYS:PARTIALMELTAND SEISMICVELOCITY
.
Maximum Dilation
0
10,979
.
Maximum Compression
Figure 4. An enlarged example resultof thefiniteelement calculation showing thehydrostatic component of strainin the unrelaxed state. Macroscopic displacement movestop of meshto the right,whilebottomremains fixed.Whiteoutlinedelineates location of melt. Darkshades showcompressire strains,lightshades showexpansive strains.Theexistence of stressgradients in the melt indicatesthat the melt is unrelaxed.Signof the strainis a functionof the inclusion orientation.
[1980])by
the boundaryby makingthe volumelarge enoughso that stressesare closeto uniform acrosseach face of the
cube,yetsmallenough sothat theporeintroduces su•1 1 •rR•d 1 2•%]. (1) cient weakness to the composite. The stresses measured •K•= K•4-4•u•V'(1 - 2v•)[12 4-14-v• on the surfacesare averagedovereachface.
The effective bulk modulus calculated from finite elUsingu, = 0.25, K, = 46.67 GPa, /•, - 28.0 GPa, ements is the simple ratio of calculatedhydrostatic R = 1, V = 1000(the solidis a 10x 10 x 10 cube),and stress to imposedstrain, Kct = 45.567 GPa, an error d = 3 yieldsKa = 45.56 GPa. of
