A SEISMOLOGICAL PERSPECTIVE ON THE LITHOSPHERE-ASTHENOSPHERE BOUNDARY BY HEATHER ANNE FORD B.S., UNIVERSITY OF MICHIGAN, 2005 M.S., BROWN UNIVERSITY, 2009 A DISSERTATION SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY IN THE DEPARTMENT OF GEOLOGICAL SCIENCES AT BROWN UNIVERSTY PROVIDENCE, RHODE ISLAND MAY 2013 ! ! ! ! ! © Copyright 2013 by Heather A. Ford This  

dissertation  

by  

Heather  

Ford  

is  

accepted  

in  

its  

present  

form  

by  

the  

Department  

of Geological  

Sciences  

as  

satisfying  

the  

dissertation  

requirements  

of the  

degree  

of  

Doctor  

of  

Philosophy. Date  

_______________  

 __________________________ Karen  

Fischer Recommended  

to  

the  

Graduate  

Council Date  

_____________  

 __________________________ Donald  

Forsyth Date  

_____________ __________________________ Greg  

Hirth Date  

_____________ __________________________ Alberto  

Saal Date  

_____________ __________________________ Geoff  

Abers Approved  

by  

the  

Graduate  

School  

and  

Research  

 Date  

_____________ ______________________________ Dr.  

Peter  

Weber Dean  

of  

the  

Graduate  

School CURRICULUM VITAE Heather Anne Ford Born: Heather Anne Whittington December 26th, 1982 Dearborn, MI, USA Education 2013 Brown University, Ph.D. in Geological Sciences 2009 Brown University, M.Sc. in Geological Sciences 2005 University of Michigan, Honors B.S. in Geological Sciences 2005 University of Michigan, Academic Minor in Physics Research Experience 2013- Postdoctoral Researcher, Yale University • Investigate possible causes of lower mantle flow through shear-wave splitting analysis and forward modeling using the spectral element method with Maureen Long. 2007-2013 Graduate Research Assistant, Brown University • Studied structure and origin of lithosphere-asthenosphere boundary and intra- lithospheric discontinuities within continents through the use of Sp receiver functions, forward modeling and extrapolation of viscosities through shear- velocities with Karen Fischer (advisor) and others. • Helped to design a semi-automated MATLAB GUI package for receiver function calculations. • Mentored undergraduate and graduate students on receiver function generation and interpretation. • Assisted in the demobilization of broadband seismic stations of the High Lava Plains seismic experiment. 2002-2005 Research Assistant, University of Michigan • Produced cross sections using MDEQ well log database, collected and analyzed groundwater samples for element chemistry with Lynn Walter. Teaching Experience 2008 Teaching Assistant, Brown University Geological Sciences 220: Physical Processes in Geology- Assisted in teaching an introductory course in geology. Tasks included leading laboratory sections and field trip groups and grading laboratory, homework and field trip assignments (~100 students). Publications Ford, H.A., K.M. Fischer, D.L. Abt, C.A. Rychert, L.T. Elkins-Tanton, The lithosphere- asthenosphere boundary and cratonic lithospheric layering beneath Australia from Sp wave imaging, Earth and Planetary Science Letters, 300, 299-310, doi:10.1016/j.epsl.2010.10.007, 2010. iv# Fischer, K.M., H.A. Ford, D.L. Abt, C.A. Rychert, The lithosphere-asthenosphere boundary, Annual Review of Earth and Planetary Sciences, 38, 551-575, doi:10.1146/annurev-earth-040809-152438, 2010. Abt, D.L., K.M. Fischer, S.W. French, H.A. Ford, H. Yuan, and B. Romanowicz, North American lithospheric discontinuity structure imaged by Ps and Sp receiver functions, J. Geophys. Res., 115, B09301, doi:10.1029/2009JB006914, 2010. Publications in progress Ford, H.A., K.M. Fischer and V. Lekic, Scattered#Seismic#Waves#Map#the#Deep# Expression#of#the#San#Andreas#Fault, submitted, Nature Ford, H.A. and K.M. Fischer, Anisotropic Layering Within the Cratonic Lithosphere Imaged by Scattered Body Waves: Towards a Better Understanding of Craton Formation, in preparation for Geophys. J. Int. Ford, H.A., K.M. Fischer and G. Hirth, A Comparison of the Seismological and Rheological Lithosphere-Asthenosphere Boundary, in preparation for Earth Planet. Sci. Lett. Conference Abstracts and Presentations Ford, H.A., K.M. Fischer and V. Lekic, Significant Variations in the Strength of the Lithosphere-Asthenosphere Boundary Across the California Plate Boundary, Abstract T53D-03, presented at 2012 Fall Meeting, AGU, San Francisco, Calif., 3-7 Dec. (talk) Fischer, K.M., H.A. Ford, V. Lekic and G. Hirth, Rheological implications of the seismological lithosphere-asthenosphere boundary, Abstract T32C-05, presented at 2012 Fall Meeting, AGU, San Francisco, Calif., 3-7 Dec. Hopper, E., H.A. Ford, K.M. Fischer, V. Lekic, and M.J. Fouch, How has magmatism in the northwest United States affected the lithosphere? Insights from Sp Receiver Functions, Abstract DI21A-2349, presented at 2012 Fall Meeting, AGU, San Francisco, Calif., 3-7 Dec. Ford, H.A., K.M. Fischer and V. Lekic, Variations in Lithosphere-Asthenosphere Boundary Strength Across the California Plate Boundary, Incorporated Research Institutions for Seismology Workshop, Boise, ID, June 13-15, 2012. (poster) Ford, H.A., K.M. Fischer, H. Yuan, and B.A. Romanowicz, Characterizing anisotropic lithospheric layering and the lithosphere-asthenosphere boundary in cratonic North America with Sp receiver functions, Abstract DI51C-02 presented at 2011 Fall Meeting, AGU, San Francisco, Calif., 5-9 Dec., 2011. (talk) Ford, H.A., K.M. Fischer, D.L. Abt, H. Yuan, and B.A. Romanowicz, The structure of the cratonic lithosphere beneath North America imaged by Sp receiver functions, Interior of the the Earth Gordon Research Conference, Mount Holyoke College, MA, 5-10 June, 2011. (poster) Yuan, H., B.A. Romanowicz, H.A. Ford and K.M. Fischer, Anisotropy in the North America Craton, Paper No. 175-7, GSA Annual Meeting, Minneapolis, MN, 9-12 October 2011. Fischer, K.M., H.A. Ford, D. Abt, H. Yuan, and B.A. Romanowicz, The lithosphere- asthenosphere boundary and cratonic lithospheric layering beneath stable North v# America, Abstract T31F-02, presented at 2010 Fall Meeting, AGU, San Francisco, Calif., 13-17 Dec., 2010. Fischer, K.M., H.A. Ford, V. Lekic, and D.L. Abt, The lithosphere-asthenosphere boundary beneath North America and Australia, Abstract D132A-06, presented at 2010 Fall Meeting, AGU, San Francisco, Calif., 13-17 Dec., 2010. Romanowicz, B.A., H. Yuan, H.A. Ford, K.M. Fischer, and D. Abt, Stratification of Azimuthal Anistropy in the North America Craton, Abstract T31F-03, presented at 2010 Fall Meeting, AGU, San Francisco, Calif., 13-17 Dec., 2010. Ford, H.A., K.M. Fischer, D.L. Abt, C.A. Rychert and L.T. Elkins-Tanton, The Lithosphere-Asthenosphere Boundary and Cratonic Lithospheric Layering Beneath Australia from Sp Wave Imaging, CIDER Workshop, UC Santa Barbara, CA, 27 June – 17 July, 2010. Ford, H.A., K.M. Fischer, D.L. Abt, C.A. Rychert and L.T. Elkins-Tanton, Imaging the Lithosphere-Asthenosphere Boundary Beneath Australia using Sp Scattered Wave Receiver Functions, IPRPI Workshop on Seismic Tomography, Troy, NY, 27-28 April, 2010. (poster) Ford, H.A., K.M. Fischer, D.L. Abt, L.T. Elkins-Tanton, and C.A. Rychert, The lithosphere-asthenosphere boundary beneath Australia imaged by Sp phases, Eos Trans. AGU, 90(52), Fall Meet. Suppl., Abstract DI11A-05, 2009. (talk) Ford, H.A, K.M. Fischer, and C.A. Rychert, Scattered wave imaging of the lithosphere- asthenosphere boundary in Australia, Earthscope National Meeting, Boise, ID, 12- 15 May, 2009. (poster) Ford, H.A, K.M. Fischer, and C.A. Rychert, Characterizing Lithospheric Thickness in Australia using Ps and Sp Scattered Waves, Eos Trans. AGU, 89(53), Fall Meet. Suppl., Abstract U43B-0059, 2008. (poster) Fischer, K.M., D.L. Abt, H.A.Ford, S.W. French and C.A. Rychert, Imaging the Lithosphere-Asthenosphere Boundary with Scattered Waves, Eos Trans. AGU, 89(53), Fall Meet. Suppl., Abstract DI21A-1740, 2008. Awards Dissertation Fellowship, Brown University Camp Davis Field Geologist Award, University of Michigan Lucile B. Conger Alumnae Group Scholarship, University of Michigan Service, Talks and Activities 2012 Session chair of “Seismic Anisotropy and Plate Tectonics: Insight Through Observations, Experiments, and Modeling”, 2012 AGU Fall Meeting 2012 "Variations in Lithosphere-Asthenosphere Boundary Amplitude Across the California Plate Boundary", presented at Solid Earth Seminar, Boston University 2010-2012 Faculty Representative, Dept. of Geological Sciences, Brown University 2009-2010 SoLED Seminar Organizer, Dept. of Geological Sciences, Brown University 2008-2009 Geoclub Treasurer, Dept. of Geological Sciences, Brown University Other: vi# • Reviewer (2) for Geophysical Journal International • Attended workshop series on Reflective Teaching Practices, Harriet W. Sheridan Center for Teaching and Learning, Brown University • Mentored undergraduate Marsella Kachingwa, and graduate student Emily Hopper Professional Memberships American Geophysical Union Seismological Society of America Skills Foreign Languages: Basic skills in German Programming: FORTRAN, MATLAB Operating Systems: UNIX, LINUX, OS X, Windows Other work experience 2005-2007 Engineer and Surveyor, Northwest Consultants, Inc. • Duties included surveying, road design, and construction inspection. • Major projects included intersection design and surveying at Detroit Metro Airport. • Helped to prepare a feasibility study analyzing conflicts between utility facilities and a potential new bridge crossing between Detroit and Windsor. vii# Preface What is the “lithosphere-asthenosphere boundary”? The lithosphere-asthenosphere boundary is the world’s most ubiquitous plate boundary system, separating the high viscosity lithosphere from the underlying convecting mantle asthenosphere. Thoroughly characterizing the properties of the lithosphere- asthenosphere boundary can help us to improve our understanding of the behavior of tectonic plates and the attendant issues surrounding plate tectonic theory, mantle dynamics and the evolution of continents and oceans. To date, there is still considerable uncertainty in our understanding of many aspects of the lithosphere-asthenosphere boundary. Debate currently exists over the depth and gradient of the boundary, as well as how the boundary varies across tectonic environments and with age. The lithosphere-asthenosphere boundary has often been defined as a thermal boundary layer, where heat transport transitions from conduction (in the lithosphere) to convection (in the asthenosphere). However, other lines of evidence suggest that the lithosphere-asthenosphere boundary can be defined as a chemical boundary layer, characterized by a change in composition and/or hydration. It has also been suggested that the difference in rheological properties may stem from a change in melt content. ! viii! Work to constrain which physical properties are responsible for the contrast in rheology is often conducted using geophysical or petrological data. Heat flow can be used to infer the location of the thermal boundary layer, discussed above, while mantle xenoliths can provide information about both the geothermal structure of a region as well as composition, and changes in electrical resistivity are often thought to be indicative of a change in hydration or melt content. Seismic wave velocities are dependent on a number of properties, including temperature, composition, hydration, melt content and grain size and are commonly used to make inferences about the physical state of the upper mantle, including the lithosphere-asthenosphere boundary, which is characterized by a transition from fast velocities in the lithosphere to slow velocities in the asthenosphere. The most fundamental questions that this dissertation explores are 1) what are the seismic properties of the lithosphere-asthenosphere boundary and 2) what do these results imply about the physical properties of the mantle at the boundary? The bulk of the work presented in this dissertation, and used to address the two questions above, was produced using Sp scattered wave receiver function analysis. This method has become an increasingly popular tool for imaging the lithosphere-asthenosphere boundary because of its sensitivity to variations in gradient structure, which can in turn be used to differentiate between candidate physical mechanisms. In chapter one, we use Sp scattered wave receiver function analysis to image and characterize the seismic and physical properties of the lithosphere-asthenosphere boundary beneath Australia. Variations in the seismic structure of the lithosphere- ! ix! asthenosphere boundary are explored across three regions of distinct tectonic age. Forward modeling of Sp receiver function analysis constrains the range of possible velocity gradients, and in the process, helps us to better understand the physical properties of the lithosphere-asthenosphere boundary beneath Australia. Chapter two addresses a different set of questions than those listed above, although the results have bearing on how we utilize Sp receiver functions to image discontinuity structure, such as the lithosphere-asthenosphere boundary. More specifically, this chapter addresses the effects of seismic anisotropy, which is thought to be an important property of the upper mantle, on Sp receiver function analysis. We find that the effect of anisotropy on Sp receiver functions is complicated, and is a function of both the orientation and type of anisotropy present in the earth and the polarization of the incident S-wave. While much of this dissertation assumes that velocity structure is a direct proxy for viscosity structure in the upper mantle, chapter three directly explores the relationship between the two beneath cratonic North America. To do this, shear wave velocity profiles from a number of surface wave tomography models are mapped into viscosity by considering the effects of temperature on velocity and then using olivine flow laws to translate temperature into viscosity. In the process, assumptions regarding cratonic lithosphere composition are explored. ! x! Chapter 4 (and the associated Supplementary Materials in Appendix A) explores systematic variations in the strength of the lithosphere-asthenosphere boundary phase across the San Andreas fault system in California. The findings discussed in this chapter have important implications not just for the nature of the lithosphere-asthenosphere boundary, but for our understanding of the extent and geometry of the San Andreas fault system within the deep mantle lithosphere. Unwrapping the secrets of the cratons Cratons are a unique feature of continents in large part because of their longevity. Although they are often described as “boring” because of their imperviousness to recent tectonic activity, their presence provides us with one of the few opportunities to look at the evolution of earth. An important secondary theme of this dissertation is what the seismic structure within, and at the base of, cratons might tell us about the formation, evolution and preservation of cratonic lithosphere. This theme was not obvious at the outset, but arose organically from work conducted in the first chapter, where evidence for discontinuity structure observed within the cratonic lithosphere prompted us to ask questions about the formation and evolution of the cratons in central and western Australia. Evidence of similar discontinuity structure in North America, in conjunction with observed changes in azimuthal anisotropy at similar depths from surface wave tomography, prompted us to ask if the discontinuity structure ! xi! observed in North America is the product of a boundary in anisotropy. The results of this analysis are discussed in chapter two. One of the most enigmatic aspects of cratonic lithosphere is its ability to remain stable, despite being cold and dense. Chemical depletion is often invoked to help explain the stability, since increased depletion reduces density. Understanding the effect these parameters have on seismic velocity has also helped to explain why velocities within the cratonic mantle lithosphere are so high. However, in several cases, such as demonstrated in chapter three, the observed velocities are larger than what is expected when accounting for variations in temperature and composition. This suggests that while our understanding of the cratons is improving with time, there are still significant gaps in our knowledge of the properties of the cratonic lithosphere. ! xii! Acknowledgements “To my teachers who became friends, to the friends who became my teachers, to all my family who have always been my best friends and teachers” -The Lithosphere, Irina Artemieva ! ! My! time! as! a! graduate! student! has! been! marked! by! a! lot! of! learning,! new! experiences!and!individual!growth.!!Luckily,!I!have!also!found!myself!surrounded!by! some!of!the!most!caring,!smart,!and!hard>working!people!imaginable.!!Without!the! support! of! my! friends,! teachers,! and! family,! graduate! school! would! have! been! a! lonely!and!pointless!journey.! ! First,!I!would!like!to!thank!my!graduate!student!family!for!making!Rhode!Island!feel! like!home.!!I!cannot!possibly!name!everyone,!but!please!know!that!I!have!enjoyed!all! of!the!conversations,!parties!and!trips!we!shared!over!the!past!six!years.!!A!special! thank!you!to!my!fellow!classmates,!including!those!who!I!matriculated!with>!Debra! Hurwitz,! Brendan! Hermalyn,! J.R.! Skok,! Li! Gao,! Susie! Theroux! and! Uly! Horodyskyj,! for! making! the! first! few! years! of! graduate! school! a! blast,! and! for! continuing! to! inspire!me!to!work!my!hardest.!!I!am!also!incredibly!grateful!to!have!gotten!to!know! Rachel! (and! Allen)! Klima,! Peter! Isaacson,! Ashley! Nagle,! Kate! Burgess,! Bronwen! Konecky,!Shannon!Loomis,!Jessica!Rodysill,!Jeff!Salacup,!Mark!Salvatore,!Leah!Cheek,! Mary!Peterson!(Go!Blue!),!Colin!Jackson!and!Tim!Goudge.!!! ! It is impossible for me to fully express my thanks to my SoLED friends. They are the ones who sat in class, discussed homework and worked in the lab with me, and have become my closest friends. I would like to acknowledge Mariela Salas, Yun Wang, Amandine Cagnioncle, Christine McCarthy and Marshall Sundberg for giving helpful advice, academic and personal, and for helping to recruit me to come to Brown. A special thank you to Linda (Chernak) Meyers, for becoming a good friend and inspiration, to David Abt for teaching me everything I know about MATLAB and to Tina Rau and McCall Burau for always being there for me. During my last few years at Brown, I have been lucky to become acquainted with Julia MacDougall and Emily Hopper, whose presence make the seismology lab the place to be, as well as Amanda Getsinger, Chris Havlin, Ved Lekic and Margarete Jadamec, who give excellent advice and support. Finally, a special, special thank you to Angela Stickle, who has become one of my best friends. ! Next,!I!would!like!to!acknowledge!the!faculty!in!the!geological!sciences!department.!! To! start! with,! I! want! to! extend! my! most! sincere! thanks! to! my! wonderful! advisor,! Karen! Fischer.! ! Karen’s! knowledge,! support! and! encouragement! have! been! indispensible!over!the!years!and!I!am!incredibly!lucky!to!have!gotten!to!know!and! ! xiii! work!with!her,!and!hope!to!grow!closer!to!her!in!the!years!to!come.!!To!my!advisory! committee,! including! Don! Forsyth! and! Greg! Hirth,! I! owe! an! incredible! debt! of! gratitude,! for! teaching! me! much! of! what! I! know! today,! and! for! being! sources! of! support! as! I! move! forward! in! my! career.! ! A! special! thanks! goes! to! Steve! Parman,! Alberto!Saal,!Marc!Parmentier,!and!Geoff!Abers,!for!providing!me!with!constructive! feedback! and! for! serving! on! my! various! committees.! ! I! would! also! like! to! thank! Terry! Tullis,! Reid! Cooper,! Jan! Tullis,! Meredith! Hastings! and! David! Goldsby,! for! providing! me! with! new! ways! of! thinking! about! my! research! as! well! as! for! the! frequent!informal!conversations!in!the!hallways.! ! I!would!be!remiss!if!I!did!not!mention!the!incredible!support!I!have!received!from! the! department! staff! over! the! years.! ! Both! Margaret! Doll! and! Bill! Collins! deserve! special! acknowledgement! for! their! ability! to! help! me! out! of! a! jam,! no! matter! how! large!(losing!everything!on!backup)!or!small!(turning!up!the!heat!in!the!lab).!!I!am! also! very! grateful! for! the! hard! work! and! dedication! of! the! staff! in! the! main! office,! including!Nancy,!Pat,!Ruth,!Caroline,!Paula,!Gloria!and!Lisa.! ! My! passion! for! geology! became! a! full>blown! obsession! during! my! undergraduate! education.! ! Those! who! helped! to! prepare! me! for! graduate! school! include! Bruce! Wilkinson,! Henry! Pollack,! Kyger! Lohmann,! and! the! late! Eric! Essene.! Lynn! Walter! and!Philip!Meyers!helped!to!mentor!and!guide!me!in!my!research!endeavors,!as!did! Jennifer! McIntosh,! and! Peter! van! Keken,! who! should! be! thanked! for! giving! me! the! courage!to!switch!career!paths!and!study!geophysics.!!I!would!also!like!to!thank!the! undergraduate!students!who!I!worked!along!side!of!in!class,!lab!and!in!the!field.!!! ! My! life! would! be! incomplete! if! I! did! not! have! friends! outside! of! academia.! ! Their! support!and!love!mean!the!world!to!me!and!I!hope!to!repay!it!back!to!them!some! day.! ! In! particular,! I! would! like! to! thank! Jackie,! Bridget! and! the! rest! of! my! high! school!cross>country!friends!for!lending!me!courage!and!strength!when!I!needed!it! the!most.!Also,!and!in!what!may!be!a!dissertation!first,!I!would!like!to!acknowledge!a! group!of!friends,!known!as!BitofGlitter,!most!of!whom!I!have!never!met!in!real!life,! but!are!always!there!to!listen!and!know!what!to!say.!!! ! I!would!like!to!acknowledge!my!amazing!family!that!I!have!been!blessed!to!be!a!part! of,! and! who! have! been! the! most! steadfast! supporters! in! my! endeavors.! ! To! start,! I! would!like!to!thank!my!in>laws,!David!and!Susan!Ford,!for!their!love,!generosity!and! words!of!support.!!I!would!also!like!the!thank!my!Grandparents,!Barbara!and!Bobby! Whittington,!and!Carolyn!Ray,!for!putting!up!with!my!absent>mindedness!and!loving! me! from! afar,! and! a! very! special! thank! you! to! the! late! Ernest! Ray,! for! his! sense! of! adventure!and!love!of!the!outdoors,!which!I!seem!to!have!inherited.!!!In!the!process! of!moving!away!from!home,!I!find!that!I!have!become!incredibly!close!to!my!sister,! Melanie!Dove,!who!always!knows!what!to!say!and!my!brother,!Daniel!Whittington,! who!has!accompanied!me!on!many!rock>hunting!(and!climbing)!adventures.! ! I!would!also!like!to!acknowledge!my!parents,!Christine!and!Terry!Whittington!and! my!husband,!David!Ford.!!Without!their!support,!I!would!have!never!embarked!upon! ! xiv! this!journey,!and!I!certainly!would!not!have!succeeded.!!My!parents!have!taught!me! many!things,!including!how!to!love,!the!importance!of!working!hard,!and!doing!what! you! love.! ! Today,! they,! along! with! Dave,! are! the! first! ones! I! go! to! for! advice! and! words!of!encouragement.!!My!husband!Dave!has!been!my!most!constant!companion,! and!I!am!so!incredibly!lucky!to!have!him!as!my!spouse!and!my!best!friend.!!Without! him,!my!life!would!be!incomplete.! ! Finally,!I!would!like!to!dedicate!this!dissertation!to!my!son,!David,!who!has!taught! me! an! incredible! amount! about! myself! and! about! life! over! the! past! two! years.!! Completing!graduate!school!with!a!child!is!never!easy,!but!it!has!certainly!made!the! process!memorable.!!! !!! ! ! ! xv! Table of Contents Curriculum Vitae ........................................................................................ iv Preface ........................................................................................................ viii Acknowledgements .................................................................................... xiii Chapter 1. The lithosphere-asthenosphere boundary and cratonic lithospheric layering beneath Australia from Sp wave imaging .............. 1 Abstract...........................................................................................................................2 1. Introduction ................................................................................................................3 1.1 Overview of Australian tectonics...............................................................................3 1.2 Results from surface wave tomography .....................................................................4 2. Method ........................................................................................................................6 2.1 Ps and Sp receiver functions ......................................................................................6 2.2 Data preparation .........................................................................................................7 2.3 Deconvolution and migration ....................................................................................8 3. Results .......................................................................................................................10 3.1 Phanerozoic Australia ..............................................................................................12 3.2 Cratonic Australia ....................................................................................................13 3.3 Large-scale correlation in depth and amplitude estimates .......................................15 4. Discussion .................................................................................................................16 4.1 Comparison to other receiver function studies ........................................................16 4.2 Estimating the gradient thickness and velocity contrast at a boundary ...................17 4.3 The lithosphere-asthenosphere boundary in eastern Australia ................................18 4.4 The lithosphere-asthenosphere boundary in central and western Australia.............21 4.5 A discontinuity within the cratonic lithosphere .......................................................22 5. Conclusions ...............................................................................................................24 6. Acknowledgements ..................................................................................................25 References .....................................................................................................................26 Figures......................................................................................................................32 Table ..............................................................................................................................42 Chapter 2. The effects of anisotropy on Sp receiver functions ............. 43 Abstract.........................................................................................................................44 1. Introduction ..............................................................................................................46 2. The interaction of Sp receiver functions with azimuthally anisotropic structures ........................................................................................................................................48 2.1 Generating synthetic seismograms ..........................................................................48 2.2 Elastic coefficients ...................................................................................................48 2.3 Calculating synthetic receiver functions ..................................................................49 3. Synthetic Receiver Function Results ......................................................................50 3.1 Comparison to Ps receiver functions .......................................................................50 3.2 Variations in anisotropy geometry ...........................................................................52 3.3 Sp receiver functions in the presence of radial anisotropy ......................................53 3.4 Coupling between P, SV and SH .............................................................................54 ! xvi! 3.5 Strategies for addressing initial particle motion dependence ............................57 4. Data example ............................................................................................................58 4.1 Craton structure and formation ................................................................................58 4.2 Deconvolution and migration ..................................................................................60 4.3 Results ......................................................................................................................60 4.4 Comparison to previously modeled receiver functions ...........................................61 4.5 Comparison to tomography results ..........................................................................62 4.6 FCC SKSp receiver functions ..................................................................................63 4.7 Model parameterization ...........................................................................................64 4.8 Translating model parameters into elastic coefficients ............................................65 4.9 Computing and comparing synthetics SKSp receiver functions..............................68 4.10 Implications............................................................................................................69 5. Conclusions ...............................................................................................................69 References .....................................................................................................................71 Figures...........................................................................................................................75 Tables ............................................................................................................................85 Chapter 3. A comparison of the seismological and rheological lithosphere-asthenosphere boundary beneath cratonic North America 89 Abstract.........................................................................................................................90 1. Introduction ..............................................................................................................92 2. Shear wave velocity profiles ....................................................................................93 2.1 Shear-wave velocity profiles....................................................................................93 2.2 Attenuation effects ...................................................................................................96 3. Estimating cratonic geotherms ...............................................................................98 3.1 Velocity into temperature ........................................................................................98 3.2 Comparison to xenolith-constrained geotherms ....................................................100 3.3 Accounting for compositional variations ...............................................................102 4. Viscosity of the cratonic lithosphere ....................................................................107 4.1 Translating temperature into viscosity ...................................................................107 4.2 Viscosity profiles – Slave craton ...........................................................................109 4.3 Viscosity profiles – Superior craton ......................................................................110 5. Implications ............................................................................................................111 5.1 Lithospheric properties ..........................................................................................111 5.2 Asthenospheric properties ......................................................................................113 5.3 A boundary in viscosity .........................................................................................114 6. Conclusions .............................................................................................................114 References ...................................................................................................................116 Figures.........................................................................................................................119 Tables ..........................................................................................................................127 Chapter 4. Localized shear in the deep lithosphere beneath the San Andreas Fault system .................................................................................132 Abstract.......................................................................................................................133 Text ..............................................................................................................................134 Methods Summary .....................................................................................................141 ! xvii! Acknowledgements ....................................................................................................142 References ...................................................................................................................143 Figures.........................................................................................................................145 Appendix A. Supplementary materials for: Localized shear in the deep lithosphere beneath the San Andreas Fault system ................................149 1. Additional information on receiver function and CCP stacking methods .......150 1.1 Optimal receiver function filtering ........................................................................150 1.2 LAB phase selection ..............................................................................................151 1.3 Influence of ray parameter .....................................................................................154 2. Modeling the relative velocity drop across the plate boundary .........................155 3. Other variations in LAB phase depth and amplitude ........................................156 References ...................................................................................................................160 Figures.........................................................................................................................162 ! xviii! Chapter 1 The lithosphere-asthenosphere boundary and cratonic lithospheric layering beneath Australia from Sp wave imaging By: Heather A. Ford1, Karen M. Fischer1, David L. Abt1,2, Catherine A. Rychert3, and Linda T. Elkins-Taunton4 1 Department of Geological Sciences, Brown University 2 ExxonMobil Exploration Company, Houston, Texas 3 Scripps Institute of Oceanography, UC San Diego 4 Dept. of Earth, Atmospheric and Planetary Science, Massachusetts Institute of Technology Published in its present form in: Earth and Planetary Science Letters (300) 299-310, 2010 Copyright Elsevier Ltd., 2010 ! 1! Abstract Sp and Ps scattered wave receiver functions were calculated for nineteen stations across Australia and the island of Tasmania in order to image the lithosphere- asthenosphere boundary and layering within the lithosphere. Within Phanerozoic eastern Australia and the eastern margin of the South Australia Craton, prominent Sp phases from a negative velocity contrast were found at depths of 61±11 km to 131±9 km, consistent with the lithosphere-asthenosphere boundary depth range from surface wave tomography. These phases imply significant velocity drops over depth ranges of 30-40 km or less, and thus cannot be explained by a lithosphere-asthenosphere boundary that is controlled by temperature alone. Rather, they imply that the asthenosphere is hydrated with respect to a drier, depleted lithosphere or contains a small amount of partial melt. The shallowest Sp phases have the largest amplitudes and occur in regions with the most recent, voluminous volcanism, strengthening the link to partial melt at the base of the lithosphere. In contrast, no significant negative Sp phases were found at the base of the thick cratonic lithosphere at stations in central and western Australia, implying that the cratonic lithosphere-asthenosphere velocity gradient is distributed over more than 50-70 km in depth. This gradient may be purely thermal in origin, although gradational changes in composition or melt content cannot be ruled out. A negative Sp phase was observed at depths of 69±8 km to 85±14 km at stations in central and western Australia, indicating the presence of a drop in velocity internal to the lithosphere. This interface within the lithosphere may be a relic of cratonic mantle formation, or the result of alteration by melt and metasomatism. ! 2! 1. Introduction Although the lithosphere is commonly defined as the strong, thermally conducting outer shell of the earth, its physical and chemical properties, its formation and evolution, and the mechanism for it’s stability, are still poorly understood, particularly beneath continents. Insight on these issues may be gained by determining the depths of the lithosphere-asthenosphere boundary (LAB) and any discontinuities internal to the lithosphere, the seismic velocity gradients associated with these discontinuities, and their variation among tectonic provinces. In this study, we focus on the continental lithosphere of Australia. 1.1. Overview of Australian tectonics The continent of Australia can be divided into three principal tectonic regions (Figure 1), which include the predominantly Archean West Australia Craton, the largely Proterozoic North and South Australia Cratons in central Australia, and the Phanerozoic accreted terranes located in the eastern portion of the continent. The West Australia Craton is the product of accretion of two Archean cratons, the Pilbara and Yilgarn, during the paleo-Proterozoic [Barley et al. 1998; Betts et al., 2002; Cawood and Tyler, 2004]. The North Australia and South Australia Cratons, located in central Australia, are amalgamations of Archean and Proterozoic terranes that assembled during the Proterozoic [Betts et al., 2002; Giles et al., 2001; Tyler, 2001]. A long-lived accretionary margin along the southern edge of the North Australia Craton resulted in ! 3! collision with the West Australia Craton and with portions of the South Australia Craton during the paleo-Proterozoic [Betts et al., 2002]. During the Phanerozoic, accretion continued along the eastern margin of the Precambrian craton(s) in a series of subduction related events [Betts et al., 2002]. The Tasman Line, first defined by the work of Hill [1951] through outcrop mapping, was thought to be a boundary dividing Proterozoic craton from Phanerozoic basement. However, more research has led to a number of different proposed locations and interpretations of the Tasman Line and a recent review has concluded that the accretion of eastern Australia onto the craton cannot be defined by a single line [Direen and Crawford, 2003]. In addition to accretion, more recent Cenozoic volcanism (perhaps associated with mantle plumes [e.g. Wellman, 1983]) has impacted much of the eastern margin (Figure 1). Magmatic ages from central volcanoes follow a well-established decrease from north to south, whereas lava field volcanism follows no known trend in age [Johnson et al., 1989; Wellman and McDougall, 1974]. Over the last 10 My, the largest volumes of magmatism have occurred at the southern and northern reaches of eastern Australia (Figure 1), with eruptions as recently as 4600 years ago within the Newer Volcanics in the south and 13,000 years ago in the north [Johnson et al., 1989]. 1.2. Results from surface wave tomography At the continental scale, surface wave tomography has indicated dramatically thicker lithosphere beneath cratonic Australia (150-250 km) than beneath Phanerozoic eastern Australia (as thin as 50 km) [Debayle and Kennett, 2000; Fichtner et al., 2009, 2010; Fishwick et al. 2005, 2008; Simons and van der Hilst, 2002, 2003; Simons et al., ! 4! 1999]. However, more detailed estimates of lithospheric thickness and its correlation with crustal age differ between these studies. Although surface wave tomography is excellent for imaging volumetric heterogeneity, it cannot distinguish between differences in vertical mantle velocity gradient thicknesses that occur over depth ranges of 50 km or less. S-to-P (Sp) and P-to-S (Ps) scattered waves are a good alternative for studying the nature of the LAB and other mantle interfaces because of their ability to more precisely constrain boundary depths and vertical velocity gradients. The LAB has been inferred in a variety of tectonic settings from Ps receiver functions [e.g., Chen et al., 2006; Collins et al., 2002; Li et al., 2000; Ozacar et al., 2008; Rychert et al., 2005; Rychert and Shearer, 2009; Snyder, 2008; Wolbern et al., 2006] and Sp receiver functions [e.g., Abt et al., 2010; Hansen et al., 2007, 2009; Heit et al., 2007; Kumar et al., 2005a, 2005b, 2007; Li et al., 2004, 2007; Mohsen et al., 2006; Oreshin et al., 2002; Sacks et al., 1979; Sodoudi et al., 2006; Vinnik et al., 2005] or a combination of both [e.g. Kawakatsu et al., 2009; Rychert et al., 2007; Wittlinger and Farra, 2007]. Reviews may be found in Rychert et al. [2010] and Fischer et al. [2010]. In Australia, a number of studies using Ps receiver functions have imaged the seismic crust-mantle boundary (Moho) [e.g. Clitheroe et al., 2000a, 2000b; Reading and Kennett, 2003; Reading et al., 2003, 2007]. A Sp receiver function study by Kumar et al. [2007] observed a negative phase at four stations in Australia and interpreted it as the LAB at depths ranging from 90 km to 207 km. In contrast, a Ps study by Rychert and Shearer [2009] identified a negative phase at five stations at depths of 71-106 km and interpreted it as the LAB off the craton or a mid-lithospheric discontinuity within the craton. Our study more completely characterizes the LAB and mid-lithospheric discontinuities ! 5! throughout Australia with Sp receiver functions at that sample all three primary tectonic regions. 2. Method 2.1. Ps and Sp receiver functions The basic premise of the receiver function method is that the deconvolution of the parent phase (e.g., P for Ps; SV for Sp) from the daughter component (e.g., SV for Ps; P for Sp) removes source and instrument effects, enhancing information regarding seismic structure beneath the station of interest. Although Ps receiver function analysis has been successfully used for many years to image the Moho and relatively deep mantle discontinuities such as those in the transition zone, interpretation of phases in terms of shallow upper mantle structure requires great care due to the presence of crustal multiples which can contaminate direct arrivals; these effects are typically most pronounced at depths of less than 200-250 km [e.g. Bostock, 1998; Rychert et al., 2005, 2007]. In Sp receiver functions, reverberations associated with the crust arrive after the direct S phase, while scattered Sp phases arrive before the direct phase. This separation prevents contamination by crustal multiples, making Sp scattered waves a useful tool for imaging upper mantle discontinuities such as the lithosphere-asthenosphere boundary. Limitations of Sp phases include possible interference with SKS and SKSp at distances greater than 75° and with P-wave phases, such as pPPP, pPPPP, and sPPPP, for events at depths greater than 300 km [Wilson et al., 2006], and the fact that they are not produced at mantle depths for distances of less than roughly 55˚ [Yuan et al., 2006]. In addition, Ps waves are often better at resolving the depth extent of vertical velocity gradients that are ! 6! distributed over less than 20-30 km, due to their typically higher frequency energy content. Moreover, because the paths of incoming Ps phases are closer to vertical than Sp paths for a given source-station distance, they sample a smaller region around the station, reducing the potential for averaging of laterally-varying discontinuity structure. As a result of the above considerations both Ps and Sp were analyzed, with distance limits of 35°-80° for Ps and 55°-75° for Sp (Figure 2), and a depth limit of less than 300 km for Sp. Receiver functions were binned by epicentral distance in order to differentiate phases of interest from crustal multiples in Ps and unwanted teleseismic phases in Sp. 2.2. Data preparation Receiver functions were calculated for nineteen permanent broadband stations located in Australia and on the island of Tasmania. An event list was compiled from the USGS National Earthquake Information Center (NEIC) global event catalog and included events with an Mw>=5.8, an epicentral distance of 35°-80° and no restrictions in depth, occurring through March 2009. Waveform data was obtained from the Incorporated Research Institutions for Seismology (IRIS) Data Management Center (DMC) and included data from networks AU, G, II and IU. The number of waveforms used varies between stations due to differences in quality as well as station operation time (Table 1). A free surface transfer matrix [Bostock, 1998; Kennett, 1991] was used to transform the recorded waveforms into the P and SV components that would have been incident on the free surface, assuming a given ray parameter and surface Vp and Vs. To optimize the free surface transformation, the surface Vp and Vs for each station were determined using ! 7! an automated procedure that is detailed in Abt et al. [2010]. Essentially, the parent phase on each waveform (P for Ps or SV for Sp) was windowed around its arrival time, and a search was performed over a range of Vp and Vp/Vs to find the value(s) that minimized the correlation of the parent phase with its corresponding window on the daughter component (SV for Ps and P for Sp). For each station, all cross-correlation surfaces with well-defined minima were stacked, and the best free-surface velocities for the station were defined as the minimum of this stack. During this process, checks were put into place to ensure that only quality waveforms were included in the analysis. These measures included rejecting waveforms with signal-to-noise ratios below a selected value (5 for Ps and 2 for Sp) and whose observed arrival-times differed from the predictions of the AK135 velocity model [Kennett et al., 1995] by more than 5 seconds for Ps or 10 seconds for Sp. 2.3. Deconvolution and migration We employed two different deconvolution techniques in order to ensure that phases observed on the receiver functions were not artifacts of a given deconvolution approach. In the first method, all waveforms for a given phase (Sp or Ps) and station were simultaneously deconvolved and migrated in the frequency domain, using a best- fitting regularization parameter (i.e. water-level) to stabilize the deconvolution [Bostock, 1998]. In the second approach, iterative time domain deconvolution [Kikuchi and Kanamori, 1982; Ligorría and Ammon, 1999] was applied to individual Sp or Ps and waveforms; waveforms for a given phase (Sp or Ps) and station were then migrated and stacked. Prior to deconvolution the waveforms are bandpass filtered (0.03-1 for Ps and ! 8! 0.03-0.5 for Sp). The 1-D velocity models used for migration vary between stations. Crustal thickness and Vp and Vs values were obtained from H-k stacking [Zhu and Kanamori, 2000] of Ps waveforms. H-k stacking results were considered robust if the estimated Moho depth fell within the error bars of the Moho phase selected from the Ps receiver function. Complicated crustal structure at FORT, where H-k stacking was unable to constrain a single best pick, required a fixed crustal model constructed using estimates of Moho depth from Clitheroe et al. [2000]. AK135 [Kennett et al., 1995] was assumed for the mantle at all stations. This choice of 1D mantle velocity model does not account for possible variations in mantle Vp/Vs and may produce systematic errors in phase depth estimation. We estimate that the resulting uncertainty in inferred mantle discontinuity depths is on the order of 10 km or less. For example, Rychert et al. [2007] demonstrated that changing Vp/Vs from 1.7 to 1.8 could vary the depth of a phase located at approximately 90 km in a Sp receiver function by roughly 6-8 km. Abt et al. [2010] migrated Sp and Ps data using two different mantle velocity models, AK135 and 1D profiles extracted from 3D Vs [Yuan and Romanowicz, 2010] and Vp [Burdick et al., 2008] models for North America, and found that differences in inferred discontinuity depths in the 50 to 115 km range were 6 km or less. Kaiho and Kennett [2000] observed as much as a 6% Vp/Vs decrease with respect to AK135 in parts of northern Australia at depths of 35-120 km. If a 6% Vp/Vs drop is assumed, Sp phases that would appear at 130 km in an AK135 mantle would instead be located at 140 km. A discontinuity at a depth of ~75 km would only be perturbed by 5 km. For a given depth in the migration, a waveform was included only if a direct, pre- critical phase was predicted to exist for the direct phase ray parameter, assuming the 1D ! 9! velocity model for the given station. This step reduces the number of events at a given depth in addition to eliminating deeper portions of some Sp receiver functions. For example, at a depth of 380 km only events with a ray parameter of 0.105 (70° epicentral distance) or smaller are included, whereas at a depth of 150 km, the largest ray parameter that can be included is 0.120 (55° epicentral distance). Ps phases do not reach post- critical incidence over the range of epicentral distances and depths included in this study. In order to understand the uncertainties associated with each receiver function, error bars were calculated using a bootstrap test in which 20% of the waveforms were randomly replaced with another random 20%. The replacement and recalculation was repeated 100 times in order to determine the mean and standard deviation. 3. Results Ps and Sp receiver functions were calculated for the nineteen stations shown in Figure 1 using the frequency- and time-domain deconvolution techniques. At fourteen stations, these methods produced consistent Sp and Ps receiver functions. For example, at station ARMA (Figure 3) the depth to the large, positive (Moho) phase is very similar on the Sp receiver functions generated by the different deconvolution methods. Moreover, the largest negative phase present on the frequency-domain receiver function is found at a depth of 93±16 km while the largest negative phase on the time-domain receiver function is found at a depth of 81±18 km; each estimate is within the error bars of the other. Receiver functions that varied greatly between the time- and frequency- domain (e.g. negative phases not within error of each other) were discarded, whereas stations where the two methods yielded receiver functions with comparable shape and ! 10! phase depths were included in the interpretation. Moho and upper mantle discontinuity depths measured from the frequency-domain receiver functions are listed in Table 1. For each station, Moho depth was found from the largest positive phase on the single-binned Ps receiver function. Measured Moho depths were in general consistent between Ps and Sp receiver functions and agree with crustal thickness estimates from H-k stacking where the latter exist (Table 1). In Phanerozoic Australia, Moho depths range from 32 ± 2 km at TAU to 40 ± 2 km at CTAO. On the craton, Moho depths vary from 32 ± 2 km at MBWA to 53 ± 2 km at MUN. Phases associated with upper mantle discontinuities were interpreted from Sp receiver functions. Most of the Ps receiver functions also contain negative energy at depths similar to the negative phases seen in Sp receiver functions. However, strong reverberations in Ps receiver functions, including apparent reverberations from intracrustral layering that exhibit negligible moveout with epicentral distance, complicate the interpretation of possible upper mantle discontinuities, leading us to report and interpret mantle phases only from Sp receiver functions. At station TOO both Ps and Sp receiver functions have strong, well-resolved negative phases (Figure 4a, b) found within error of each other. It is possible that the negative phase in Ps is the result of an upper mantle discontinuity; however, influence from a midcrustal discontinuity cannot be ruled out. At station FORT the presence of shallow crustal layering produces reverberations that could potentially interfere with upper mantle phases within the Ps receiver function (Figure 4c). The strong negative arrival observed at 79 ± 6 km on the Sp receiver function for FORT does not suffer from reverberation contamination (Figure 4d). ! 11! At mantle depths, Sp receiver functions are in general dominated by a single negative phase over the depth range imaged (300 km). We interpret this phase and report its depth in Table 1 when it is reasonably consistent across epicentral distance bins at a given station (e.g. Figure 5) and between frequency- and time-domain decovolution results. Although smaller arrivals also appear in the Sp receiver functions, they are less consistent between stations and often also less consistent as a function of epicentral distance. In order to determine whether our negative phase pick from the Sp receiver function is the LAB or some other upper mantle discontinuity we compared the depth of our pick for a given station to the corresponding absolute shear wave velocity profile (from 75 to 300 kilometers) constructed from the updated shear wave velocity model of Fishwick et al. [2008a, 2008b]. For the negative phase pick of the Sp receiver function to be interpreted as the LAB, it must fall within the depth range between the minimum Vs beneath the Moho (which presumably lies within the asthenosphere) and the next maximum in Vs in the upward direction (which presumably lies within the lithosphere). At three stations in eastern Australia (COEN, CTAO, and TOO) the minimum Vs lies at 75 km, which corresponds to the top of the Fishwick et al. [2008a, 2008b] model; in these cases the potential LAB depth range is assumed to be from the Moho to 75 km. If the depth of the negative phase is above the LAB depth range, then it is considered to be a mid-lithospheric discontinuity (MLD). Station BBOO is the one exception, as discussed in section 3.2. 3.1. Phanerozoic Australia Of the nine stations within Phanerozoic Australia, six were found to have consistent ! 12! Sp receiver functions from time- and frequency-domain deconvolution methods (e.g. station ARMA, Figure 3). In regard to the other three, at MOO and TAU constraints on Moho depth were obtained from Ps receiver functions and H-k stacking (Table 1); at CAN, even crustal structure was not well constrained. Station COEN (Figure 5) is representative of other stations located along the eastern margin of the continent. The Sp receiver function contains a strong negative phase at 67±8 km that lies within the potential LAB depth range and is interpreted as the lithosphere-asthenosphere boundary. Sp receiver functions for the six Phanerozoic stations show well-resolved negative phases that range in depth from 61±11 km at TOO to 93±16 km at ARMA (Figure 6b and Table 1). All of these phases fall within the potential LAB depth range defined from the shear wave model of Fishwick et al. [2008a, 2008b] and are interpreted as the LAB. 3.2. Cratonic Australia Within the Proterozoic and Archean portions of the Australian continent, Ps and Sp receiver functions were calculated for ten stations. Of these stations eight have interpretable mantle Sp phases, while at BLDU and MUN only crustral structure was obtained. Figure 6 (c and d) shows two profiles of the Sp receiver functions in southern and northern Australia that intersect Archean, Proterozoic and Phanerozoic terranes. Stations STKA and BBOO are located near the eastern edge of cratonic Australia (Figure 1). At STKA, the largest significant negative phase in Sp is located at a depth of 104±9 km and is consistent between bins of epicentral distance (Figure 5). The negative phase falls within, but near the edge of, the potential LAB depth range for the absolute Vs profile beneath STKA and is interpreted to be the LAB (Figure 6c). To the west of ! 13! STKA, station BBOO’s deepest, most statistically significant negative phase is found at 131±9 km (Figure 6c). This depth falls outside of the potential LAB depth range estimated from the Vs profile and the phase could be interpreted to represent a mid- lithospheric discontinuity (MLD). However, the negative phase lies only 10 km outside of potential LAB depth range. In addition, its depth is considerably greater than the MLD phases observed at other cratonic stations. Thus we interpret the negative phase to be the LAB, while acknowledging the ambiguity of this choice. If the negative phases at STKA and BBOO are interpreted to be LAB, then an argument could be made for a dipping boundary that increases in depth to the west (Figure 6c). At the remaining stations within cratonic Australia, a prominent negative Sp phase is observed within the cratonic lithosphere (a MLD), but no clear arrival is observed from depths comparable to the LAB (Figure 6c and 6d). Station WRAB is located well within the interior of the North Australia Craton. A large negative phase at a depth of 81±14 km appears in the single-binned receiver function and is also consistent between bins of epicentral distance (Figure 5). The inferred LAB depth range at WRAB from Fishwick et al. [2008a, 2008b] is 175-200 km, which is significantly deeper than the largest negative phase on the Sp receiver function, leading to the latter’s interpretation as a MLD. No significant negative phase occurs in the LAB depth range at WRAB. Elsewhere in the craton along the northern profile (Figure 6d), the largest negative phases at FITZ and MBWA are similar in depth to the negative phase at WRAB. They also lie above the potential LAB depth range and each is inferred to be a MLD. Along the southern margin of the craton at stations FORT, KMBL, and NWAO (Figure 6c) the negative Sp phase depths are too shallow to be the LAB when compared to the respective ! 14! shear wave velocity profile and are each taken to be a MLD. No significant negative phase occurs within the LAB depth range at these stations. It should also be noted that although the receiver functions are displayed to depths of 250 km, both Ps and Sp receiver functions were calculated to 400 km for all stations and were examined to determine whether significant, well-constrained phases existed at greater depths. In Sp receiver functions, depths greater than approximately 250-300 km are often noisy and poorly constrained, with large differences in time and frequency- domain deconvolution methods. This is due to the small number of waveforms included as a result of our removing waveforms from depths for which no direct, pre-critical phase was predicted to exist. Although we cannot rule out the possibility of a significant negative Sp phase at depths of 250-400 km, such an arrival would lie well below the potential LAB depth range from the Fishwick et al. [2008a, 2008b] velocity model. 3.3. Large-scale correlation in depth and amplitude estimates In order to better visualize how the LAB and MLD vary between stations, two plots of smoothed Sp conversion points, color-coded as a function of depth and amplitude, are shown in Figure 7 (see figure caption for smoothing details). From this point forward, stations whose negative phase is interpreted to be the LAB will be referred to as stations of eastern Australia. It should be noted that these stations sample mantle beneath both cratonic (STKA and BBOO) and Phanerozoic (COEN, CTAO, EIDS, ARMA, YNG, TOO) regions. Stations that are referred to as being located in central and western Australia are cratonic stations and have a negative phase interpreted to be a MLD. ! 15! In eastern Australia, a large variation occurs in the depth and amplitude of the LAB phase (Figure 7), and depth and amplitude are negatively correlated. Where the LAB is found to be relatively shallow, its amplitude is larger (Figure 8). Intriguingly, the locations of the shallowest, strongest amplitude, negative phases fall in or near regions of the most recent, voluminous volcanism (see Figure 1). The relationship between these observations is further discussed in section 4.3. In contrast, the depths and amplitudes of the MLD phases observed in central and western Australia are more tightly clustered (Figure 8) and no observable variation exists between cratonic blocks. 4. Discussion 4.1. Comparison to other receiver function studies Previous receiver function studies in Australia have used Ps scattered waves to determine crustal structure [Clitheroe et al., 2000; Reading and Kennett, 2003; Reading et al., 2003, 2007]. A continent-wide study of crustal structure by Clitheroe et al. [2000] found crustal thickness values at/near several stations in our study, including CTAO, ARMA, CAN, KMBL, NWAO, and WRAB [see figure 8, Clitheroe et al., 2000]. At these locations the Moho depth from our results and the results of Clitheroe et al. [2000] are within error of each other. More recent work by Reading et al. [2003] found that the crustal thickness of the Yilgarn Craton varies slightly from west to east, but is approximately 35 to 40 km, which agrees well with our results at stations NWAO and KMBL. Reading and Kennett [2003] found that Moho depth is 30±2 km in the Pilbara Craton, in agreement with our measured Moho depth of 32±2 km at station MBWA. ! 16! Turning to the mantle, a Sp receiver function study by Kumar et al. [2007] interpreted the depth of the LAB at stations NWAO, WRAB, and STKA to be 164 km, 180 km, and 207 km respectively. In contrast, we did not find any significant negative phase at these depths at NWAO and WRAB. At STKA, a small amount of negative energy does appear around 200 km, but its amplitude is much smaller than the phase we observe at 104 ± 9 km. Our results are more consistent with the conclusions of Rychert and Shearer [2009] based on Ps receiver functions. Rychert and Shearer [2009] found discontinuity depths of 71 to 106 km for three of the same stations as our study, both on and off the craton. At NWAO they observe a negative phase at 71 kilometers, which comes close to the negative phase at 81±8 kilometers observed from our Sp receiver function. At stations WRAB and CTAO, Rychert and Shearer [2009] identified negative phases at 106 and 86 km, respectively, which are 11 and 9 km outside of the error bars for the Sp arrivals we observe at those stations. 4.2. Estimating the gradient thickness and velocity contrast at a boundary To better understand the mechanism(s) that may responsible for producing a boundary in seismic velocity, it is important to constrain the velocity gradient parameters that describe the boundary. In the case of the negative Sp phases considered in this study, the parameters are the velocity drop and the gradient thickness (the depth range over which the velocity drop is distributed). To more completely constrain the gradient parameters, detailed inverse modeling is needed [e.g., Rychert et al., 2005; Rychert et al., 2007]. However, simple forward modeling can still be useful in constraining the gradient and in differentiating between possible mechanisms. ! 17! For each modeled receiver function, velocity gradient parameters were systematically varied to determine the parameter ranges that provide acceptable fits to the observed LAB or MLD. The tested velocity models contained the specific crustal structure for the given station and an LAB or MLD described by velocity drop and gradient thickness values. Gradient thicknesses of 0 km to 90 km were tested at 10 km increments. Velocity decreases of up to 10% were considered; 10% is greater than the total shear velocity drop from lithosphere to asthenosphere typically seen in surface wave models [Gaherty et al., 1999; Nettles and Dziewonski, 2008; Romanowicz, 2009; Yuan and Romanowicz, 2010]. A propagator matrix method [Keith and Crampin, 1977] was used to generate synthetic waveforms that were turned into receiver functions using the same processes that were applied to the observed seismograms. In the data, the dominant period of the incident S waveform was 10-11 seconds, and a source-time function with a dominant period of 10.5 seconds was employed in the synthetics. In addition, the number and distribution of ray parameters in the synthetic seismograms duplicated the distribution of events in the observed receiver functions. The results of the modeling are discussed in sections 4.3, 4.4. and 4.5. 4.3. The lithosphere-asthenosphere boundary in eastern Australia The negative Sp phase observed at depths of approximately 61±11 km to 131±9 km beneath stations in eastern Australia is interpreted to be the LAB based on comparison with shear wave velocity structure [Fishwick et al., 2008a, 2008b]. To constrain the velocity gradient associated with the LAB, the LAB Sp phases were modeled for station COEN which has the largest amplitude LAB phase, ARMA which ! 18! has the smallest amplitude LAB phase on the eastern margin of Australia, and BBOO which has the smallest amplitude LAB phase overall. For station COEN, the best fitting velocity gradients vary from a 7% velocity decrease over 0 km, to a 10% decrease in velocity over 20 km (Figure 9), but all models with gradient thicknesses of 40 km or more or velocity drops of 4% or less failed to match the observed receiver function to within two standard deviations. At station ARMA, the best fitting gradients vary from a 4% velocity drop over 0 km, to a 7% velocity drop over 10 km and a 10% drop in velocity over 20-30 km (Figure 9). At ARMA and BBOO, a gradient thickness of 40 km for a 10% velocity drop grazes the lower amplitude two standard deviation limit of the observed Sp phases. Overall, for stations representing the range of LAB Sp phase amplitudes in eastern Australia, gradient thicknesses of 40 km or less are required, and smaller gradient thicknesses provide better fits. In geodynamical models for cratonic lithosphere and surrounding continental margins [Cooper et al., 2004; King and Ritsema, 2000; Korenaga and Jordan, 2002] temperature gradients between the lithosphere and asthenosphere occur over at least 50- 70 km. In contrast, the modeling of eastern Australia LAB Sp phases in this study rules out velocity gradient thicknesses of more than 40 km. We thus conclude that the LAB in eastern Australia cannot be the result of a change in temperature alone. Another possible explanation for the velocity drop is that the lithosphere is more dehydrated and depleted relative to a hydrated and fertile asthenosphere [Hirth and Kohlstedt, 1996; Karato and Jung, 1998]. Mg numbers for the mantle of the Phanerozoic Australia lithosphere lie in the range of 90-91 [Gaul et al., 2000], relatively close to ! 19! values expected for the asthenosphere (88-89), with the result that depletion alone could reduce velocities by less than 1% [Lee, 2003]. However, hydration in the asthenosphere could create a drop in velocity of roughly 4.5% across the LAB [Rychert and Shearer, 2009], and thus could explain the eastern Australia LAB Sp phases, with or without depletion effects, assuming that the velocity gradient occurs over 10 km or less. Alternatively, a small fraction of melt in the asthenosphere could produce a large drop in velocity. The exact percent melt needed to produce a given percent change in seismic wave speed depends on melt geometry [e.g., Hammond and Humphreys, 2000; Takei, 2002; Takei and Holtzman, 2009], but 1-2% appears to be sufficient [Hammond and Humphreys, 2000; Kawakatsu et al., 2009]. The xenolith-based southeast Australia (SEA) geotherm [O’Reilly et al., 1997] is thought to reflect the present-day geotherm in high heat flow areas like north Queensland (e.g., COEN), east-central Queensland (e.g., CTAO) and western Victoria (e.g., TOO) [O’Reilly et al., 1997]. O’Reilly et al. [1997] inferred that the SEA geotherm crossed the dry peridotite solidus at a depth of approximately 120±20 km, significantly deeper than the negative velocity gradient indicated by the Sp phases. However, reduction of the asthenospheric solidus temperature due to hydration [e.g. Grove, 2006; Hirschman et al., 2009] creates the possibility that the eastern Australia LAB coincides with the damp solidus and that the LAB velocity gradient reflects partial melt within the asthenosphere. Alternatively, small degrees of carbonatite melt may exist in the asthenosphere [Dasgupta and Hirschmann, 2007]. A link appears to exist between eastern Australia LAB properties, lithospheric temperature, and present-day partial melt. The highest heat flow areas [O’Reilly et al., ! 20! 1997] coincide with regions of recent and voluminous volcanism (Figure 1) and with stations COEN, CTAO, and TOO, where the lithosphere appears to be thinner and the amplitude of the LAB phase is larger than at other eastern Australia stations (Figure 7). An interesting question is whether the zones of thinner lithosphere created [e.g., Ebinger and Sleep, 1998], or were created by, focused mantle flow and melting. 4.4. The lithosphere-asthenosphere boundary in central and western Australia In the craton, surface wave tomography models indicate that the transition from seismically fast lithosphere to slow mantle asthenosphere occurs at depths of approximately 150-250 km in Australia [Fishwick et al., 2008a, 2008b]. The absence of significant negative energy in Sp receiver functions at LAB depths is a striking feature of stations in central and western Australia. At cratonic stations near the coast, it is possible that the negative phase found at depths of 69±8 km to 85±14 km, and interpreted to be a MLD, is actually the LAB and that the edge of the thick lithosphere is further inland than resolved by surface wave tomography. However, this explanation is highly unlikely at stations in the continental interior, such as WRAB. It is also possible that the lack of a negative phase at potential LAB depths could reflect the lack of an asthenospheric layer that contains isotropically slow velocities [Gaherty and Jordan, 1995; Pedersen et al., 2009; Revenaugh and Jordan, 1991]. Assuming that an asthenospheric low velocity zone does exist beneath the Australian craton [Cammarano and Romanowicz, 2007; Fishwick, 2005; Romanowicz, 2009], modeling with synthetic Sp receiver functions for periods of 10.5 s, comparable to dominant periods in the observed receiver functions, indicates that increasing gradient ! 21! thickness from 0 km to 50 km produces a 70% reduction in the amplitude of the Sp phase produced by the discontinuity. For velocity drops of no more than 10%, gradient thicknesses of 60 km or more could easily prevent the observation of Sp phases when typical noise levels are taken into account. Sp receiver functions for stations WRAB, FITZ, MBWA and NWAO were re-calculated decreasing the low-pass filter that increased the average dominant period of the individual waveforms from ~10.5 seconds to ~20 seconds. Even with these larger dominant periods, no Sp phases from LAB depths were observed. We conclude that if a significant reduction in velocity occurs at the base of the cratonic lithosphere, it must be distributed over 60-90 km or more. Such velocity gradients could be produced by geotherms typical of models where no change in composition or melt content occurs, although comparably gradual vertical variations in composition or melt cannot be ruled out (Figure 10). 4.5. A discontinuity within the cratonic lithosphere The negative Sp phase at 69±8 km to 85±14 km imaged throughout central and western Australia is interpreted to be a negative velocity contrast internal to the lithosphere (a MLD). This feature correlates with a layer of low shear wave velocity in surface wave models, seen in central Australia at 75 km depth by Fishwick et al. [2005] and more broadly in cratonic Australia by Lekic and Romanowicz [submitted]. A decrease in velocity within the cratonic lithosphere at similar depths has also been observed in receiver functions from North America [Abt et al., 2010] and globally in Ps receiver functions [Rychert and Shearer, 2009], surface wave tomography [Romanowicz, 2009] and long-range seismic profiles [Thybo, 2006 and references therein]. ! 22! Discontinuities in the 70-100 km depth range in the continental lithosphere have sometimes been associated with the Hales discontinuity. However, the original definition of this discontinuity was a velocity increase with depth [Hales, 1969], and subsequent citings have included both positive discontinuities [Revenaugh and Jordan, 1991] and anisotropic boundaries [Bostock, 1998; Fuchs, 1983; Levin and Park, 2000; Mercier et al., 2008]. To test whether the negative isotropic velocity gradients associated with MLDs represent averages of an azimuthally anisotropic boundary, the Sp receiver functions were binned in 60° increments as a function of back-azimuth. No consistent back- azimuthal patterns in timing or amplitude were observed for the MLD phases (or for LAB Sp phases at eastern Australia stations), although a lack of events in a number of back- azimuth bins hindered the analysis. We conclude that anisotropic layering is not evident, but it cannot be ruled out. Binning by back azimuth was also used to look for variations in Sp phase depth at stations near regions of apparent rapid transitions in LAB depth (e.g. COEN, STKA, and BBOO) inferred from surface wave tomography [Fishwick et al., 2008]. However, robust trends in Sp phase depth internal to the data for a single station were not resolved. Using synthetic Sp receiver function modeling, we determined that a single negative phase similar in amplitude and depth to the observed MLD phases can be produced by isotropic models for the cratonic lithosphere that contain either a thin low velocity layer (< 6 km in vertical extent) or a localized drop in velocity (for example a 0 km thick gradient) followed by a gradual increase in velocity of equal magnitude (>50 km thick gradient). Given that aspects of these structures are very localized in depth, ! 23! they are likely related to factors such as composition or grain size and fabric, rather than vertical variations in temperature. Partial melt could produce sufficiently sharp vertical boundaries, but temperatures estimated for the cratonic lithosphere in Australia [O’Reilly et al., 1997] lie beneath the peridotite solidus, even allowing for the presence of water [Grove et al., 2006]. Layering in composition or texture could date to the formation of the cratonic lithosphere, perhaps related to imbrication of originally thinner lithosphere during mantle accretion. Stacking of lithospheric layers has been suggested in the Canadian shield [Bostock, 1998; Chen et al., 2009; Mercier et al., 2008; Snyder, 2008], although the discontinuity dips apparent in these studies are not obviously consistent with the relatively uniform depth of the MLD observed in Australia. Another possibility is that the MLD represents the top of a melt cumulate layer (such as a low-velocity pyroxenite [Behn and Kelemen, 2006]) emplaced in the cratonic lithosphere during an earlier time when lithospheric temperatures were higher. 5. Conclusions A strong, coherent negative Sp phase at 61±11 km to 131±9 km in eastern Australia is interpretable as the lithosphere-asthenosphere boundary (LAB). The drop in velocity required to produce the observed phases is too localized in depth (< 40 km) to be produced by models in which seismic velocities depend solely on temperature. Rather, the asthenosphere must be made weak relative to the lithosphere by other properties, for example greater water content or a small amount of partial melt. The strongest and shallowest LAB Sp phases correlate with regions of high heat flow and the most recent ! 24! large volume magmatic eruptions in Australia, suggesting a link between LAB topography and mantle melting processes. The absence of a negative Sp phase associated with the LAB in cratonic Australia implies that the velocity drop associated with the cratonic LAB is very gradual (distributed over 60-90 km or more) or very weak. Such a gradual boundary could be produced by temperature alone, although gradual variations in mantle composition or melt content cannot be ruled out. The only significant negative Sp phase for stations in central and western Australia is a mid-lithospheric discontinuity found at depths of 69±8 km to 85±14 km. This boundary could represent vertical variations in mantle composition, grain size or fabric, for example a low velocity melt cumulate layer. 6. Acknowledgements We thank Scott French for contributions to the analysis and modeling codes, Stewart Fishwick for the Australia shear wave velocity model, Greg Hirth for conversations about interpretation, and Rainer Kind for discussion regarding Sp receiver functions. Data were obtained from the IRIS Data Management System. This work was funded by NSF Geophysics award EAR-0538155. ! 25! References Abt, D.L., Fischer, K.M., French, S.W., Ford, H.A., Yuan, H., Romanowicz, B., 2010. North American lithospheric discontinuity structure imaged by Ps and Sp receiver functions. J. Geophys. Res., in press. Barley, M.E., Loader, S.E., McNaughton, N.J., 1998. Calc-alkaline volcanism in the McPhee Dome and Kelly Belt, and growth of the eastern Pilbara Craton. Precambr. Res., 88, 3-23. Behn, M.D., Kelemen, P.B., 2006. Stability of arc lower crust: Insights from the Talkeetna arc section, south central Alaska, and the seismic structure of modern arcs. J. Geophys. Res., 111, B11207, doi:10.1029/2006JB004327. Betts, P.G., Giles, D., Lister, G.S., Frick, L.R., 2002. Evolution of the Australian lithosphere, Aust. J. Earth Sci., 49, 661-695. Bostock, M.G., 1998. Mantle stratigraphy and evolution of the Slave province. J. Geophys. Res., 103, 21183-21200. Burdick, S., Li, C., Martynov, V., Cox, T., Eakins, J., Mulder, T., Astiz, L., Vernon, F.L., Pavlis, G.L., van der Hilst, R.D., 2008. Upper mantle heterogeneity beneath North America from travel time tomography with global and USArray Transportable Array data. Seis. Res. Lett., 79, 384-392. Cammarano, F., Romanowicz, B., 2007. Insights into the nature of the transition zone from physically constrained inversion of long-period seismic data. Proc. Natl. Acad. Sci., 104, 9139-9144. Cawood, P.A., Tyler, I.M., 2004. Assembling and reactivating the Proterozoic Capricorn Orogen: lithotectonic elements, orogenies, and significance. Precambr. Res., 128, 201-218. Chen, L., 2009. Lithospheric structure variations between the eastern and central North China Craton from S- and P-receiver function migration. Earth. Planet. Sci. Lett., 173, 216-277. Chen, L., Zheng, T.Y., Xu, W.W., 2006. A thinned lithospheric image of the Tanlu Fault Zone, eastern China: Constructed from wave equation based receiver function migration. J. Geophys. Res., 111, B09312. doi:10.1029/2005JB003974. Clitheroe, G., Gudmundsson, O., Kennett, B.L.N., 2000a. The crustal thickness of Australia. J. Geophys. Res., 105, 13697-13713. Clitheroe, G., Gudmundsson, O., Kennett, B.L.N., 2000b. Sedimentary and upper crustal structure of Australia from receiver functions. Aust. J. Earth Sci., 47, 209-216. Collins, J.A., Vernon, F.L, Orcutt, J.A., Stephen, R.A., 2002. Upper mantle structure beneath the Hawaiian swell: Constraints from the ocean seismic network pilot experiment. Geophys. Res. Lett., 29, 1522, doi:10.1029/2001GL013302. Cooper, C.M., Lenardic, A., Moresi, L., 2004. The thermal structure of stable continental lithosphere within a dynamic mantle. Earth. Planet. Sci. Lett., 222, 807-17. Dasgupta, R., Hirschmann, M.M., 2007. Effect of variable carbonate concentration on the solidus of mantle peridotite. Am. Miner., 92, 370-379. Debayle, E., and Kennett, B.L.N., 2000. Anisotropy in the Australasian upper mantle from Love and Rayleigh waveform inversion. Earth Planet. Sci. Lett., 184, 339-351. Direen, N.G., Crawford, A.J., 2003. The Tasman Line: Where is it, and is it Australia’s Rodinian breakup boundary?. Aust. J. Earth Sci., 50, 491-502. ! 26! Ebinger, C.J., Sleep, N.H., 1998. Cenozoic magmatism throughout east Africa resulting from impact of a single plume. Nature, 395, 788-791. Fischer, K.M., Ford, H.A., Abt, D.L., Rychert, C.A., 2010. The lithosphere- asthenosphere boundary. Ann. Rev. Earth Planet. Sci., 38, 551-575. Fishwick, S., Heintz, M., Kennett, B.L.N., Reading, A.M., Yoshizawa, K., 2008. Steps in lithospheric thickness within eastern Australia, evidence from surface wave tomography. Tectonics, 27, TC4009. doi:10.1029/2007TC002116. Fishwick, S., Reading, A.M., 2008. Anomalous lithosphere beneath the Proterozoic of western and central Australia: a record of continental collision and intraplate deformation?. Precambrian Res., 166, 111-121. Fishwick, S., Kennett, B.L.N., and Reading, A.M., 2005. Contrasts in lithospheric structure within the Australian craton--insights from surface wave tomography. Earth Planet. Sci. Lett., 231, 163-176. Fichtner, A., Kennett, B.L.N., Igel, H., Bunge, H.-P., 2009. Full seismic waveform tomography for upper-mantle structure in the Australasian region using adjoint methods. Geophys. J. Int., 179, 1703–1725. Fichtner, A., Kennett, B.L.N., Igel, H., Bunge, H.-P., 2010. Full waveform tomography for radially anisotropic structure: New insights into present and past states of the Australasian upper mantle. Earth Planet. Sci. Lett., 290, 270–280. Fuchs, K. 1983. Recently formed elastic anisotropy and petrological models for the continental subcrustal lithosphere in southern Germany. Phys. Earth Planet. Int., 31, 93–118. Gaherty, J.B., Jordan, T.H., 1995. Lehmann discontinuity as the base of an anisotropic layer beneath continents. Science, 26, 1468-1471. Gaherty, J.B., Kato, M., Jordan, T.H., 1999. Seismological structure of the upper mantle: a regional comparison of seismic layering. Phys. Earth Planet. Int., 110, 21–41. Gaul, O.F., Griffin, W.L., O’Reilly, S.Y., Pearson, N.J., 2000. Mapping olivine composition in the lithospheric mantle. Earth Planet. Sci. Lett., 182, 223-235. Giles, D., Betts, P.G., Lister, G.S., 2001. A continental backarc setting for Early to Middle Proterozoic basins of north-eastern Australia. Geol. Soc. Aust. Abstr. 64, 55- 56. Grove, T.L., Chatterjee, N., Parman, S.W., Médard, E., 2006. The influence of H2O on mantle wedge melting. Earth Planet. Sci. Lett., 249, 74-89. Hales, A.L., 1969. A seismic discontinuity in the lithosphere. Earth Planet. Sci. Lett., 7, 44–46. Hammond, WC, Humphreys ED. 2000. Upper mantle seismic wave velocity: Effects of realistic partial melt geometries. J. Geophys. Res., 105, 10,975-10,986. Hansen, S.E., Nyblade, A.A., Julià, J., Dirks, P.H.G.M., Durrheim, R.J., 2009. Upper- mantle low-velocity zone structure beneath the Kaapval craton from S-wave receiver functions. Geophys. J. Int., 178, 1021-1027. Hansen, S.E., Rodgers, A.J., Schwartz, S.Y., Al-Amri, A.M.S., 2007. Imaging ruptured lithosphere beneath the Red Sea and Arabian Peninsula. Earth Planet. Sci. Lett., 259, 256-265. Heit, B., Sodoudi, F., Yuan, X., Bianchi, M., Kind, R., 2007. An S receiver function analysis of the lithospheric structure in South America. Geophys. Res. Lett., 34, L14307. doi:10.1029/2007GL030317. ! 27! Hill, D., 1951. Geology, in: Mack, G. (Eds.), Handbook of Queensland, Australian Association for the Advancement of Science, Brisbane. Hirschmann, M.M., Tenner, T., Aubaud, C., Withers, A.C., 2009. Dehydration melting of nominally anhydrous mantle: The primacy of partitioning. Phys. Earth Planet. Int., 176, 54-68. Hirth, G., Kohlstedt, D.L., 1996. Water in the oceanic upper mantle: implications for rheology, melt extraction and the evolution of the lithosphere. Earth Planet. Sci. Lett., 144, 93- 108. Johnson, R.W., Knutson, J. & Taylor, S.R., 1989. Intraplate Volcanism in Eastern Australia and New Zealand. Cambridge University Press, Melbourne. Karato, S., Jung, H., 1998. Water, partial melting and the origin of the seismic low velocity and high attenuation zone in the upper mantle. Earth Planet. Sci. Lett., 157, 193–207. Kawakatsu, H., Kumar, P., Takei, Y., Shinohara, M., Kanazawa, T., Araki, E., Suyehiro, K., 2009. Seismic evidence for sharp lithosphere-asthenosphere boundaries of oceanic plates, Science, 324, 499–502. Keith, C. M., Crampin, S., 1977. Seismic body waves in anisotropic media: Synthetic seismograms, Geophys. J. R. Astron. Soc., 49, 225–243. Kennett, B.L.N., 1991. The removal of free surface interactions from three-component seismograms. Geophys. J. Int., 104, 153-163. Kennett, B.L.N., Engdahl, E.R., Buland, R., 1995. Constraints on seismic velocities in the Earth from travel-times. Geophys. J. Int., 122, 108-124. Kikuchi, M., Kanamori, H., 1982. Inversion of complex body waves. Bull. Seismol. Soc. Am., 72, 491-506. King, S.D., Ritsema, J., 2000. African hot spot volcanism: Small-scale convection in the upper mantle beneath cratons. Science. 290, 1137-40. Korenaga, J., Jordan, T.H., 2002. On the state of sublithospheric upper mantle beneath a supercontinent. Geophys. J. Int., 149, 179–89. Kumar, P., Kind, R., Hanka, W., Wylegalla, K., Reigber, C., Yuan, X., Woelbern, I., Schwintzer, P., Fleming, K., Dahl-Jensen, T., Larsen, T.B., Schweitzer, J., Priestley, K., Gudmundsson, O., Wolf, D., 2005b. The lithosphere-asthenosphere boundary in the North-West Atlantic region. Earth Planet. Sci. Lett., 236, 249-257. Kumar, P., Yuan, X., Kind, R., Kosarev, G., 2005a. The lithosphere-asthenosphere boundary in the Tien Shan-Karakoram region from S receiver functions: Evidence for continental subduction. Geophys. Res. Lett., 32, L07305. doi:10.1029/2004GL022291. Kumar, P., Yuan, X.H., Kumar, M.R., Kind, R., Li, X.Q., Chadha, R.K., 2007. The rapid drift of the Indian tectonic plate. Nature, 449, 894-897. Lee, C.T.A., 2003. Compositional variation of density and seismic velocities in natural peridotites at STP conditions: Implications for seismic imaging of compositional heterogeneities in the upper mantle. J. Geophys. Res., 108. doi:10.1029/2003JB002413. Lekic and Romanowicz, submitted to Nature Geoscience. Tectonic regionalization without a priori information: a cluster analysis of tomography. Li, X., Kind, R., Priestley, K., Sobolev, S.V., Tilmann, F., Yuan, X., Weber, M. 2000. ! 28! Mapping the Hawaiian plume conduit with converted seismic waves. Nature, 405, 938-941. Li, X., Kind, R., Yuan, X.H., Wolbern, I., Hanka, W., 2004. Rejuvenation of the lithosphere by the Hawaiian plume. Nature, 427, 827-829. Li, X., Yuan, X., Kind, R., 2007. The lithosphere-asthenosphere boundary beneath the western United States. Geophys. J. Int., 170, 700-710. Ligorria, J.P. Ammon, C.J., 1999. Iterative deconvolution and receiver-function estimation. Bull. Seismol. Soc. Am., 89, 1395-1400. Levin V., Park J., 2000. Shear zones in the Proterozoic lithosphere of the Arabian Shield and the nature of the Hales discontinuity. Tectonophysics 323, 131-48. Mercier, J.-P., Bostock, M.G., Audet, P., Gaherty, J.B., Garnero, E.J., Revenaugh, J., 2008. The teleseismic signature of fossil subduction: Northwestern Canada. J. Geophys. Res., 113, B04308. doi:10.1029/2007JB005127. Mohsen, A., Kind R., Sobolev, S.V., Weber, M., DESERT Group, 2006. Thickness of the lithosphere east of the Dead Sea Transform. Geophys. J. Int., 167, 845-852. Nettles, M. Dziewonski, A.M., 2008. Radially anisotropic shear velocity structure of the upper mantle globally and beneath North America. J. Geophys. Res., 113, 61-67. O’Reilly, S.Y., Griffin, W.L., Gaul, O., 1997. Paleogeothermal gradient in Australia: key to 4-D lithosphere mapping. J. Aust. Geol. and Geophys., 17, 63-72. Oreshin, S., Vinnik, L., Peregoudov, D., Roecker, S., 2002. Lithosphere and asthenosphere of the Tien Shan imaged by S receiver functions. Geophys. Res. Lett., 29, 1191. doi:10.1029/2001GL014441. Ozacar, A.A., Gilbert, H., Zandt, G., 2008. Upper mantle discontinuity structure beneath East Anatolian Plateau (Turkey) from receiver functions. Earth Planet. Sci. Lett., 269, 426-434. Pedersen, H.A., Fishwick, S., Snyder, D.B., 2009. A comparison of cratonic roots through consistent analysis of seismic surface waves. Lithos, 109, 81-95. Reading, A.M., Kennett, B.L.N., Goleby, B., 2007. New constraints on the seismic structure of West Australia: Evidence for terrane stabilization prior to the assembly of an ancient continent?,.Geology, 35, 379-382. Reading, A.M, Kennett, B.L.N., 2003. Lithospheric structure of the Pilbara Craton, Capricorn Orogen and northern Yilgarn Craton, Western Australia, from teleseismic receiver functions. Aust. J. Earth Sci., 50, 439-445. Reading, A.M., Kennett, B.L.N., Dentith, M.C., 2003. Seismic structure of the Yilgarn Craton, Western Australia, Aust. J. Earth Sci., 50, 427-438. Revenaugh, J.S., Jordan, T.H., 1991. Mantle layering from ScS reverberations: 3. The upper mantle. J. Geophys. Res., 96, 19781–19810. Romanowicz, B., 2009. The thickness of tectonic plates. Science, 324, 474–476. Rychert, C.A., Fischer, K.M., Rondenay, S., 2005. A sharp lithosphere-asthenosphere boundary imaged beneath eastern North America. Nature, 436, 542-545. Rychert, C.A., Rondenay, S., Fischer, K.M., 2007. P-to-S and S-to-P imaging of a sharp lithosphere-asthenosphere boundary beneath eastern North America. J. Geophys. Res., 112, B08314. doi:10.1029/2006JB004619. Rychert, C.A., Shearer P.M., 2009. A Global View of the Lithosphere-Asthenosphere Boundary, Science, 324, 495-498. Rychert, C.A., Shearer, P.M., Fischer, K.M., 2010. Scattered Wave Imaging of the ! 29! Lithosphere-Asthenosphere Boundary. Lithos, in press. Sacks, I.S., Snoke, J.A., Husebye, E.S., 1979. Lithosphere thickness beneath the Baltic shield. Tectonophysics, 56, 101-110. Simons, F.J., van der Hilst, R.D., 2002. Age-dependent seismic thickness and mechanical strength of the Australian lithosphere. Geophys. Res. Lett., 29, 1529–1533. Simons, F.J., van der Hilst, R.D., 2003. Seismic and mechanical anisotropy and the past and present deformation of the Australian lithosphere. Earth Planet. Sci. Lett., 211, 271-286. Simons, F.J., Zielhuis, A., van der Hilst, R.D., 1999. The deep structure of the Australian continent from surface-wave tomography. Lithos, 48, 17-43. Snyder, D.B., 2008. Stacked uppermost mantle layers within the Slave craton of NW Canada as defined by anisotropic seismic discontinuities. Tectonics, 27, TC4006. doi:10.1029/2007TC002132. Sodoudi, F., Yuan, X., Liu, Q., Kind, R., Chen J., 2006. Lithospheric thickness beneath the Dabie Shan, central eastern China from S receiver functions. Geophys. J. Int., 166, 1363-1367. Takei, Y., 2002. Effect of pore geometry of Vp/Vs: From equilibrium geometry to crack. J. Geophys. Res., 107, B22043 doi:10.1029/2001JB005850. Takei, Y., Holtzman, B.K., 2009. Viscous constitutive relations of solid-liquid composites in terms of grain boundary contiguity: 1. Grain boundary diffusion control model. J. Geophys. Res., 114, B06205. doi:10.1029/2004JB002965. Thybo, H., 2006. The heterogeneous upper mantle low velocity zone. Tectonophysics, 416, 53- 79. Tyler, I.M., 2001. Collisional orogeny during the Paleoproterozoic in Western Australia. Geol. Soc. Aust. Abstr., 64, 187-188. Vinnik, L., Kurnik, E., Farra, V., 2005. Lehmann discontinuity beneath North America: No role for seismic anisotropy. Geophys. Res. Lett., 32, L09306. doi:10.1029/2004GL022333. Wellman, P., 1983. Hotspot volcanism in Australia and New Zealand: Cainozoic and mid- Mesozoic. Tectonophysics, 96, 225-243. Wellman, P., 1998. Mapping of geophysical domains in the Australian continental crust using gravity and magnetic anomalies, in: Braun, J., Dooley, J., Goleby, B., van der Hilst, R., Klootwijk, C. (Eds.), Structure and Evolution of the Australian Continent. American Geophysical Union Geodynamics Series, 26, 59-71. Wellman, P., McDougall, I., 1974. Cainozoic igneous activity in eastern Australia, Tectonophysics, 23, 49-65. Wilson, D., Grand, S., Angus, D., Ni, J., 2006. Constraints on the interpretation of S-to-P receiver functions. Geophys. J. Int., 165, 969-980. Wittlinger, G., Farra, V., 2007. Converted waves reveal a thick and layered tectosphere beneath the Kalahari super craton. Earth Planet. Sci. Lett., 254, 404-415. Wolbern, I., Jacob, A.W.B., Blake, T.A., Kind, R., Li, X., Yuan, S., Duennebier, F., Weber, M., 2006. Deep origin of the Hawaiian tilted plume conduit derived from receiver functions. Geophys. J. Int., 166, 767-781. Yuan H., Romanowicz, B., 2010. Lithospheric layering in the North American Craton. Nature, submitted. ! 30! Yuan, X., Kind, R., Li, X., Wang, R., 2006. The S receiver functions: Synthetics and data example. Geophys. J. Int., 165, 555-564. Zhu, L., Kanamori, H., 2000. Moho depth variation in southern California from teleseismic receiver functions. J. Geophys. Res., 105, 2969-2980. ! 31! 120° E 150° E Panel A Legend Archean COEN COEN Proterozoic 15° S Phanerozoic FITZ WRAB CTAO CTAO Panel B Legend MBWA Volume (cubic km) EIDS < 100 EIDS 100-500 500-1000 30° S BLDU ARMA > 1000 FORT ARMA MUN KMBL STKA STKA NA NWAO BBOO YNG YNG Average Age (Ma) CAN < 10 TOO TOO 10-20 20-30 30-40 MOO 40-50 TAU >50 A B 45° S Figure 1. (a) Overview of major tectonic provinces within Australia; solid black lines outline the seven major geophysical domains [Wellman, 1998]. Region with inclined lines is the West Australia Craton (WC); horizontal lines indicate the North (NC) and South Australia (SC) Cratons; dots indicate the presence of Phanerozoic basement inferred from the Tasman Line (Gunn et al., 1997). The nineteen stations for which Ps and Sp receiver functions are calculated are marked with an inverted triangle. Black inverted triangles indicate the receiver functions are included in later interpretation. (b) Volcanism related to central and lava field volcanism (from Johnson et al. [1989]). Volcanic age is indicated by color and volume is shown by size of the symbol. Star and oval show location of most recent volcanism. ! 32! WRAB Figure 2. Distribution of events used to calculate Ps (blue circles) and Sp (red circles) receiver functions for station WRAB. A total of 523 events from epicentral distances of 35°-80° were used for Ps at WRAB and 121 events from epicentral distances of 55°-75° were used for Sp at WRAB. ! 33! ARMA (AU) Single Sp Frequency Domain RF 112 Events 0.2 Fraction of SV Component Amplitude 0.15 0.1 0.05 0 −0.05 −0.1 −0.15 −0.2 0 100 200 300 Depth (km) ARMA (AU) Single Sp Time Domain RF 112 Events 0.1 Fraction of SV Component Amplitude 0.05 0 −0.05 −0.1 0 100 200 300 Depth (km) Figure 3. Comparison of single binned Sp receiver functions for station ARMA calculated using frequency-domain (a) and time-domain (b) deconvolution methods. The thick, solid black line is the mean calculated receiver function from the bootstrap test. The mean is virtually identical to the single-stacked receiver function. The solid grey lines on either side of the mean are the two standard deviations. The receiver functions are similar in form, with a strong positive phase at a depth of 34±5 km in the frequency domain and 33±5 km in the time domain, which is interpreted as the Moho. At 93±16 km (a) and 81±18 km (b) a robust negative phase corresponds to the estimated LAB depth from the shear wave velocity model of Fishwick et al. [2008a, 2008b]. A phase is only considered robust if there is good agreement between the time and frequency domains. ! 34! TOO (AU) Single Ps FORT (AU) Single Ps Frequency Domain RF Frequency Domain RF a) 216 Events c) 169 Events 0.08 0.5 Fraction of P Component Amplitude Fraction of P Component Amplitude 0.06 0.04 0.02 0 0 −0.02 −0.04 −0.06 −0.08 −0.5 0 50 100 150 200 250 0 100 200 300 Depth (km) Depth (km) TOO (AU) Single Sp FORT (AU) Single Sp Frequency Domain RF Frequency Domain RF b) 91 Events d) 44 Events 0.2 0.1 Fraction of SV Component Amplitude Fraction of SV Component Amplitude 0.15 0.1 0.05 0.05 0 0 −0.05 −0.1 −0.05 −0.15 −0.1 −0.2 0 50 100 150 200 250 0 100 200 300 Depth (km) Depth (km) Figure 4. Example of single bin Ps and Sp receiver functions calculated for stations TOO (a and b) and FORT (c and d). The Ps receiver function for TOO (a) exhibits a strong, well-defined positive phase at 34±3 km indicating a velocity increase with depth, interpreted as the Moho. A negative phase is observed at 56±2 km and is consistent with the negative phase in Sp (b) at a depth of 61±11 km. The Ps receiver function at FORT (c) has complicated crustal structure with numerous reverberations, in contrast to the Sp receiver function at FORT (d), which exhibits a clear signal, with a strong negative phase indicating a velocity decrease with depth at 79±6 km. This negative phase is not apparent on the Ps receiver function (c), which may be due to interference from reverberations. ! 35! WRAB (II) Single Sp STKA (AU) Single Sp COEN (AU) Single Sp Frequency Domain RF Frequency Domain RF Frequency Domain RF 121 Events 97 Events 70 Events 0.15 0.15 0.2 Fraction of SV Component Amplitude Fraction of SV Component Amplitude Fraction of SV Component Amplitude 0.1 0.15 0.1 0.1 0.05 0.05 0.05 0 0 0 −0.05 −0.05 −0.05 −0.1 −0.1 −0.1 −0.15 −0.15 −0.2 −0.15 0 100 200 300 0 100 200 300 0 100 200 300 Depth (km) Depth (km) Depth (km) WRAB (II) Binned Sp STKA (AU) Binned Sp COEN (AU) Binned Sp Frequency Domain RF Frequency Domain RF Frequency Domain RF 121 Events 97 Events 70 Events Epicentral Distance (degrees) 80 Epicentral Distance (degrees) 80 80 Epicentral Distance (degrees) 75 75 75 70 70 70 65 65 65 60 60 60 55 55 55 50 50 50 0 100 200 300 0 100 200 300 0 100 200 300 Depth (km) Depth (km) Depth (km) Figure 5. Example of Sp receiver functions for three tectonically distinct regions. The thick, solid black line is the mean calculated receiver function from bootstrapping, where the solid grey lines on either side of the mean are the two standard deviations. (a-b) WRAB is located well within the interior of the Proterozoic North Australia Craton and is characterized by a negative phase at 81±14 km that is interpreted to be mid-lithospheric discontinuity. (c-d) STKA is found along the ambiguously defined Proterozoic margin. The Sp receiver function is characterized by a well-constrained negative phase at 104±9 km inferred to be the LAB. (e-f) The results at station COEN are typical of many receiver functions for stations along the eastern margin of the continent. A negative at 67±8 km is interpreted to be the LAB. ! 36! Distance (km) 120° E 150° E A 0 500 1000 1500 2000 2500 3000 A’ a) b) TOO YNG ARMA EIDS CTAO COEN A’ 0 15° S COEN C FITZ C’ 50 WRAB CTAO MBWA Depth (km) EIDS 100 30° S FORT BNWAOKMBL STKA ARMA BBOO B’ 150 YNG TOO A 200 250 Distance (km) Distance (km) B 0 500 1000 1500 2000 2500 3000 B’ C 0 500 1000 1500 2000 2500 3000 C’ c) NWAO KMBL FORT BBOO STKA YNG d) MBWA FITZ WRAB CTAO 0 0 50 50 Depth (km) Depth (km) 100 100 150 150 200 200 250 250 Figure 6. (a) Lines demonstrating location of cross-sections of Sp receiver functions A- A’, B-B’ and C-C’. (b thru d) The mean of the bootstrapped receiver functions are shown in solid black in cross-sections A-A’, B-B’ and C-C’. The receiver functions are plotted to the same scale in each cross section and the statistically significant portions are represented in either red (positive) or blue (negative). Depths are shown to 250 km, but individual receiver functions were examined to 400 km to ensure that no identifiable phases existed at greater depths. Black horizontal lines mark the location of the largest significant negative phase, and the surrounding grey box indicates the two standard deviations. Solid grey and black circles indicate the range of potential LAB depths from surface wave tomography [Fishwick et al., 2008a, 2008b]. If the negative phase from the receiver function falls within the potential LAB depth range, the phase is interpreted to be the LAB; otherwise it is characterized as a MLD. ! 37! 120° E 150° E a) 15° S COEN FITZ MBWA WRAB CTAO MLD LAB EIDS 30° S ARMA KMBL FORT STKA NWAO BBOO YNG TOO 45° S 50 60 70 80 90 100 110 120 130 Depth (km) 120° E 150° E b) 15° S COEN FITZ MBWA WRAB CTAO MLD LAB EIDS 30° S ARMA KMBL FORT STKA NWAO BBOO YNG TOO 45° S −0.12 −0.11 −0.1 −0.09 −0.08 −0.07 −0.06 −0.05 −0.04 −0.03 Amplitude (Relative to Parent Phase) Figure 7. Smoothed conversion depths (a) and amplitudes (b) of selected negative Sp phases. Negative phases are selected based on the largest, statistically significant negative peak within a given receiver function. The smoothed conversion point locations are created by first calculating the theoretical piercing points of the selected negative phase for each station using the respective H-k stacking crustal model combined with AK135 [Kennett et al., 1995] for the mantle. The piercing points were then placed onto a 0.05° x 0.05° grid and were averaged if more than one piercing point fell onto a single element. The piercing points were then smoothed with a circular filter, which produced a new spacing of 0.3°. Upside-down black triangles indicate seismic station location. The thick dashed black line is used to graphically illustrate the separation between stations imaging the LAB and stations imaging a (MLD). ! 38! 0.16 0.14 |Amplitude| Relative to the Parent Phase 0.12 0.1 0.08 0.06 0.04 LAB 0.02 MLD 0 0 50 100 150 Depth (km) Figure 8. Plot of amplitude versus depth of the negative phase interpreted to be the LAB (solid black squares) and the MLD (black and white squares) shown with two standard deviation error bars. Black line shows the best fit through the LAB phases. ! 39! a) 0.25 b) 0.25 10% Velocity Drop 0 km 10% Velocity Drop 0 km 0.2 COEN 10 km ARMA 10 km 20 km 0.2 20 km 30 km 30 km 0.15 40 km 40 km 50 km 0.15 50 km 0.1 Amplitude Amplitude 0.1 0.05 0.05 0 0 −0.05 −0.1 −0.05 −0.15 −0.1 0 50 100 150 0 50 100 150 Depth (km) Depth (km) Figure 9. (a) Solid black line is the mean of the single-binned receiver functions for station COEN and the dashed lines are two standard deviations. The colored lines are the synthetic receiver functions calculated for a 10% velocity dropped over a gradient thickness of 0 to 50 km. (b) Same as part (a) but for station ARMA. ! 40! W E 100 1 3 Mid-lithospheric discontinuity Sharp LAB Depth (km) 200 300 2 Gradual LAB 400 Figure 10. Schematic cross-section through Australia intersecting Archean, Proterozoic and Phanerozoic terranes. LAB depth estimated from surface wave tomography [Fishwick et al., 2005] and Sp receiver functions (this study). ! 41! Events Crustal Velocity Model Phases and Interpretation Negative LAB depth Net STA #Ps #Sp Model Moho Vp Vp/Vs Moho Phase range Boundary AU ARMA 199 112 H-k Stack 34 6.42 1.72 36±3 93±16 <150 km LAB AU BBOO 192 58 H-k Stack 41 6.55 1.77 42±3 131±9 150-175 km LAB AU BLDU 13 2 H-k Stack 35 6.42 1.73 36±2 - - - AU COEN 266 70 H-k Stack 37 6.47 1.7 38±2 67±8 Absent Lid LAB AU EIDS 209 112 H-k Stack 35 6.51 1.7 37±3 76±12 <150 km LAB AU FITZ 279 59 H-k Stack 41 6.45 1.71 44±5 81±8 125-225 km MLD AU FORT 169 44 Fixed 40 6.45 1.85 40±9* 79±6 125-200 km MLD AU KMBL 294 64 H-k Stack 37 6.46 1.73 38±1 85±14 125-225 km MLD AU MOO 84 39 H-k Stack 33 6.27 1.68 33±3 - - - AU MUN 151 34 H-k Stack 51 6.42 1.87 53±2 - - - AU STKA 333 97 H-k Stack 43 6.47 1.68 46±3 104±9 100-175 km LAB AU TOO 216 91 H-k Stack 33 6.43 1.77 34±3 61±11 Absent Lid LAB AU YNG 184 79 H-k Stack 34 6.57 1.83 33±2 70±8 <150 km LAB II TAU 313 120 H-k Stack 31 6.18 1.74 32±3 - - - II WRAB 523 121 H-k Stack 48 6.58 1.74 49±2 81±14 175-200 km MLD IU CTAO 619 251 H-k Stack 39 6.51 1.7 40±2 73±6 Absent Lid LAB IU MBWA 349 129 H-k Stack 31 6.51 1.69 32±2 69±8 100-200 km MLD IU NWAO 372 176 H-k Stack 38 6.46 1.78 41±3 81±8 100-175 km MLD Table 1. Summary of receiver function results at the nineteen stations used in the study. The velocity model used for the migration for all but FORT was created using a crustal model from Ps H-k stacking for the respective station and a mantle model from AK135. The Moho phase depths listed in the table were obtained from Ps receiver functions, except for FORT* which was obtained from the Sp receiver function. The largest statistically significant negative phase for each station was selected from Sp receiver functions. Phase depth uncertainties correspond to the maximum range of depth in between where the bootstrap two standard deviations have the same amplitude as the phase peak of the bootstrap mean. A negative phase is not listed for stations where the time domain and frequency domain deconvolution methods yielded significantly different receiver function results. The potential “LAB depth range” was determined from the shear wave velocity model [Fishwick et al., 2008a, b] and corresponds to the depth range between the minimum Vs beneath the Moho up to the next local maximum in velocity. Stations where there is no high velocity lid present beneath 75 km are designated with the label “Absent Lid” and the LAB depth range is assumed to be between the Moho and 75 km. If the negative phase pick from the Sp receiver functions falls within the LAB depth range, the phase is interpreted as the LAB otherwise it is interpreted to be a mid-lithospheric discontinuity (MLD). BBOO is the one exception. See section 3 for a more detailed discussion. ! 42! Chapter 2 The effects of anisotropy on Sp receiver functions By: Heather A. Ford1 and Karen M. Fischer1 1 Department of Geological Sciences, Brown University ! 43! Abstract Sp receiver functions have become a commonly used method for imaging upper mantle structure. However, there has been little work done to understand the effects of anisotropy on Sp receiver functions, despite the evidence for anisotropic discontinuities in the upper mantle. Using synthetic seismograms calculated for a number of simple anisotropic upper mantle models, we illustrate the effects of layered azimuthal and radial anisotropy on Sp receiver functions and compare the Sp results to Ps receiver functions. We find that vertical variations in azimuthal anisotropy produce distinct back-azimuthal patterns in Sp receiver functions, which could be useful in identifying anisotropy in observed Sp phases. However, we also document a strong dependence of Sp receiver function on the initial polarization of the direct-S phase, a complexity that does not exist with Ps receiver functions. Resulting Sp receiver functions that include radial anisotropy with a vertical symmetry axis do not have the same dependence on back azimuth and incident polarizaition; instead a coherent positive or negative phase is produced, depending on whether the vertical axis is fast or slow. Sp receiver functions observed within the Canadian craton, including stations newly analyzed in this study, contain evidence for a mid-lithospheric discontinuity, which appears to have a negative amplitude (using the Ps polarity convention) indicative of a velocity decrease with depth if only isotropic structure is considered. At station FCC, located within the craton, observed SKSp receiver functions show a negative phase at depths of ~125 km. The back-azimuthal signature of this phase is matched by synthetic SKSp receiver functions for an anisotropic model based on the SAWum_NA2 long- ! 44! period waveform tomography of Yuan et al. (2011). In this model, the negative SKSp phase is produced by a vertical change in the orientation and strength of azimuthal anisotropy, rather than a decrease in isotropic velocity. This result opens the door to an interpretation in which the mid-lithospheric discontinuity observed widely in the craton represents layering in azimuthal anisotropy, as opposed to the top of mid-lithospheric low velocity layer. ! 45! 1. Introduction Receiver function analysis is a seismic imaging technique that is used to identify discontinuity structure in the crust, upper mantle and mantle transition zones (e.g., Vinnik, 1977). The technique uses the partial scattering of a primary P- or S- seismic phase at a boundary in seismic properties, into an S- or P- wave, respectively, to infer information about the scatterer. The term “receiver function” arises from the process of removing source effects, through deconvolution, to isolate information about seismic structure beneath the receiver. Assuming a velocity structure, receiver functions can be migrated and the depth of the scatterer can be determined. Additional analysis, such as modeling of the receiver functions, can constrain information regarding the seismic velocity gradient associated with the boundary (e.g., Rychert et al., 2007; Ford et al., 2010) and binning receiver functions as a function of back-azimuth can aid in the identification of azimuthal anisotropy. The behavior of Ps receiver functions in the presence of anisotropy has been well documented (e.g., Bostock, 1998; Fredericksen and Bostock, 2000; Levin and Park, 1997; Savage, 1998) and used to confirm the existence of anisotropic structure in a number of tectonic settings (e.g., Bostock, 1998; Mercier et al., 2008; Ozacar and Zandt, 2009; Savage, 1998; Schulte-Pelkum et al., 2005; Snyder, 2008; Wirth and Long, 2012). Receiver functions have a distinct advantage over shear-wave splitting methods, which rely on the path-integrated observations of anisotropy and therefore cannot directly image boundaries in anisotropy (e.g., Silver and Chan, 1991; Silver, 1996; Savage, 1999). ! 46! Similarly, while surface wave tomography has the ability to delineate variations in anisotropy with depth (e.g., Yang and Forsyth, 2006), its vertical resolution of sharp velocity gradients (velocity interfaces) is not as precise as receiver function analysis, which typically has an uncertainty of ±10 km in depth (Lekic et al., 2011). One large limitation of Ps receiver function imaging is the appearance of crustal reverberations, which often arrive in the same time window as upper mantle phases. In contrast, the converted P-wave in Sp receiver functions is separated from the reverberations by the primary S-phase. As a result, Sp receiver function analysis does not suffer from issues of reverberation contamination, making the method ideal for imaging upper mantle structure, such as the lithosphere-asthenosphere boundary (LAB) (e.g., Li et al., 2007; Rychert et al., 2007; Abt et al., 2010; Miller and Eaton, 2011; Lekic et al., 2011; Kind et al., 2012; Kumar et al., 2012; Levander and Miller, 2012). Despite the rising popularity of Sp receiver function analysis, and the observation of boundaries in anisotropy in the upper mantle (e.g. Bostock, 1998), the interaction of Sp phases and receiver functions with anisotropic structure has only begun to be explored (Bourguignon, 2009). To improve our understanding of these effects, we calculated synthetic Sp receiver functions for a range of simple anisotropic models, testing the impact of anisotropy strength and geometry, and comparing these results to Ps receiver functions. Using the information gained from our simple modeling, we then evaluated whether the negative mid-lithospheric phase commonly observed in Sp receiver functions within the Canadian craton can be explained by invoking vertical variations in ! 47! anisotropy. To test this hypothesis, we calculated synthetic SKSp receiver functions for a velocity model that represents a simplified version of the azimuthally and radially anisotropic SAWum_NA2 model (Yuan et al., 2011) for station FCC, located within the Canadian craton, and we compared the synthetics to the observed SKSp receiver function for station FCC. 2. The interaction of Sp receiver functions with anisotropic structure 2.1 Generating synthetic seismograms In order to characterize the behavior of Sp receiver functions in the presence of anisotropy, we first calculated synthetic seismograms for a range of simple models using a propagator matrix method (Keith and Crampin, 1977), which allows for wave propagation through an n-layered, horizontally stratified medium, with anisotropic and/or isotropic elastic coefficients. The synthetic seismograms were calculated for cases at distances of 75° and 110° (corresponding to epicentral distances of S- and SKS- phase arrivals respectively), depths of 0-300 km, and back-azimuths of 0-360°. Since the initial polarization direction for the S-phase is dependant on focal mechanism, which varies between earthquakes, we also vary the amount of incident SV-SH energy present. 2.2 Elastic coefficients Although the upper mantle is primarily composed of olivine and orthopyroxene, both of which are orthorhombic, studies of naturally- and experimentally-deformed olivine ! 48! aggregates indicate that assuming a hexagonal symmetry is reasonable (Mainprice and Silver, 1993; Bystricky et al., 2000; Mehl et al., 2003; Michibayashi et al., 2006). Therefore, we employ hexagonal anisotropy, which requires only five independent elastic coefficients, in all of the simple models examined in this section. (A more general anisotropic symmetry is employed when modeling structure beneath station FCC in Canada in the next section.) The elastic coefficients employed in this section (listed in Table 1 as “SPB – full” ) were taken from Schulte-Pelkum and Blackman (2003) and are thought to be representative of upper-most mantle velocities. The amount of P anisotropy in SPB is set at 10% and the S anisotropy is set at 9.5%. This value is considered to be at the high-end of the possible range in upper mantle anisotropy inferred from observations (e.g., Smith and Ekström, 1999), and is used in this paper for illustrative purposes only. In the elastic coefficient case labeled “SPB – half”, the amount of anisotropy is reduced by half by diluting the original elastic coefficients with the equivalent isotropic coefficients that have the same Voigt-averaged velocity (Babuska and Cara, 1991). The strength of anisotropy in the second case is closer to values that have been measured in the upper mantle (e.g. Forsyth et al. 1998; Lin et al., 2010; Yuan and Romanowicz, 2010; Yuan et al., 2011). Table 1 lists the elastic coefficients, thickness of layers, and the azimuth and plunge of the fast-axis in layers with anisotropy. For both SPB-full and SPB-half, the axis of symmetry is the fast axis. In section 3.3, we explore the effects of a slow symmetry axis on radial anisotropy. 2.3 Calculating synthetic receiver functions ! 49! For this study, we calculated Ps, Sp and SKSp receiver functions. Waveforms were filtered with a band-pass of 0.03 to 0.5 Hz, windowed around the parent phase (S for Sp; SKS for SKSp; P for Ps), and rotated into P, SV and SH components using a free surface transfer matrix (Kennett, 1991) in which the values of Vp and Vs are set to match the velocities of the crustal layer used to construct the synthetic seismograms (Table 1). The component containing the scattered phase was then deconvolved by the corresponding parent phase using a simultaneous frequency domain technique (Bostock, 1998) and migrated to depth. For Ps receiver functions, the P component is deconvolved from the SV component. For the SKSp cases and most of the Sp examples, the SV component was deconvolved from the P component (“radial component” receiver functions). However, in a few cases “transverse component” Sp receiver functions were also generated (SH deconvolved from P). The polarity of the Sp and SKSp receiver functions was flipped in order to match the sign of the Ps receiver functions. 3. Synthetic receiver function results 3.1 Comparison to Ps receiver functions Our first modeling case, Model 1 (Table 1b), uses the full strength hexagonal anisotropy coefficients from Schulte-Pelkum and Blackman (2003) for the mantle layer and the mantle half-space. However, the fast-axis direction between the two layers is offset by 90°, so that the fast-axis in the top mantle layer is oriented N-S, and the fast-axis in the mantle half-space is oriented E-W (Fig. 1). For the Sp case here (Model 1) and in ! 50! Models 2-4, incident S motion is SV and the Sp receiver functions are shown for SV deconvolved from P. Both the Ps (Fig. 2a) and Sp (Fig. 2b) receiver functions have a positive Moho phase at 36 km and a phase located at 100 km due to the interface in anisotropy in Model 1 (Fig. 1). The Ps receiver function also exhibits a clear set of crustal reverberations at depths of ~150 km and ~190 km, which demonstrates the potential difficulty with using Ps receiver functions to image upper mantle structure. For both Ps and Sp, there is a periodicity in the mantle phase at 100 km every 180° in back-azimuth. The periodicity in Ps receiver function amplitudes has been well documented by others (e.g., Levin and Park, 1998; Savage, 1998), and is often referred to as a sin(2θ) pattern. An isotropic dipping interface can also produce variations in back- azimuth. However, careful comparison, which includes the use of the P-SH receiver functions, can help to distinguish between boundaries in anisotropy and dipping interfaces (e.g., Wirth and Long, 2012). The second modeling case, Model 2 (Table 1b) uses the same geometry as Model 1, but with elastic coefficients that have been diluted (SPB – half); incident S particle motion remains SV. The overall behavior in Ps and Sp is similar to that with Model 1, however the change in amplitude and depth is more subtle in both Ps and Sp receiver functions (Fig. 2c and 2d). In particular, the back-azimuthal variation in depth is significantly ! 51! smaller in Sp, and may be difficult to observe in a real world data example. The amplitude variation in Ps although smaller, is still pronounced. 3.2 Variations in anisotropy geometry In the next two modeling cases, Models 3 and 4, we used SPB - half and SV incident S particle motion, while varying the azimuth and plunge of the fast-axis (Fig. 4). In Model 3, the top mantle layer retains the same fast-axis direction of N-S, while the mantle half- space now has a fast-axis direction of 45° (NE-SW). Both layers maintain a fast-axis plunge of 0°. The resulting Ps receiver function (Fig. 4a) has a pattern similar to what is observed in Models 1 and 2, with the difference being that there is a back-azimuthal shift in the zero crossing of the sin(2θ) amplitude pattern. In contrast, there is a striking change in the Sp receiver functions (Fig. 4b) when the fast- axis azimuth offset of between the two layers is changed from 90˚ to 45˚. The new Sp receiver function pattern in back-azimuth is characterized by a rapid change in the sign of the phase at back-azimuths of 45° and 225°, and more gradual changes between 90- 130° and 270-310°. This variation in amplitude can be thought of as a cot(2θ) pattern. The difference in behavior between Models 1-2 and Model 3, demonstrates the potential utility in using Sp receiver functions to differentiate between models with varying fast- axis azimuth. The assumption of anisotropy with a plunging symmetry axis in the upper mantle is a reasonable one, since plunging anisotropy has been invoked to explain observations in a ! 52! number of tectonic settings (e.g. Bostock, 1998; Mercier et al., 2008; Snyder, 2008; MacDougall et al., 2012). In Model 4, the azimuths of the fast-axes are the same as in Model 3, with the only difference being that the plunge of the fast-axis is changed from 0° to 20° in the top mantle layer, with the plunge downward to the north. The resulting difference in both Ps and Sp receiver functions (Fig. 4c and 4d) is the loss of the 180° periodicity, which was observed in Models 1 through Model 3. In both Ps and Sp cases, the amplitude of the phase from the interface at 100 km is reduced in the direction of a- axis plunge. We do not consider the possibility of a dipping interface, either isotropic or anisotropic, but synthetic examples can be seen in numerous studies for Ps receiver functions (e.g. Fredericksen and Bostock, 2000; Mercier et al., 2008; Wirth and Long, 2012) and in Bourguignon (2009) for Sp receiver functions. 3.3 Sp receiver functions in the presence of radial anisotropy In Models 5 and 6 we investigate the effects of radial anisotropy with the axis of symmetry oriented in the vertical direction. Radial anisotropy exists in the upper layer, and In contrast to Models 1-4, the half space of Models 5 and 6 correspond to the isotropic Voigt-averaged velocities of the anisotropic coefficients from the same model. Additional epicentral distances (57.5-77.5°) are also included. In Model 5, the elastic coefficients of SPB – half are used for the anisotropic layer, and the Voigt-averaged coefficients (SP1 – iso) are used in the half space. In Model 6, the anisotropic coefficients correspond to the “hexagonal slow” coefficients given in Schulte-Pelkum and ! 53! Blackman (2003), and diluted by 50% (Fig. 1). We label these coefficients as SPB – slow, since the axis of symmetry is the slow direction, and the plane perpendicular to the slow direction is fast. The isotropic half space in Model 6 is the equivalent Voigt- averaged velocities (SP2 – iso). The coefficients for SP1 – iso, SP2 – iso, and SPB – slow are listed in Table 1. In Models 5 and 6, and in contrast to Models 1-4, there is no change in amplitude with back azimuth (Fig. 3) as is expected for radial anisotropy with a vertical symmetry axis. For the case of an all-SV incident phase in Model 5, where the elastic coefficients in the anisotropic layer are oriented so that the fast axis is vertical, the observed phase at 100 km depth is positive. The corresponding phase observed in Model 6 is negative (Fig. 3), due to the differing form of anisotropy (slow axis vertical). 3.4 Coupling between P, SV and SH In a purely isotropic system with non-dipping interfaces, energy is coupled between P and SV, but not with SH. However, when anisotropy (or a dipping interface) is introduced, P-SV-SH coupling can occur. Ps receiver functions have utilized this interaction by examining the SH/Transverse component of Ps receiver functions (e.g., Levin and Park, 1997; Bostock, 1998; Fredericksen and Bostock, 2000). In contrast, while the interaction between incident SV and SH particle motion in Sp receiver functions was noted by Bourguignon (2009), it has not been as extensively explored. In this study we consider the anisotropic interaction between P, SV and SH in Sp synthetic receiver functions by determining the impact that differing amount of incident SV and ! 54! SH energy can have on Sp receiver functions, and by calculating both radial and transverse component Sp receiver functions. The amount of SV and SH particle motion in a given Sp phase varies according to the source mechanism radiation pattern, although anisotropy along the path, on both the receiver- and the source- side, can also introduce variations in the ratio of SV to SH energy. To understand the effect that incident SH versus SV has on Sp receiver functions we tested two end-member cases, where the amount of incident SV to SH particle motion varies from all SV, to all SH, as well as an intermediate case where the incident SV and SH amplitudes are equal. These cases are all calculated using Model 3, where the offset in fast-axis direction between the two mantle layers is 45°, but where the fast-axes are aligned horizontally. Striking differences are observed as a function of the ratio of incident SV and SH particle motion. In the all-SV incidence case for the radial receiver function (SV deconvolved from P) (Fig. 5a), the phase at ~100 km is most prominent (larger amplitudes over a larger range of back-azimuths) when negative. It is also “concave up” in the sense that the depth of the phase decreases with increasing back-azimuth over every 180˚. In contrast, the phase in the all-SH incidence radial receiver function (Fig. 5b) case appears most prominent when positive, and is “concave down.” The change in sign also occurs at different back-azimuths in the all-SH case relative to the all-SV case. In the all-SH case, the component used in the deconvolution (SV) is opposite of the incident (SH) particle motion, and hence represents energy that was wholly polarized from SH to SV by the ! 55! interaction of the wave with anisotropic structure. The pattern in back-azimuth for the case of equal parts SV and SH (Fig. 5c) is again very different, and in this case manifests 180˚ periodicity in which the negative and positive polarities occupy comparable ranges of back-azimuth. Bourguignon (2009) also tested cases of all-SV and all-SH incident polarization direction in Sp receiver functions, deconvolved with the radial component, and found systematic differences in the resulting receiver functions. The right-hand column in Figure 5 illustrates the cases for Model 3 in which the SH component is used in the deconvolution from P (transverse receiver functions). Each row again assumes a different incident particle motion: SV (Fig. 5b), SH (Fig. 5d) and equal SV and SH (Fig. 5f). In all three cases, the transverse receiver functions show systematic differences from their radial receiver function counterparts. In the case of incident SV energy, the phase at 100 km appears to be dominated by positive amplitudes (Fig. 5b), in contrast to the predominantly negative phase observed in the corresponding radial receiver function (Fig. 5a); the dominant phase on the transverse receiver function is also deeper. In the case of incident SH energy, the dominant negative polarity of the phase from 100 km on the transverse receiver function is also reversed from the corresponding (positive) radial receiver function phase although an argument could be made for a multi-lobed positive-negative-positive phase (Fig. 5d). In the all-SH case the dominant arrivals on the radial and transverse receiver functions are at comparable depths. In the case of the equal SH and SV particle motion, the transverse receiver function shows alternating negative and positive phases at ~100 km (as does the radial receiver function for this case) but the polarity flips occur at different back-azimuths. ! 56! Significant back-azimuthal variations in the depth and amplitude of the Moho phase exist in the all-SV and all-SH transverse receiver functions, and in the polarity of the Moho phase for the equal SH and SV transverse receiver function. The effect of the ratio of SV to SH incidence energy was also tested on the radially anisotropic models, Models 5 and 6. Because SH incident waves do not convert to P at the base of an anisotropic layer with this geometry of anisotropy, the cases with a pure SV incident polarization and mixed SV and SH incident polarization yield identical Sp receiver functions (Fig. 3). Because of this lack of dependence on initial particle motion, Sp conversions from the base of a layer of radial anisotropy with a slow vertical symmetry axis will in general produce a negative phase at all back-azimuths, and a layer of radial anisotropy with a fast vertical symmetry axis will produce a positive phase. 3.5 Strategies for addressing initial particle motion dependence Although the interplay between SH, SV and P at an anisotropic boundary has the potential to significantly complicate Sp receiver functions, it also produces unique patterns that could be useful in characterizing anisotropic structure. One caveat though, is the difficulty in measuring the amount of incident SV and SH particle motion in a real phase scattering at a boundary in anisotropy. Without the ability to measure this ratio directly beneath the discontinuity, our next best option is to approximate the initial polarization direction. One option would be to calculate the polarization pattern emitted from the moment tensor and to bin phases by their SV/SH ratio. Another option would be to employ a phase, such as SKS, which is naturally polarized SV. The limitation to ! 57! both methods however, is that any anisotropy along the path encountered prior to the receiver-side upper mantle discontinuity structure can change the ratio of SV to SH energy. In this paper, we explore the use of SKS phases to constrain initial polarization direction. To test the feasibility of using SKSp to image anisotropic structure, we generate synthetics for an epicentral distance range of 110°, which is within the SKS phase distance range. The results of the synthetic SKSp radial receiver function are shown in Figure 6b, which are compared to the corresponding Sp radial receiver function generated using the same model (Model 3) (Fig. 6a). Both Sp and SKSp are calculated using incident SV polarization. Generally, there is very good agreement between the Sp and SKSp receiver functions. Differences in the amplitude of the Moho and mantle phase are due to differences in incidence angle between the direct-S and SKS phases. As expected, the more vertically oriented SKS phase generates scattered phases with smaller amplitudes. 4. Data example 4.1 Craton structure and formation In global tomography models (e.g., Cammarano and Romanowicz, 2007; Kustowski et al., 2008; Lebedev and van der Hilst, 2008; Nettles and Dziewonski, 2008; Dalton et al., 2009; Romanowicz, 2009; Lekic and Romanowicz, 2011; Ritsema et al., 2011) the lithosphere of cratons is often found to be seismically faster and/or thicker than off-craton ! 58! continental or oceanic lithosphere. In a number of recent Sp receiver function studies of North America, little evidence is found for a phase associated with the boundary between the cratonic lithosphere and underlying asthenosphere (Kumar et al., 2012; Levander and Miller, 2012; Abt et al., 2010), although some studies have argued for a weak arrival from the lithosphere-asthenosphere boundary depth range (Miller and Eaton, 2010). Instead, the most prominent feature is a negative phase observed at depths thought to be within the mantle lithosphere. Evidence for discontinuity structure within the lithosphere has also been observed in Ps receiver functions (Bostock, 1998; Snyder, 2008; Mercier et al., 2008; Rychert and Shearer, 2009) as well as in long-range seismic profiles (Thybo, 2006). In some models of surface wave tomography, variations in azimuthal anisotropy at mid-lithospheric depths are suggestive of layered structure within the lithosphere (Yuan and Romanowicz, 2010; Yuan et al., 2011). These variations often fall within the same depth range as mid-lithospheric discontinuity phases imaged by Sp receiver functions (Abt et al., 2010). Such discontinuities have been interpreted to be a result of craton formation, which could either be explained by a more dehydrated, depleted cratonic core underlain by a less depleted mantle lithosphere (Yuan and Romanowicz, 2010) or by the stacking of lithospheric layers (Bostock, 1998; Mercier et al., 2008; Snyder, 2008). It has alternatively been suggested that the layer may reflect a melt cumulate resulting from melt migration or metasomatism (Ford et al., 2010). Despite these conjectures, the origin of the intra-lithospheric layering remains uncertain. ! 59! 4.2 Deconvolution and migration Building on the earlier work of Abt et al. (2010) we calculated single station Sp receiver functions using data from the Canadian National Seismograph Network (CN), POLARIS Network (PO) and the Global Seismograph Network (II), totaling 12 stations (Fig. 6). Events were requested for all back-azimuths, epicentral distances of 55-80° and depths of less than 300 km. As in Abt et al. (2010) and Ford et al. (2010), waveforms were filtered with a band-pass of 0.03 to 0.5 Hz, windowed around the S phase, and rotated into P and SV components using a free surface transfer matrix (Kennett, 1991). All the waveforms for a given station were then simultaneously deconvolved and migrated in the frequency domain, using a best fitting water-level to stabilize the deconvolution (Bostock, 1998). At each station, crustal parameters used in the migration were calculated using an H-k stacking method (Zhu and Kanamori, 2000). At stations SCHQ and EDM, the H-k stacking method failed to determine reasonable values of crustal thickness and/or Vp/Vs ratio, therefore a fixed model was used (Table 2). For mantle parameters, the migration model at each station was constructed by calculating 1D average profiles from the three dimensional models of Vp from Burdick et al. (2008) and Vs from Yuan and Romanowicz (2010) (Table 2). Station RES fell outside of the region modeled by Burdick et al. (2008), and as a result, the one-dimensional earth model, ak135 (Kennett et al., 1995) was used instead. 4.3 Results Moho depth, calculated using Ps receiver functions, was found to range from 36-44 km beneath the 12 stations where Sp receiver functions were also calculated. A negative ! 60! phase, observed on the Sp receiver functions, is found between 62 and 109 km depth (Table 2; Fig. 6), and corresponds to a velocity decrease with increasing depth, if assumed to be caused by isotropic changes in velocity. At each of these stations, the depth of the negative phase falls short of the potential lithosphere-asthenosphere depth range, which is broadly defined as the depth range from the first peak in shear wave velocity beneath the Moho, downwards to the largest minima in shear wave velocity, determined using the shear-wave velocity model of Yuan and Romanowicz (2010). As a result, we interpret the negative phases to be mid-lithospheric discontinuities. 4.4 Comparison to previously modeled receiver functions If we assume that the mid-lithospheric discontinuity is the result of an isotropic change in velocity, simple forward modeling for similar structure in Australia suggests that the negative phase may be the result of a thin low velocity layer, or the result of a rapid drop in velocity followed by a gradual increase (Ford et al., 2010). Another possibility is that the negative phase is a product of a boundary in anisotropy. Models 1 and 2 presented earlier in this paper produce robust negative phases in both the single-binned radial receiver function, and in back azimuth binned radial receiver functions (Fig. 2). Models 3 and 4 generate a negative phase in the single-binned receiver function, however polarity changes exist in the back azimuth binned receiver functions (Fig. 4). Model 6, a radial anisotropy model with a slow symmetry axis, also produces a negative phase in both the single- and back azimuth- binned transverse receiver function, and in this case no dependence on back-azimuth or incident polarization exists, suggesting radial anisotropy might be particularly easy to observe in Sp receiver functions (Fig. 3). Evidence of ! 61! vertical slow radial anisotropy in the mantle lithosphere beneath continents (e.g., Gaherty, 2004; Gaherty and Jordan, 1995; Nettles and Dziewonski, 2008; Yuan et al., 2011) further suggests radial anisotropy may be a plausible explanation. However, one caveat regarding radial anisotropy with a slow vertical symmetry axis is that creating such a layer through crystallographic alignment of olivine crystals would require deformation scenarios not typically considered plausible for the continental lithosphere (Rychert et al., 2007). Taken together, these modeling results suggest that the mid- lithospheric discontinuity observed in Sp receiver functions within the North American craton could be the result of a boundary in azimuthal or radial anisotropy. 4.5 Comparison to tomography results At several of the stations, the depth range over which negative phase energy is observed is significant and persists over depth ranges of more than 50 km (Fig. 6). Most interestingly, these broad regions of apparent velocity decrease often overlap with the change in fast-axis direction observed within the craton in the model of azimuthal anisotropy (Fig. 6), SAWum_NA2, constructed through the joint inversion of long-period surface waves and SKS splitting measurements (Yuan et al., 2011). This overlap suggests that the observed mid-lithospheric discontinuity in the Sp receiver functions may be the result of a vertical change in anisotropy. To test the hypothesis we calculated SKSp receiver functions, binned by back-azimuth, for station FCC and compare these results to a set of synthetic receiver functions produced from a simplified version of the SAWum_NA2 model for the same station. SKSp receiver functions are used in order to minimize incident SH particle motion, since, as we have shown, the ratio of incident SV ! 62! and SH particle motions has a significant effect on Sp receiver functions. 4.6 FCC SKSp receiver functions Station FCC (CN Network) is found along the western edge of Hudson Bay, well within cratonic North America (Fig. 6). A total of 1102 individual waveforms, from epicentral distances of 90-120°, were used. The rotation, deconvolution and migration operations performed on the observed waveforms were identical to those described in the previous section for the observed Sp receiver functions. Receiver functions were binned according to back-azimuth, in 5° increments (Fig. 8a). The distribution of events is limited to a relatively narrow back-azimuth range of 240- 360°. This range is large enough however, to observe variations in back-azimuth. The most notable feature, other than the presence of a Moho phase, is a negative phase at ~125 km depth whose apparent depth decreases from back-azimuths of ~250˚ to ~360˚. This phase is prominent in both the back-azimuth binned, and the single binned, SKSp receiver functions. A negative phase at ~55 km depth, and a positive phase at ~85 km depth are present as well. To determine the robustness of the SKSp receiver functions for station FCC, we bootstrapped each back-azimuth bin (Fig. 7b). The negative phase at ~125 km depth in the FCC SKSp receiver function is still a robust feature (negative even at two standard deviations) when bootstrapped. The negative and positive phases at depths of 55 and 85 ! 63! km are present at certain back-azimuths but are no longer significant features across all back-azimuth bins (Fig. 10). To determine whether or not the phase information contained in our SKSp receiver function for station FCC is related to anisotropic structure, we compute synthetic receiver functions for the same location, using information from the radially and azimuthally anisotropic model SAWum_NA2 (Yuan et al., 2011) to guide us in our model creation. 4.7 Model parameterization The SAWum_NA2 model (Yuan et al., 2011) is a model of three-dimensional variations in shear-wave velocity and anisotropy in the upper mantle beneath North America obtained by inverting a combination of long-period waveforms and SKS splitting measurements. Figure 9 (black line) illustrates one-dimensional variations in fast-axis azimuth (ψ), percent azimuthal anisotropy (G), radial anisotropy (ξ) and isotropic shear- wave velocity (Vs), from SAWum_NA2 beneath station FCC. To produce synthetic receiver functions, we discretized the SAWum_NA2 model for FCC into four mantle layers (including a half-space), and an isotropic crust (Table 3 and red lines in Figure 9). The interface between mantle layers 1 and 2 and between layers 2 and 3 correspond to steps in Vs and ξ, while the azimuthal anisotropy properties remain constant. Between layers 3 and 4, there is a change in azimuthal anisotropy but not in Vs or ξ. Three significant differences exist between SAWum_NA2 and its simplified parameterization used for the synthetic SKSp receiver function calculations. First, the ! 64! depths of the first two interfaces shown in our preferred model do not align with SAWum_NA2 gradients in Vs. Rather, the low velocity zone with a maximum at ~100 km in SAWum_NA2 is shifted to a shallower depth to better match the negative and positive phases observed in the FCC SKSp receiver functions at ~55 km and ~85 km, respectively, under the assumption that receiver functions can better constrain the depth of an interface than can tomography. Second, the percent of azimuthal anisotropy (G) was doubled, in agreement with the conclusion of Yuan et al. (2011) that SAWum_NA2 underestimates azimuthal anisotropy at mantle depths Third and finally, the drop in shear velocity and slight rotation of azimuthal anisotropy below 150 km depth in the SAWum_NA2 model are not included in our simplified model for FCC. These velocity gradients are thought to be representative of the lithosphere-asthenosphere boundary. Since no phase is observed in Sp receiver functions from the LAB depth range, either at FCC or in our results from cratons in general (Abt et al., 2010; Ford et al., 2010), velocity gradients in this depth range are likely to be as gradual in depth as they appear in the SAWum_NA2 model. Even if these gradual gradients were explicitly added to the simplified model, they would not significantly alter the predicted Sp receiver functions. Overall, our simplified FCC model should be thought of as a hybrid structure which utilizes information from both receiver function analysis and tomography. 4.8 Translating model parameters into elastic coefficients In order to translate the SAWum_NA2 model (as described by Vs, ξ, ψ and G) into the 6x6 tensor of elastic coefficients, Cij, used to generate our synthetic seismograms, we employ the empirical relationships assumed in the SAWum_NA2 inversion to scale ! 65! between Vs, ρ (density) and Vp, where Vs and Vp are the Voigt average isotropic ! !!" !! ! velocities, and to scale between radial anisotropy (! = ! !!" ), ϕ!(= !!" ! ), and η (= !!!! ) !" (Montagner and Anderson, 1989). ! !" !! ! !" !! = 0.5 (1) ! !" ! ! !" !! = 0.33 (2) ! !" ! ! !" ! = −2.5 (3) ! !" ! ! !" ! = −1.5 (4) These parameters were then used to solve for the Love coefficients A, C, F, L and N, which are independent of azimuth (Love, 1927). ! ! = !!!" (5) ! ! = !!!" (6) ! ! = !!!" (7) ! ! = !!!" (8) ! 66! ! = !(! − 2!) (9) The azimuthal anisotropy parameters in the SAWum_NA2 model (G and ψ) are derived from the 2θ azimuthal terms, GC and Gs, of Montagner and Nataf (1986). ! = ! !!! + !!! (10) !! = 0.5 tan!! (!! !! ) (11)! If the 4θ terms for azimuthal anisotropy are assumed to be zero, and the fast-axis direction, ψ, used in the 2θ terms, is set to 0°, Cij can be described as: ! + !" ! − 2! ! + !" !!!!!!! !!!!!!!!!! !!!!!0 ! − !" ! − !" !!!!!! !!!!!!!!!! !!!!!0 ! !!! !!!!!!!!!! !!!!!0 (12) ! − !" 0 ! ! + !" ! ! The other 2θ azimuthal terms in Montagner and Nataf (1986), BC,S and HC,S, were not solved for in the SAWum_NA2 inversion. Therefore, we obtained values for AB and HF by scaling them from the SAWum_NA2 LG values using the relationships between LG, AB and HF found in the Schulte-Pelkum and Blackman (2003) coefficients used earlier in this paper. However, any other reasonable set of coefficients could be used. The ! 67! coefficients marked as “X” were set to zero, since no information on these coefficients was obtainable (Montagner and Nataf, 1986). The complete set of elastic coefficients for our model of FCC can be seen in Table 3. 4.9 Computing and comparing synthetic SKSp receiver functions Similar to the simple synthetic cases computed earlier, we use the elastic coefficients from our FCC model to compute synthetic seismograms using a propagator-matrix method (Keith and Crampin, 1977). To model SKSp, we assume initially SV incident particle motion and use an epicentral distance of 110°, which falls within the range of SKS. The results of this synthetic example are shown in Figure 10. Beneath the Moho, the two isotropic changes in Vs, a decrease in Vs at 55 km, and an increase in Vs at 80 km, are only faintly perceptible, in large part due to interference with the Moho phase, which is located at 42 km depth and was constrained from H-k stacking of Ps receiver functions for station FCC. These phases are not as pronounced as the negative/positive phase pair observed in the SKSp receiver function for station FCC. The most striking feature of the synthetic receiver function is observed at a depth of ~125 km, which corresponds to the change in azimuthal anisotropy found between layers 3 and 4 (Table 3). Here, a negative phase is seen to steadily increase in amplitude before rapidly transitioning to a positive phase at back-azimuths of ~140° and ~320°. The back- azimuthal variation of this synthetic phase agrees reasonably well with the negative phase observed in the SKSp receiver function for FCC (Fig. 8). An interesting facet of this result is that it shows that the negative phase observed on the stack over all back- ! 68! azimuths at FCC (Fig. 7) can be explained by a model in which the discontinuity is produced by azimuthal anisotropy with no decrease in average isotropic velocity (Fig. 9). 4.10 Implications The coincidence in the back-azimuth patterns between the synthetic model and the calculated SKSp receiver functions for station FCC suggests that the boundary in anisotropy observed in SAWum_NA2, may be related to the pervasive mid-lithospheric discontinuity imaged in Sp receiver functions throughout cratonic North America (Kumar et al., 2012; Levander and Miller, 2012; Miller and Eaton, 2010; Abt et al., 2010). However, this initial conclusion will remain speculative until robust comparisons of the SAWum_NA2 and observed SKSp or Sp receiver functions are completed at other stations in the North American craton. In addition, the implications of SKSp and Sp receiver functions also need to be reconciled with anisotropic structures implied by Ps receiver functions (e.g., Bostock, 1998; Snyder, 2008; Mercier et al., 2008). Regardless, modeling of azimuthal anisotropy in SKSp and Sp receiver functions holds promise for better understanding the internal structure of the cratonic lithosphere and hence the mechanisms responsible for craton formation. 5. Conclusions Sp receiver functions depend on the orientation and strength of anisotropy, as has been observed in Ps receiver functions. An important difference, however, is that Sp receiver functions also vary significantly with the polarization of the incident S phase. For ! 69! researchers interested in interrogating observed Sp receiver functions for signs of anisotropy, SKSp receiver functions should minimize complications due to incident polarization direction, although use of the direct S ray parameter range may be feasible with careful binning as a function of incident polarization. Numerous stations in the Canadian craton (and globally) manifest a negative Sp phase at mid-lithospheric depths when data from different back-azimuths are averaged. We have found that back-azimuthally-averaged negative phases are produced by a variety of anisotropic models, including radial anisotropy with a vertical slow axis, and, with the right back-azimuth and incident polarization sampling, certain azimuthally anisotropic models. SKSp receiver functions observed at station FCC in the Canadian craton contain a negative phase whose depth varies from 140 km to 100 km as a function of back- azimuth. The SKSp phase is matched by synthetic SKSp phases generated by a rapid vertical change in the orientation and strength of azimuthal anisotropy, unaccompanied by a comparable gradient in isotropic velocity, as suggested by the SAWum_NA2 model. A goal for future analyses will be to test whether mid-lithospheric layering in azimuthal anisotropy can in general explain the negative mid-lithospheric discontinuity seen in Sp and SKSp receiver functions. Understanding this relationship should shed light on models for the formation of the cratonic mantle lithosphere. ! 70! References Abt, David L., et al. "North American lithospheric discontinuity structure imaged by Ps and Sp receiver functions." J. geophys. Res 115.B09301 (2010): B09301. Babuška, Vladislav, and Michel Cara. Seismic anisotropy in the Earth. Vol. 10. Springer, 1991. Bostock, M. G. "Mantle stratigraphy and evolution of the Slave province."Journal of Geophysical Research 103.B9 (1998): 21183-21. Burdick, Scott, et al. "Upper mantle heterogeneity beneath North America from travel time tomography with global and USArray transportable array data."Seismological Research Letters 79.3 (2008): 384-392. Bourguignon, Sandra. Lithospheric Deformation at the South Island Oblique Collision, New Zealand. Diss. Victoria University of Wellington, Wellington, New Zealand, 2009. Bystricky, M., et al. "High shear strain of olivine aggregates: rheological and seismic consequences." Science 290.5496 (2000): 1564-1567. Cammarano, Fabio, and Barbara Romanowicz. "Insights into the nature of the transition zone from physically constrained inversion of long-period seismic data." Proceedings of the National Academy of Sciences 104.22 (2007): 9139-9144. Dalton, Colleen A., Göran Ekström, and Adam M. Dziewonski. "Global seismological shear velocity and attenuation: A comparison with experimental observations." Earth and Planetary Science Letters 284.1 (2009): 65-75. Ford, Heather A., et al. "The lithosphere–asthenosphere boundary and cratonic lithospheric layering beneath Australia from Sp wave imaging." Earth and Planetary Science Letters 300.3 (2010): 299-310. Forsyth, Donald W., et al. "Phase velocities of Rayleigh waves in the MELT experiment on the East Pacific Rise." Science 280.5367 (1998): 1235-1238. Frederiksen, A. W., and M. G. Bostock. "Modelling teleseismic waves in dipping anisotropic structures." Geophysical Journal International 141.2 (2002): 401-412. Gaherty, James B. "A surface wave analysis of seismic anisotropy beneath eastern North America." Geophysical Journal International 158.3 (2004): 1053-1066. Gaherty, James B., and Thomas H. Jordan. "Lehmann discontinuity as the base of an anisotropic layer beneath continents." Science 268 (1995): 1468-1471. Keith, C.M., Crampin, S., 1977. Seismic body waves in anisotropic media: synthetic seismograms. Geophys. J. R. Astron. Soc. 49, 225–243. Kennett, B. L. N., E. R. Engdahl, and R. Buland. "Constraints on seismic velocities in the Earth from traveltimes." Geophysical Journal International122.1 (1995): 108-124. Kennett, B. L. N. "The Removal of Free Surface Interactions From Three‐Component Seismograms." Geophysical journal international 104.1 (2007): 153-154. Kind, Rainer, Xiaohui Yuan, and Prakash Kumar. "Seismic receiver functions and the lithosphere-asthenosphere boundary." Tectonophysics (2012). Kumar, P., et al. "The lithosphere-asthenosphere boundary observed with USArray receiver functions." Solid Earth 3 (2012): 149-159. Kustowski, B., G. Ekström, and A. M. Dziewoński. "Anisotropic shear-wave velocity structure of the Earth's mantle: a global model." Journal of Geophysical Research 113.B6 (2008): B06306. ! 71! Lebedev, Sergei, and Rob D. Van Der Hilst. "Global upper‐mantle tomography with the automated multimode inversion of surface and S‐wave forms."Geophysical Journal International 173.2 (2008): 505-518. Lekić, V., and B. Romanowicz. "Inferring upper‐mantle structure by full waveform tomography with the spectral element method." Geophysical Journal International 185.2 (2011): 799-831. Lekic, Vedran, Scott W. French, and Karen M. Fischer. "Lithospheric thinning beneath rifted regions of southern California." Science 334.6057 (2011): 783-787. Levander, Alan, and Meghan S. Miller. "Evolutionary aspects of lithosphere discontinuity structure in the western US." Geochemistry Geophysics Geosystems 13.null (2012): Q0AK07. Levin, Vadim, and Jeffrey Park. "P‐SH conversions in a flat‐layered medium with anisotropy of arbitrary orientation." Geophysical Journal International 131.2 (1997): 253-266. Levin, Vadim, and Jeffrey Park. "P-SH conversions in layered media with hexagonally symmetric anisotropy: a cookbook." Pure and applied geophysics151.2-4 (1998): 669-697. Li, Xueqing, Xiaohui Yuan, and Rainer Kind. "The lithosphere–asthenosphere boundary beneath the western United States." Geophysical Journal International 170.2 (2007): 700-710. Lin, Fan-Chi, et al. "Complex and variable crustal and uppermost mantle seismic anisotropy in the western United States." Nature Geoscience 4.1 (2010): 55-61. Love, A. E. H. The mathematical theory of elasticity. 1927. MacDougall, Julia G., Karen M. Fischer, and Megan L. Anderson. "Seismic anisotropy above and below the subducting Nazca lithosphere in southern South America." Journal of Geophysical Research: Solid Earth (1978–2012)117.B12 (2012). Mainprice, David, and Paul G. Silver. "Interpretation of SKS-waves using samples from the subcontinental lithosphere." Physics of the Earth and Planetary Interiors 78.3 (1993): 257-280. Mehl, Luc, et al. "Arc-parallel flow within the mantle wedge: Evidence from the accreted Talkeetna arc, south central Alaska." Journal of geophysical research108.B8 (2003): 2375. Mercier, J. P., et al. "The teleseismic signature of fossil subduction: Northwestern Canada." J. geophys. Res 113 (2008): B04308. Michibayashi, Katsuyoshi, et al. "Seismic anisotropy in the uppermost mantle, back-arc region of the northeast Japan arc: Petrophysical analyses of Ichinomegata peridotite xenoliths." Geophysical research letters 33.10 (2006): L10312. Miller, Meghan S., and David W. Eaton. "Formation of cratonic mantle keels by arc accretion: evidence from S receiver functions." Geophysical Research Letters 37.18 (2010): L18305. Montagner, Jean-Paul, and Don L. Anderson. "Petrological constraints on seismic anisotropy." Physics of the earth and planetary interiors 54.1 (1989): 82-105. Montagner, Jean-Paul, and Henri-Claude Nataf. "A simple method for inverting the azimuthal anisotropy of surface waves." Journal of Geophysical Research91.B1 (1986): 511-520. ! 72! Nettles, Meredith, and Adam M. Dziewoński. "Radially anisotropic shear velocity structure of the upper mantle globally and beneath North America."Journal of Geophysical Research 113.B2 (2008): B02303. Ozacar, A. Arda, and George Zandt. "Crustal structure and seismic anisotropy near the San Andreas Fault at Parkfield, California." Geophysical Journal International 178.2 (2009): 1098-1104. Ritsema, J., et al. "S40RTS: a degree‐40 shear‐velocity model for the mantle from new Rayleigh wave dispersion, teleseismic traveltime and normal‐mode splitting function measurements." Geophysical Journal International 184.3 (2011): 1223-1236. Romanowicz, Barbara. "The thickness of tectonic plates." Science 324.5926 (2009): 474- 476. Rychert, Catherine A., and Peter M. Shearer. "A global view of the lithosphere- asthenosphere boundary." Science 324.5926 (2009): 495-498. Rychert, Catherine A., Stéphane Rondenay, and Karen M. Fischer. "P-to-S and S-to-P imaging of a sharp lithosphere-asthenosphere boundary beneath eastern North America." Journal of Geophysical Research 112.B8 (2007): B08314. Savage, Martha K. "Lower crustal anisotropy or dipping boundaries? Effects on receiver functions and a case study in New Zealand." Journal of geophysical research 103.B7 (1998): 15069-15. Savage, M. K. "Seismic anisotropy and mantle deformation: what have we learned from shear wave splitting." Rev. Geophys 37.1 (1999): 65-106. Schulte‐Pelkum, Vera, and Donna K. Blackman. "A synthesis of seismic P and S anisotropy." Geophysical Journal International 154.1 (2003): 166-178. Schulte-Pelkum, Vera, et al. "Imaging the Indian subcontinent beneath the Himalaya." Nature 435.7046 (2005): 1222-1225. Silver, Paul G., and W. Winston Chan. "Shear wave splitting and subcontinental mantle deformation." J. geophys. Res 96.16 (1991): 429-16. Silver, Paul G. "Seismic anisotropy beneath the continents: probing the depths of geology." Annual Review of Earth and Planetary Sciences 24 (1996): 385-432. Smith, Gideon P., and Göran Ekström. "A global study of Pn anisotropy beneath continents." Journal of geophysical research 104.B1 (1999): 963-980. Snyder, David B. "Stacked uppermost mantle layers within the Slave craton of NW Canada as defined by anisotropic seismic discontinuities." Tectonics 27.4 (2008): TC4006. Thybo, Hans. "The heterogeneous upper mantle low velocity zone."Tectonophysics 416.1 (2006): 53-79. Wirth, Erin A., and Maureen D. Long. "Multiple layers of seismic anisotropy and a low- velocity region in the mantle wedge beneath Japan: Evidence from teleseismic receiver functions." Geochemistry Geophysics Geosystems 13.null (2012): Q08005. Yang, Yingjie, and Donald W. Forsyth. "Rayleigh wave phase velocities, small-scale convection, and azimuthal anisotropy beneath southern California." J. geophys. Res 111 (2006): B07306. Yuan, Huaiyu, et al. "3‐D shear wave radially and azimuthally anisotropic velocity model of the North American upper mantle." Geophysical Journal International 184.3 (2011): 1237-1260. ! 73! Yuan, Huaiyu, and Barbara Romanowicz. "Lithospheric layering in the North American craton." Nature 466.7310 (2010): 1063-1068. Vinnik, L. P. "Detection of waves converted from P to SV in the mantle."Physics of the Earth and Planetary Interiors 15.1 (1977): 39-45. Zhu, Lupei, and Hiroo Kanamori. "Moho depth variation in southern California from teleseismic receiver functions." Journal of Geophysical Research 105 (2000): 2969- 2980. ! 74! Model 1 Model 2 Model 3 Model 4 N-S N-S N-S N-S x x x SPB-Half x SPB-Full SPB-Half SPB-Half dip 20N E-W E-W NE-SW NE-SW SPB-Full SPB-Half SPB-Half SPB-Half Model 5 Model 6 SPB-Half SPB-Slow SP1-Iso SP2-Iso Figure 1. Schematic illustrations of the models used to generate synthetic receiver functions. Model numbers correspond to those listed in Table 1. ! 75! ! a. b. Back Azimuth (degrees) Back Azimuth (degrees) Depth (km) Depth (km) c. d. Back Azimuth (degrees) Back Azimuth (degrees) Depth (km) Depth (km) Amplitude Figure 2. Comparison of Ps and Sp back-azimuth-binned & single-binned synthetic receiver functions. (a) and (b) are Ps and Sp receiver functions, respectively, calculated using Model 1 (Table 1), and (c) and (d) are Ps and Sp receiver functions calculated using Model 2 (Table 1). For the Sp cases, incident particle motion is SV, and the receiver functions correspond to SV deconvolved from P. The positive (red) phase at a depth of 36 km in (a)-(d) is the Moho. The negative (blue) phase observed at a depth of ~100 km in (b) and (d) is the boundary between the two anisotropic layers. The corresponding phase in the Ps receiver functions ((a) and (c)) changes sign as a function of back- azimuth. The amplitude of the phase is diminished in (c) relative to (a). The positive and negative phases observed in (a) and (c), at a depth of 150 km and 190 km, respectively, are crustal reverberations. The discordances observed at 75-90° and 270-285° in (b) and (d) are likely to due to interference between the Moho phase (red) and the negative phase at 100 km. ! 76! a. b. 0.1 0.1 0 0 350 350 300 300 Back Azimuth (degrees) Back Azimuth (degrees) 150 150 100 100 50 50 0 0 0 50 100 150 0 50 100 150 Depth (km) Depth (km) c. d. 0.1 0.1 0 0 350 350 300 300 Back Azimuth (degrees) Back Azimuth (degrees) 150 150 100 100 50 50 0 0 0 50 100 150 0 50 100 150 Depth (km) Depth (km) 0 0.1 Amplitude Figure 3. Variations in varying types of radial anisotropy in Sp receiver functions, with varying amounts of SH energy. In (a) the symmetry axis is fast and oriented vertically; the input energy is SV only. The resulting phase that is produced at the boundary in anisotropy is positive. In (b) the symmetry axis is slow, and oriented vertically; the input energy is SV only. The resulting phase is negative. The elastic coefficients and symmetry are identical between (c) and (a) and (d) and (b). In (c) and (d), the input energy is equal parts SV and SH, indicating SH has little effect on the resulting receiver functions. ! 77! a. b. 0.1 0.1 0 0 350 350 300 300 Back Azimuth (degrees) Back Azimuth (degrees) 150 150 100 100 50 50 0 0 0 50 100 150 0 50 100 150 Depth (km) Depth (km) c. d. 0.1 0.1 0 0 350 350 300 300 Back Azimuth (degrees) Back Azimuth (degrees) 150 150 100 100 50 50 0 0 0 50 100 150 0 50 100 150 Depth (km) Depth (km) 0 0.1 Amplitude Figure 4. Comparison of effects of variations in fast-axis direction on synthetic Ps and Sp receiver functions. In (a)-(d), the difference in the horizontal offset in fast-axis between the two layers of anisotropy is 45° (Model 3 & 4). For the Sp cases, incident particle motion is SV, and the receiver functions correspond to SV deconvolved from P. The Ps receiver function in (a) shows a pattern of behavior similar to Figure 1a & 1c, except that the minimum and maximum amplitudes for the phase at ~100 km is offset relative to the earlier case. The pattern in back-azimuth is significantly different for the Sp receiver function (b). For cases in which the top layer has a fast axis which plunges at 20° to the N (0˚)(Model 4), the 180° periodicity in back-azimuth no longer exists for either (c) or (d). The offset in Moho depth and the large negative phase directly beneath the Moho in (b) and (d) are likely due to interference between the Moho phase and the energy associated with the boundary in anisotropy. ! 78! a. 0.1 0 b. 0.05 0 350 350 300 300 Back Azimuth (degrees) Back Azimuth (degrees) 150 150 100 100 50 50 0 0 0 50 100 150 0 50 100 150 Depth (km) c. 0.1 0 d. 0.05 0 350 350 300 300 Back Azimuth (degrees) Back Azimuth (degrees) 150 150 100 100 50 50 0 0 0 50 100 150 0 50 100 150 Depth (km) Depth (km) 0.05 e. 0 f. 0.01 0 350 350 300 300 Back Azimuth (degrees) Back Azimuth (degrees) 150 150 100 100 50 50 0 0 0 50 100 150 0 50 100 150 Depth (km) Depth (km) 0 0.1 Amplitude Figure 5. Synthetic Sp receiver functions calculated using a range of incident SV to SH polarizations, and varying the component (SV or SH) used in the deconvolution. In (a), (c) and (e) the SV component is deconvolved from P (radial). In (b), (d), and (f), the SH component is deconvolved from P (transverse). The incident polarization direction is all- SV motion in (a) and (b). In (c) and (d), the incident polarization direction is all-SH, while in (e) and (f), the amounts of incident SV to SH energy are equal. All six Sp receiver functions were calculated using Model 3. The offset in Moho depth and the large negative phase directly beneath the Moho in (a) and (c) are likely due to interference between the Moho phase and the energy associated with the boundary in anisotropy. ! 79! a. b. 0.1 0 0 350 350 300 300 Back Azimuth (degrees) Back Azimuth (degrees) 150 150 100 100 50 50 0 0 0 50 100 150 0 50 100 150 Depth (km) Depth (km) 0 0.1 Amplitude Figure 6. Synthetic Sp (a) and SKSp (b) receiver functions calculated for Model 3 (Table 1). Overall patterns in back-azimuth appear to be similar between (a) and (b), however the depth range of the phase at ~100 km appears to be slightly broader for the SKSp receiver functions (b). The amplitudes in the single binned receiver function are also larger in (a) than in (b), owing to the difference in incidence angle. Incident particle motion is SV, and the receiver functions correspond to SV deconvolved from P. ! 80! Distance (km) 0 500 1000 1500 2000 2500 3000 3500 4000 a. b. A A’ EYMN YKW3 UTMT JFWS JERN PLAL SIUC ULM FCC SLM FFC 0 RES Depth (km) 50 A JERN FRB YKW3 100 IVKQ FCC SCHQ 150 FFC EDM ULM 200 SADO 250 H | Su pe rior c. A | R a e /H e a r ne | T |M c R | Ya va /M A’ a ve aza ma y| S l A’ Wo p |G 0 100 50 60 70 80 90 100 110 120 m) 200 th (k Depth (km) 300 Dep 55° 50° 45° 60° 65° Latitude 40° 400 35° Deviation from APM 0º 50º Figure 7. (a) Gives an overview of lithosphere-asthenosphere boundary and mid- lithospheric discontinuity depths from Sp receiver functions in North America (modified from Abt et al., 2010). Black inverted triangles indicate that the negative corresponds to the LAB, while a white triangle indicates that the negative phase is interpreted to be a MLD. Grey triangles indicate uncertainty in the interpretation. (b) Profile through Sp receiver functions shown as line A-A’ in (a). Positive phases correspond to a velocity increase with depth, and negative phases correspond to an apparent velocity decrease with depth. Magenta line marks the depth pick of the Moho, and the cyan box indicates the pick (± 2σ) of the LAB or MLD. Phase interpretation is dependent on the potential LAB depth range indicated by the gray box (Yuan and Romanowicz, 2010). Phases found at depths shallower than the potential LAB depth range are interpreted to be a MLD. (c) Model of azimuthal anisotropy from within the North American craton from Yuan and Romanowicz (2010). The change in fast-axis direction observed at depths of ~80-150 km, which is thought to be within the cratonic lithosphere, often occur in a similar depth range to negative phase energy associated with MLD phases observed in (b). ! 81! 0 0 50 50 100 100 Depth (km) Depth (km) 150 150 200 200 250 250 0 50 100 150 200 250 300 350 150 175 200 225 250 275 300 325 350 Back Azimuth (degrees) Back Azimuth (degrees) Figure 8. (a) Data example of back-azimuth binned SKSp receiver functions from station FCC. A strong, consistent negative phase is observed at a depth of ~130-150 km, and appears to either diminish in amplitude and/or flip in polarity at 310 °. The Moho phase, characterized by positive energy, is observed at a depth of approximately 40 km. Additional negative and positive phases, found at ~55 km and ~85 km depth respectively, are also consistent across back-azimuth bins. (b) Plot of bootstrapped, SKSp receiver functions for station FCC. Only the portions of the receiver function that can be resolved (± 2σ) are shown. The negative phase at depths of ~125 km, thought to be associated with a boundary in anisotropy, is well resolved ! 82! Azimuth Vs (m/s) Xi G depth (km) 1 1 2 3 Figure 9. Model of (a) azimuthal anisotropy fast-axis direction, (b) isotropic shear-wave velocities, (c) radial anisotropy and (d) strength of anisotropy, shown for station FCC. Solid black lines in (a)-(d) indicate model values from Yuan et al. (2011). Red lines in (a)-(d) indicate the simplified model used in this study (Table 3). ! 83! Back Azimuth (degrees) Depth (km) Figure 10. Synthetic SKSp receiver functions generated from our simplified model (Table 3) for station FCC. A rapid change in polarity is observed at ~130 km depth and at back- azimuths of 140° and 320°, with a rapid shallowing of the negative phase prior to become positive. Incident particle motion is SV, and the receiver functions correspond to SV deconvolved from P. ! 84! SPB - SPB - SPB - i j crust full half slow SP1 - iso SP2 - iso 1 1 121.0 241.3 226.8 212.7 212.2 227.8 2 2 121.0 197.6 204.9 234.6 212.2 227.8 3 3 121.0 197.6 204.9 234.6 212.2 227.8 4 4 39.9 60.5 65.0 69.5 69.4 65.0 5 5 39.9 73.9 71.7 62.8 69.4 65.0 6 6 39.9 73.9 71.7 62.8 69.4 65.0 1 2 41.2 71.7 72.5 98.2 73.3 97.8 1 3 41.2 71.7 72.5 98.2 73.3 97.8 1 4 0.0 0.0 0.0 0.0 0.0 0.0 1 5 0.0 0.0 0.0 0.0 0.0 0.0 1 6 0.0 0.0 0.0 0.0 0.0 0.0 2 3 41.2 76.6 75.0 95.7 73.3 97.8 2 4 0.0 0.0 0.0 0.0 0.0 0.0 2 5 0.0 0.0 0.0 0.0 0.0 0.0 2 6 0.0 0.0 0.0 0.0 0.0 0.0 3 4 0.0 0.0 0.0 0.0 0.0 0.0 3 5 0.0 0.0 0.0 0.0 0.0 0.0 3 6 0.0 0.0 0.0 0.0 0.0 0.0 4 5 0.0 0.0 0.0 0.0 0.0 0.0 4 6 0.0 0.0 0.0 0.0 0.0 0.0 5 6 0.0 0.0 0.0 0.0 0.0 0.0 Table 1a. List of elastic coefficients used in models 1-6 (Table 1b, Fig. 1), in the form of the 6x6 elastic tensor Cij. Units sre in GPa. ! 85! mantle half- crust mantle layer 1 space H H Model coeff (km) coeff (km) az dip coeff az 1 crust 36 SPB - full 64 0 0 SPB - full 90 SPB - 2 crust 36 half 64 0 0 SPB - half 90 SPB - 3 crust 36 half 64 0 0 SPB - half 45 SPB - 4 crust 36 half 64 0 20 SPB - half 45 SPB - 5 crust 36 half 64 NA 90 SP1 - iso NA SPB - 6 crust 36 slow 64 NA 90 SP2 - iso NA Table 1b. Models used to calculate synthetic seismograms. See Figure 1 for cartoon of models. Table 1a for corresponding elastic coefficients. ! 86! Velocity Model Phase NET STA Crust Moho Vp Vp/Vs Mantle Vp Mantle Vs MLD Error # RFs CN EDM Fixed 39 6.47 1.75 MIT_NA_Vp UCB_NA_Vs 85 11 130 CN FCC H-k Stack 42 6.46 1.72 MIT_NA_Vp UCB_NA_Vs 91 38 292 CN FRB H-k Stack 42 6.42 1.78 MIT_NA_Vp UCB_NA_Vs 65 7 299 CN RES AK135 NA NA NA AK135 AK135 62 36 461 CN SADO H-k Stack 39 6.49 1.77 MIT_NA_Vp UCB_NA_Vs 101 8 116 CN SCHQ Fixed 24 6.5 1.77 MIT_NA_Vp UCB_NA_Vs 109 13 404 CN ULM H-k Stack 34 6.42 1.73 MIT_NA_Vp UCB_NA_Vs 105 12 83 CN YKW3 H-k Stack 36 6.46 1.74 MIT_NA_Vp UCB_NA_Vs 79 18 261 II FFC H-k Stack 41 6.52 1.72 MIT_NA_Vp UCB_NA_Vs 89 18 271 PO IVKQ H-k Stack 35 6.42 1.77 MIT_NA_Vp UCB_NA_Vs 76 19 60 PO JERN H-k Stack 37 6.46 1.71 MIT_NA_Vp UCB_NA_Vs 72 8 72 PO SILO H-k Stack 38 6.46 1.73 MIT_NA_Vp UCB_NA_Vs 92 25 85 Table 2. Summary of receiver function results for the twelve stations used in the study. See text for discription of migration models. The depth and uncertaintly of the mid-lithospheric discontinuity phase was determined by selecting the largest, statisically significant phase. Error bars are two standard deviations calculated using bootstrap analysis ! 87! i j crust layer 1 layer 2 layer 3 layer 4 (hs) 1 1 121.000 246.214 247.619 264.657 258.877 2 121.000 232.486 234.181 250.503 256.283 3 3 121.000 223.420 220.490 228.380 228.380 4 4 39.916 74.183 74.183 76.481 78.254 5 5 39.916 78.395 78.395 80.823 79.050 6 6 39.916 79.410 79.410 84.230 84.230 1 2 41.168 80.530 82.798 89.120 89.120 1 3 41.168 76.777 80.391 86.274 85.618 1 4 41.168 0.000 0.000 0.000 0.000 1 5 0.000 0.000 0.000 0.000 0.000 1 6 0.000 0.000 0.000 0.000 0.000 2 3 0.000 75.219 78.865 85.324 85.324 2 4 0.000 0.000 0.000 0.000 0.000 2 5 0.000 0.000 0.000 0.000 0.000 2 6 0.000 0.000 0.000 0.000 0.000 3 4 0.000 0.000 0.000 0.000 0.000 3 5 0.000 0.000 0.000 0.000 0.000 3 6 0.000 0.000 0.000 0.000 0.000 4 5 0.000 0.000 0.000 0.000 0.000 4 6 0.000 0.000 0.000 0.000 0.000 5 6 0.000 0.000 0.000 0.000 0.000 AZ (°) NA 80 80 80 -42 G (%) NA 2.76 2.76 2.76 0.51 Vs (m/s) 3760 4770 4715.5 4807 4807 H km 42 13 25 45 NA ξ NA 1.0409 1.0587 1.0709 1.0709 Table 3. List of elastic coefficients in the form of the elastic tensor Cij. Units are in GPa. The corresponding azimuth (AZ), percent azimuthal anisotropy (G), shear wave velocity (Vs), layer thickness (H) and radial anisotropy (ξ) for each layer are listed as well. ! 88! Chapter 3 A comparison of the seismological and rheological lithosphere-asthenosphere boundary beneath cratonic North America By: Heather A. Ford1, Greg Hirth1 and Karen M. Fischer1 1 Department of Geological Sciences, Brown University ! 89! Abstract From a rheological point of view, the lithosphere-asthenosphere boundary is an upper mantle discontinuity characterized by a rapid change in viscosities, which allows for differential movement to occur between the overriding lithosphere and the convecting mantle asthenosphere. Seismic imaging has shown that a drop in velocity often exists at depths similar to where the rheologically defined lithosphere asthenosphere boundary is thought to be located. Whether the seismically imaged lithosphere-asthenosphere boundary is the same as the rheological boundary is yet unclear. Using a number of global and regional seismic tomography models for the Slave and Superior cratons in North America, we translated shear velocity into temperature, and then to viscosity. We find that when velocities are translated into temperature, the misfit between our results and xenolith-constrained geotherms are significant. Including a compositional correction for Mg# (which reflects depletion) and volume percent garnet helps to reduce the misfit in temperatures (and thereby viscosity) significantly. While our goal was to constrain the rheology of the lithosphere-asthenosphere boundary using velocity, the large discrepancies encountered in temperature suggest that the fidelity of seismic models needs to be addressed. In assuming a dislocation creep flow law with a strain rate of 10-15 1/s, we find that the average maximum viscosity of the mantle lithosphere in the Slave craton is approximately 5x1024 Pa*s, and 3x1024 Pa*s beneath the Superior craton. Minimum asthenospheric viscosity is 3x1021 Pa*s and 6x1021 Pa*s beneath the Slave and Superior ! 90! cratons. The value of 1021 Pa*s agrees well with estimates of asthenospheric viscosities from glacial glacial isostatic adjustment. The contrast between lithospheric and asthenospheric viscosities implied by our scaling of seismic velocity models suggests that the seismologically defined lithosphere-asthenosphere boundary corresponds to the rheological lithosphere-asthenosphere boundary. Interestingly, the average value of viscosity observed in the sub-cratonic asthenosphere is two orders of magnitude larger than what is observed in the western U.S., which implies that the asthenosphere does not maintain a constant viscosity everywhere. ! 91! 1. Introduction The lithosphere is commonly described as a rigid plate, which translates across the top of the convecting, lower-viscosity, asthenosphere. As such, the boundary between the two is typically defined rheologically. The contrast in physical properties at the lithosphere- asthenosphere boundary, and their effects on the rheology of the lithospheric plate and convecting mantle are not as well understood. Arguments have been made for contributions from composition (changes in fertility and hydration), temperature, grain size and melt content (e.g. Hirth and Kohlstedt, 1996; Karato and Jung, 1998; Fischer et al., 2010). The lithosphere is typically characterized as being a region of high seismic velocities in global tomography models (e.g., Cammarano and Romanowicz, 2007; Kustowski et al., 2008; Lebedev and van der Hilst, 2008; Nettles and Dziewonski, 2008; Dalton et al., 2009; Romanowicz, 2009; Lekic and Romanowicz, 2011; Ritsema et al., 2011) and is often found to be thicker, and faster, beneath cratons than in other tectonic settings. In contrast to the fast velocities of the lithosphere, the asthenosphere is thought to coincide with the location of the low velocity zone imaged in global tomography models. Using a variety of seismic methods, including receiver function analysis, researchers can also image and constrain the velocity gradient at the lithosphere-asthenosphere boundary. This information can then be used to differentiate between the physical properties of the lithosphere and the asthenosphere (Fischer et al., 2010 and references therein). ! 92! Improved understanding of the lithosphere-asthenosphere boundary is an important science goal, and is one of the scientific targets of the NSF-funded Earthscope Program (Williams et al., 2010). However, it is unclear whether the seismically imaged lithosphere and asthenosphere reflect the rheologically defined lithosphere and asthenosphere. In this paper, we interrogate this question by using shear wave velocity models of the upper mantle to obtain viscosity profiles of the lithosphere and asthenosphere. In the process, we consider the effects of attenuation and composition on velocity, temperature and viscosity. We limit our discussion to the Slave and Superior cratons of North America, where global and regional velocity models show that the lithosphere is considerably faster and thicker than in regions off the craton (Nettles and Dziewonski, 2008; Bedle and van der Lee, 2009; Yuan et al., 2011; French et al., 2012). The rest of this paper is broken into four sections, in which we discuss the shear wave velocity models and corrections for attenuation (1), translating velocity into temperature and the necessary compositional corrections (2), using dislocation creep flow laws to translate temperature into viscosity (3), and the implications of our results (4). 2. Shear wave velocity profiles 2.1 Shear wave velocity profiles We compare five shear wave velocity models that sample the Slave and Superior cratons of North America (Fig. 1). These models include two global-scale models, SEMum2 ! 93! (French et al., 2012) and ND08 (Nettles and Dziewonski, 2008); two North American models, SAWum_NA2 (Yuan et al., 2011) and NA07 (Bedle and Van der Lee, 2009); and one model specific to the Slave craton (Chen et al., 2007). The global tomography models SEMum2 and ND08 use long period waveforms, and forward modeling using the spectral element method (SEM), to invert for shear velocities and radial anisotropy, and in the case of SEMum2, azimuthal anisotropy. SAWum_NA2 uses a data set of long period waveforms and SKS splitting to produce a model of shear velocity that incorporates variations in both radial and azimuthal anisotropy. Chen et al. (2007) account for variations in phase velocity and azimuthal anisotropy by inverting Rayleigh waves using a two-plane wave approach (Forsyth and Li, 2005), while NA07 is based on an inversion of body-S waves and Rayleigh waves for isotropic variations in shear wave velocities. For all of the models, excluding the model of Chen et al. (2007), we compute average regional 1D shear velocity profiles within, or immediately surrounding, the Slave and Superior cratons (Fig. 1). This averaging is done to minimize small-scale variations that while in some cases are resolvable, are not necessary representative of the craton as a whole. In general, all of the models exhibit evidence for a relatively thick, seismically fast, lithospheric lid, characteristic of many tomographically imaged cratons (Fig. 2). There is considerable variation in the averaged velocities for the Slave and Superior cratons and neither craton is consistently faster or slower than the other, when compared across all ! 94! models. Within the Slave craton, there is good agreement between SAWum_NA2 and SEMum2. The model of Chen et al. (2007) also agrees with SAWum_NA2 and SEMum2 at certain depth ranges, but is distinguished by a very pronounced low velocity zone at ~100 km depth. The amplitude of the low velocity zone, however, is not considered to be a robust feature of the model (Chen et al., 2007). There is also reasonable agreement between NA07 and ND08 within the Slave craton, except at shallow depths. In the Superior craton, there are larger differences between NA07 and ND08, although the agreement between these two models is better than their agreement with the SAWum_NA2 and SEMum2 models. SAWum_NA2 and SEMum2 still show considerable consistency within the Superior craton, however some variation between them does exist at depths greater than 225 km. Between the model sets, two distinct groups of profiles emerge. The first group, composed of NA07 and ND08, has the largest wave speeds at shallow depths, followed by a quasi-linear decrease in velocities down to the low-velocity zone found at approximately 150-225 km. The second group, which includes SAum_NA2, SEMum2, and Chen, exhibits a trend of increasing velocities from 50 to 150 km, reaching a peak at 150 km, before eventually decreasing in velocities from 150 km to 250 km depth. Differences between these tomographic models exist for a number of reasons, including differences in the amount, distribution and type of data used, damping and variations in starting model. SEMum2 and SAWum_NA2 both use starting models that depend on a mineralogically meaningful 1D velocity model (Cammarano et al., 2005), while ND08, ! 95! NA07 and Chen et al. (2007) use variations of the 1D velocity models PREM (Dziewonski et al., 1981), MC35 (van der Lee and Nolet, 1997) and ak135 (Kennett et al., 1995), respectively. Treatment of the crust can significantly impact results at shallow depths, and may explain the anomalously large values observed in model ND08 at depths of less than 75 km in the Superior craton (Fig. 2). The differences in the velocity profiles, particularly between the two groups, is largest at 150 km depths, where the first group (ND08 and NA07) has an average velocity that is 4% smaller than the average velocity of the second group (Chen, SEMum2, and SAWum_NA2). The discrepancies in velocity between the two sets of models are significant, and have important implications for our interpretation of the physical properties of the mantle beneath cratonic North America. However, we do not attempt to reconcile the differences between the two sets of models in this paper. Instead, we use these differences to understand the range of possible viscosities within the mantle lithosphere, understanding that future work should be done to better ascertain the correct velocity structure. It should also be noted that while the variations in velocity between the models are not trivial, they are small when compared to the average change in velocity between cratons and tectonically active regions, such as the western U.S. (e.g., Nettles and Dziewonski, 2008). 2.2 Attenuation effects ! 96! The effect of seismic attenuation on seismic velocities has been well documented and is an important feature at upper mantle depths (e.g. Dalton et al., 2008). To account for effects of attenuation we correct the shear velocity profiles (Fig. 2) using the following relationship, taken from Karato and Jung (1998): (1) where V0(T,P) is the elastic seismic velocity, dependent on temperature (T) and pressure (P). Q-1(ω,T,P,C) is seismic attenuation, and is dependent on the frequency (ω), temperature, pressure and water content (C). Although water content can also have an effect on elastic wave speeds (e.g., Jacobsen et al., 2008), the effect is small at low weight percent water, and is not considered in this paper. The cratonic lithosphere is thought to be relatively dry on the basis of conductivity profiles (Hirth et al., 2000; Evans et al., 2011), although the evidence from xenoliths is more equivocal (Lee et al., 2011). Attenuation effects are corrected for in the velocity profiles by using a Precambrian shield averaged Q-1 taken from Dalton et al. (2008). However, a concern with this model is that it attains Q values as low as ~100 within the depth range of the cratonic lithosphere. An assessment of other shear attenuation models suggests that Q ≈ 230 - 520 (Q-1 = 0.0043 – 0.0019) within the mantle lithosphere (Resovsky et al., 2005), while surface wave attenuation estimates find that Q = 200 (Q-1 = 0.005) within young oceanic lithosphere, increasing to a range of Q = 200-1100 (Q-1 = 0.005-0.9*10-3) beneath old oceanic lithosphere (Yang et al., 2007). These studies suggest that Q is at least 200 ! 97! within the lithosphere. Therefore, we do not allow Q-1 to exceed 0.005 (Q less than 200) in our model. Additionally, we assume α = 0.25. This assumption is consistent with seismic attenuation models which estimate α = 0.2-0.3 in the upper mantle low-velocity zone (Gaherty et al., 1999; Romanowicz, 1995) and α = 0.1-0.3 (Karato and Spetzler, 1990). The attenuation-corrected shear velocity profiles for the Slave and Superior cratons are shown in Figure 2. On average, the increase in velocity due to attenuation is ~0.3 km/s, but does not exceed ~0.5 km/s. The effect of attenuation on velocity is the smallest in the uppermost mantle lithosphere (<75-100 km depth) where attenuation estimates from Dalton et al. (2008) are the lowest. 3. Estimating cratonic geotherms The seismic velocity of a material cannot be directly translated into viscosity. Therefore, in this section we discuss the methods and assumptions needed to translate seismic velocity into temperature. In the next section, we calculate viscosities from the computed geotherms. 3.1 Velocity into temperature By correcting for attenuation in the last section, we are able to limit our discussion to the elastic shear velocity, which is dependent on pressure, temperature and composition. ! 98! Assuming a pyrolitic composition in thermodynamic equilibrium, and a typical upper mantle range of pressure and temperature, Stixrude and Lithgow-Bertelloni (2003) computed a model that relates isotropic shear wave velocities to temperature and pressure. This expression is written below, where Vs is shear velocity, P is pressure, z is depth, and T is temperature: (2) Within this formulation both the phase equilibria and the physical properties are considered jointly, with emphasis placed on thermodynamic self-consistency (Stixrude and Lithgow-Bertelloni, 2005). Seismic velocities for the individual phase components (of the pyrolite model) were computed along a 100 Ma geotherm, and Voigt-Reuss-Hill averaging was used to compute the shear velocity of the aggregate material. Estimated uncertainties for shear velocity using this expression average approximately 0.006 km/s (Stixrude and Lithgow-Bertelloni, 2005). Using the relationship between Vs, P, and T established in equation (2), we compute a geotherm for each of the velocity models in the Slave and Superior cratons (Fig. 3). Within the Slave craton, we find that temperatures increase steadily in models NA07 and ND08, before hitting maximums at a depth of 200-250 km, which suggests a transition from the base of a conductive thermal boundary layer to a convecting mantle adiabat. The profiles for ND08 and NA07 follow a similar trend in the Superior craton (Fig. 3), with a transition to adiabatic-like trends occurring at depths of 150-200 km. The ! 99! temperature profiles of Chen, SAWum_NA2 and SEMum2 are more complicated and require additional explanation. In the Chen profile (Slave craton, Fig. 3a), a jump in temperature at ~100 km depth is followed by a rapid decrease in temperature at ~150 km depth. This feature is due to the low velocity zone present in Figure 2. The magnitude of the feature in velocity is not well constrained (Chen et al., 2007) and therefore, we do not consider the amplitude of this thermal feature to be robust. Both the SAWum_NA2 and SEMum2 models display a decrease in temperatures at depths similar to the decrease in temperature observed in in the Chen model in both the Slave craton, and to a smaller degree, in the Superior craton (Fig. 3). These deviations in temperature/velocity fall well within the high velocity lithospheric lid, which implies that these thermal excursions exist in the thermally conductive lithospheric mantle. Such an excursion is difficult to explain in the lithosphere of these cratonic regions, where any thermal events capable of producing such gradients in temperature would have likely diffused away. Therefore, it seems more likely that the thermal anomaly that is being observed is due to incorrectly mapping a decrease in shear velocity into temperature. Possible reasons for a drop in velocity at shallow depths within cratonic lithosphere include a change in composition (Ford et al., 2010) or a boundary in anisotropy (Yuan et al., 2011). 3.2 Comparison to xenolith-constrained geotherms ! 100! To aid us in the interpretation of our computed geotherms, we superimpose geotherms calculated for the Slave and Superior cratons using a heat production model developed by Hasterok and Chapman (2011) (Fig. 3). The geotherms calculated by Hasterok and Chapman (2011) are derived from a partitioned heat production model in which the best- fitting model can be described by a ratio of basal heat flow and upper crustal radiogenic heat generation. The heat production models of Hasterok and Chapman (2011) are further constrained using thermal isostasy, which is corrected for compositional variations, and estimates of pressure and temperature from mantle xenoliths. The P-T estimates from mantle xenoliths for the Slave craton were obtained from approximately 110 garnet- peridotite xenoliths, taken from 6 studies. For the Superior craton, approximately 45 garnet-periodite xenoliths were obtained from two studies (Hasterok and Chapman, 2011). A significant feature of all velocity-derived models, and of the SAWum_NA2 and SEMum2 models in particular, is that their implied temperatures are significantly lower than the xenolith-constrained geotherms. Dalton and Faul (2010) reached a similar conclusion for shear velocities from the fastest cratonic regions in model S362ANI (Kustowski et al., 2008) relative to xenolith-constrained cratonic geotherms of (Pollack and Chapman; 1977; Artemieva, 2006). This discrepancy can be explained in three ways: 1) the velocity models overestimate absolute shear wave velocities, 2) the xenolith geotherms overestimate present-day temperature or 3) the velocities in the craton are higher than expected because of a composition different than what is assumed. For the time being, we consider only the ! 101! effects of composition on shear velocity and assume that the absolute velocities are correct and that the geotherms computed from the heat production model (Hasterok and Chapman, 2011) are also correct. Discussion regarding points 1) and 2) is saved for section 5.1. 3.3 Accounting for compositional variations The effect of varying composition on density and seismic velocities is important when trying to properly ascribe changes in velocity to changes in temperature and pressure. Using 133 samples of natural spinel- and garnet- peridotites, Lee (2003) formulated a relationship betwee Mg# and velocity by calculating the velocities of the natural peridotite samples, with a range of Mg# compositions, at standard temperature and pressure, using Hashin-Strickman averaging. The relationship between Mg# and velocity is given below: (3) The positive correlation between increasing velocity and increasing Mg# suggests that increased velocities in the mantle lithosphere beneath craton may be, in part, due to depletion. In North America, depleted cratonic mantle xenoliths have been found to have an Mg# (= Mg/(Mg+Fe) x 100) as high as 93 (Griffin et al. 2004). ! 102! To account for variations in composition we incorporate the Mg# correction factor derived by Lee (2003) into our velocity models. To assess the degree of increased depletion needed to reconcile the velocity-derived and xenolith-derived geotherms, we first calculate velocity residuals by translating the Hasterok and Chapman (2011) (HC) geotherms into velocity (HC-Vs), and then calculated the difference in velocity between HC-Vs and the attenuation-corrected, shear velocity models (Figs. 4a and 4b). Using the velocity residuals for each model, we can then determine the change in Mg# (dMg#) required to explain the difference between HC-Vs and the various shear velocity models (Figs. 4c and 4d). Theoretically, the range in dMg# can extend from 0 to 100. However, the pyrolitic model used in Stixrude Lithgow-Bertelloni (2003) has an Mg# of ~89, putting a cap of increasing dMg# at +11. Observations of depleted mantle lithosphere find that Mg# does not typically exceed 94, even in the most depleted cratonic lithosphere (Bernstein et al., 2007). This puts an upper bound on the Mg# and suggests that we cannot use a correction factor of greater than dMg# = +5. With such an upper bound, it is clear that melt depletion alone cannot reconcile the velocity models from tomography with the Hasterok and Chapman (2011) geotherms. Dalton and Faul (2010) reached a similar conclusion for cratonic regions with particularly high shear velocities. To provide reasonable upper bound for the influence of Mg# on velocity, we correct the shear velocity models with dMg# = +5 in the mantle lithosphere and a dMg# = 0 within the low velocity zone, assuming that the low velocity zone corresponds to the more fertile ! 103! mantle asthenosphere. Between the craton and asthenosphere depths, a gradient 75 km wide (in depth) was used to grade between dMg# = 5 and dMg# = 0. The location of this gradient was centered on the midpoint between the velocity maximum and minimum. As previously shown, if dMg# is capped at +5, velocities observed in the shear velocity models are still too high relative to the computed velocities of the geotherm HC-Vs. Another way of increasing velocity via composition is through the addition of garnet, whose seismic velocity is ~10% greater than the velocity of olivine (at a depth of 150 km) (Stixrude and Lithgow-Bertelloni, 2003). To include the addition of garnet into our velocity models, we derive a relationship between changing garnet content (dGt%) and velocity. dVs/dGt% is determined by the Hacker and Abers (2004) Excel macro, which calculates the densities and seismic velocities for a rock with mineral proportions described in terms of volume percent of end-member minerals. In our derivation of dVs/dGt%, we assume P-T conditions from the HC geotherm for the Slave craton. Shear velocities were computed for three different starting composition models, which were formulated using the mode percent and Mg# of the three garnet- peridotite samples listed in Lee (2003). Garnet was then added to each of the samples, with the volume percent added ranging between 0% and 30% (Table 1) (Fig. 5), and the corresponding velocities, calculated using Hashin-Shtrikman averaging (Hashin and Shtrikman, 1962), were recorded (Table 2). Differences between the Mg# of the starting and garnet-enhanced samples were corrected for (Lee, 2003), so that dVs/dGt% reflects the change in velocity of a material for which the Mg# of the bulk material is the same as ! 104! the Mg# of the garnet. There is also a weak dependence on the magnitude of dVs/dGt% with respect to depth (Fig. 5), with a correction of ~5.5 m/s/% at depths of 40 km up to ~7.5 m/s/% at 250 km, with 7.5 m/s/% being roughly half as large as the correction of Vs for dMg# = 1. Using our dVs/dGt%, which is depth dependent, we calculate the amount of additional garnet needed to explain the difference between the Mg#-corrected shear velocity profiles and HC-Vs (Figs. 4e and 4f). Within the Slave craton, the maximum amount of additional garnet needed to explain the difference between the velocity models and HC- Vs is ~ 50%. Considering the amount of garnet already present in the pyrolite model used to calculate temperature, the total amount of garnet present approaches 65%. At a minimum the additional amount needed to explain the largest velocity residual is only ~25% (total ~40%). In the Superior craton, the maximum amount needed is roughly 60% (total 75%), but the minimum amount needed to completely eliminate the residuals is a garnet composition close to 40% (total 55%). Although garnet-peridotites can contain significant amounts of garnet, a reasonable range of values for mantle lithosphere is 2-15% (Lee, personal communication). Since the amount of garnet already included in the pyrolite model ranges from 10-15% at lithospheric depths (Stixrude and Lithgow-Bertelloni, 2005),, we cap the contribution of additional garnet at 10%, which would put the total garnet content at a generous 25%,. Similar to the change in Mg#, we linearly decrease the amount of added garnet until we have reached the low-velocity zone. ! 105! For independent confirmation of the amount of garnet needed, we produce temperature profiles for each of the velocity models using the Hacker and Abers (2004) Excel macro and a composition of the most depleted garnet-peridotite sample listed in Table 1 (Table 1 of Lee, 2003). These new profiles are shown superimposed on plots of temperature calculated using the Stixtude and Lithgow-Bertelloni (2005) formulation (Fig. 3c,d). In both the Slave and Superior cratons, the greatest discrepancy exists between the calculated temperature profiles and the xenolith-constrained geotherm (Fig. 3c,d) at a depth of 150 km. The minimum amount of additional garnet needed to explain these residuals varies from 35% (for a total of ~37%) in the Slave craton to 48% (total of ~50%) in the Superior craton. These estimates agree well with the estimates established using the method of Stixrude and Lithgow-Bertelloni (2005). At 150 km depth the computed density (Hacker and Abers, 2004) of a garnet-peridotite containing 35-50% garnet is between 3460 and 3510 g/m3, which is significantly larger than the density of 3367 g/m3 estimated for the asthenospheric mantle (PREM, Dziewonski and Anderson, 1981). The cratonic lithosphere is often inferred to be neutrally buoyant, with increases in density due to colder temperatures being offset by a decrease in density due to a more depleted composition (e.g., Jordan, 1988; Lee, 2003). The effect of adding 35-50% garnet would make it significantly more difficult to maintain neutral buoyancy, and would likely result in the presence of large geoid anomalies in the vicinity of cratons, which are not observed (Mooney and Vidale, 2003). Assuming an Mg# of 100, the difference in densities would be reduced by approximately ! 106! 100 g/m3 (Lee, 2003), which is large enough to reconcile the smallest difference in calculated densities. However, such a depleted composition is not observed in cratonic lithosphere, which evidence suggests does not exceed an Mg# of 94 (Bernstein et al., 2007). Furthermore, the extreme melting temperature required to produce Mg# of 100 is inconsistent with thermal models, After accounting for the compositional effects on seismic velocity, we translate our new velocity models into temperature, using equation (2) for the Slave and Superior cratons (bold lines in Fig. 6). These new profiles are superimposed on the temperature profiles calculated prior to the compositional corrections. Within the Slave craton, all five profiles now fall within the HC geotherm at depths of less than 100 km. At greater depths the geotherms of NA07 and ND08 fall on the edge of the temperature range of the HC geotherm. The misfits between Chen et al. (2007), SAWum_NA2 and SEMum2 and the HC geotherm in the Slave craton are not eliminated, but are significantly reduced at depths of 100-200 km. Similarly, while all four models fall outside of the HC geotherm within the Superior craton, the discrepancies have been significantly reduced. 4. Viscosity of the cratonic lithosphere 4.1 Translating temperature into viscosity The rheological behavior of rocks is frequently described in terms of a power law behavior in which the strain rate (ė) is dependent on the differential stress (σ). In this paper, we assume a dislocation creep flow law, in the form written below: ! 107! (4) The effective viscosity (η) profile can then be calculated by assuming a strain rate (ė) and solving for stress (σ): (5) where A is a constant, n is the stress exponent, COH is the water concentration, r is the water concentration exponent, E* is activation energy, P is pressure, V* is activation volume, R is the gas constant, and T is temperature. The decision to use a dislocation creep flow law is due, in part, to the abundant evidence for the presence of seismic anisotropy at upper mantle depths, and because modeling results indicate that dislocation creep is the dominant deformation mechanism in the upper mantle (e.g. Behn et al., 2009). The values for the parameters A, E*, V* , r and n used in this paper are listed in Table 3 and were taken from Hirth and Kohlstedt (2003), where constraints on these parameters are discussed in detail. The water concentration assumed when producing the viscosity profiles shown in this paper is 50 H/106Si, which is reasonable considering the cratonic lithosphere is thought to be relatively dry (e.g., Hirth et al., 2000). ! 108! In order to accurately quantify the range of cratonic viscosities, we compute uncertainties in the Slave viscosity profiles for a range of activation energies (E*) and volumes (V*) (Fig. 7). The range of values considered for E* (=440 – 520 kJ/mol) and V* (=7 – 15 x10-6 m3/mol) are taken from the uncertainties discussed in Hirth and Kohlstedt (2003). The uncertainties estimated from E* are greatest at depths of ~125 km, before falling off at greater depths. For activation volume, V*, the uncertainties increase with depth. The peak uncertainty due to V* is less than an order of magnitude in viscosity, but is greater than the uncertainty present due to E*. Throughout the rest of this paper, these ranges of E* and V* are used to produce error bars on the computed viscosity profiles. 4.2 Viscosity profiles – Slave craton Viscosity profiles were calculated using two different, geologically reasonable, strain rates, ė = 10-15 1/s (a value determined to be suitable in the mantle beneath the Mojave Desert by Freed et al. (2012)) and ė = 10-14 1/s (Fig. 8). A smaller strain rate (ė = 10-15 1/s) results in larger viscosities than the larger strain rate (ė = 10-14 1/s), while the misfit between the various viscosity models remains the same, regardless of strain rate. Once the uncertainty in E* and V* is accounted for, the difference between the individual viscosity profiles is relatively small, and there is significant overlap between the models at most depths. The range of maximum lithospheric viscosities, inclusive of error, for a strain rate of 10-15 1/s is 8x1022 Pa*s to 4.2x1025 Pa*s, with an average maximum of 4.7*1024 Pa*s (Fig. 8). At a larger strain rate (ė=10-14 1/s) the average peak viscosity is ! 109! 6x1023 Pa*s. The Chen et al. (2007), SAWum_NA2 and SEMum2 models produce a peak in viscosity at ~125-150 km, which falls outside of the calculated viscosity range for the HC geotherm, while the NA07 and ND08 viscosity profiles follow the trend of the HC geotherm at lithospheric depths. The presence of a low-velocity layer within the uppermost mantle lithosphere of the Slave craton is observed at various depths in the models of Chen et al. (2007), SAWum_NA2 and SEMum2 (Fig. 2). This feature translates into a low-viscosity zone within the mantle lithosphere (Fig. 8). The correlation in depth between this low velocity zone and the location of a change in azimuthal anisotropy within the SAWum_NA2 model (Yuan et al., 2011) raises the possibility that velocity anisotropy may be partially responsible for the low velocity layer. The average minimum viscosity found in the asthenosphere, assuming a strain rate of ė=10-15 1/s, is 3x1021 Pa*s (with a total range of 3.2x1020 to 2x1022 Pa*s) beneath the Slave craton. At asthenospheric depths, the SEMum2 and SAWum_NA2 both fall within the anticipated range of viscosities assumed by the HC geotherm. For a larger strain rate of 10-14 1/s, the average minimum viscosity decreases to 5x1020 Pa*s beneath the Slave. 4.3 Viscosity profiles – Superior craton Similar to profiles of the Slave craton, viscosity profiles for the Superior craton were calculated using strain rates of 10-15 and 10-14 1/s. ! 110! The misfit between the viscosity profiles calculate from velocity models and the viscosity profile derived from the HC geotherm is in general larger than in the case of the Slave craton. At all but the shallowest depths (50-75 km) there is no overlap between the velocity-derived viscosity models and those from the HC geotherm (Fig. 8). There is however, significant overlap between the individual velocity-derived viscosity models at many depths, although the ND08 viscosity profile provides a noteable exception at the shallowest depths, with an estimate peak viscosity of 7x1033 Pa*s. This viscosity appears to be unreasonable and is not considered in the remaining discussion. On average, for a strain rate of 10-15 1/s, the peak in maximum lithospheric viscosities (excluding ND08) is 3x1024 Pa*s, and the minimum viscosity in the mantle asthenosphere is 5.6x1021 Pa*s. Increasing strain rate (10-14 Pa*s) reduces the average lithospheric viscosity to 5.6x1023 Pa*s, and 1x1021 Pa*s in the asthenosphere. 5. Implications 5.1 Uncertainties in the lithosphere Despite invoking compositional variations, there is still mismatch between the velocity- derived viscosities and the HC geotherm at lithospheric depths. This mismatch is most apparent at depths of 100-150 km in the Slave craton with the second group of viscosity models (those derived from SAWum_NA2, SEMum2, and Chen), and for all models at ! 111! depths greater than 75 km in the Superior craton. One possible explanation is that the attenuation correction is overestimated in the lithosphere, resulting in velocities that are too high. This explanation however, does not fully account for the viscosity misfit. Another possibility is that the original velocity models overestimate velocity. Within the Slave craton, it is certainly possible that one group of models, and therefore certain techniques, is better at determining velocities than the other. However, the systematic offset between the velocity models and the xenolith-constrained geotherm within the Superior craton suggests that additional confounding factors, either related to the velocity models themselves, or to other issues, exist. One such possibility is that the geotherms calculated by Hasterok and Chapman (2011) are too warm. Large uncertainties exist in our understanding of heat generation within the crust and mantle lithosphere, and it is possible that this uncertainty translates into errors within the HC model. A comparison to other steady-state geotherm families (Artemieva and Mooney, 2001; Chapman, 1986; Rudnick and Nyblade, 1999) suggests that temperatures could be overestimated by as much as 150°C at depths of 150 km in the HC model (Hasterok and Chapman, 2011) if a heat flow of 40 mW/m2 is assumed. Errors in xenolith P-T estimates can also impact the computed HC geotherms. Uncertainties due to assumptions regarding phase equilibrium could produce an error of at most 100°C (Saal, personal communication). While the errors associated with the construction of the HC geotherms are not small, they cannot fully explain the disparity in temperature, and ultimately viscosity. Reconciling the differences between velocity models and xenolith P-T estimates is not a trivial matter, and would help us to significantly improve our understanding the physical and ! 112! mechanical properties of the upper mantle. So while we have begun to broach the topic in this paper, further work is needed. 5.2 Asthenospheric properties The range of mantle viscosities inferred for the asthenosphere vary between models, but average ~1021 Pa*s (for ė=10-15 1/s) in both the Superior and Slave cratons. This value is in good agreement with studies of glacial isostatic adjustment which suggest that the viscosity of the asthenospheric mantle beneath cratons is on the order of 1021 Pa*s (Sella et al., 2007). The misfit between most of the velocity-derived viscosity models and the HC viscosity profile is significant, particularly in the Superior craton. If real, this discrepancy suggests that the asthenosphere beneath the North American craton is colder and stronger than implied by a 1300°C mantle adiabat. Asthenospheric viscosities in the western U.S., estimated to be 1019 Pa*s (Freed et al., 2012), are two orders of magnitude smaller than the averaged asthenospheric viscosities of our models. If taken at face value, this suggests that the asthenosphere is colder and more viscous beneath the North American craton than beneath young continental lithosphere. The processes that could have produced and maintained such lateral variations in asthenospheric viscosity are something of a puzzle since the average rate of North American plate motion would prevent the craton from remaining in place long enough to significantly cool the asthenosphere directly beneath it. Conversely, the low viscosities beneath the western U.S. could be lower than expected due to the presence of ! 113! melt or water. However, it is important to note that while the strain rate used to calculate the viscosity structure is constant, it is variable within the earth, due to changes in pressure. In the western U.S., the presence of asthenospheric material at shallower depths (lower pressure) will result in higher strain rates, and therefore smaller viscosity values. 5.3 A boundary in viscosity One of the most robust features of the viscosity models for the Slave and Superior cratons is the systematic difference in viscosities between the lithosphere and asthenosphere (as defined seismically). Within the Slave craton, the largest possible change is four orders of magnitude (excluding ND08; errors included) in SEMum2, and the smallest possible change is two orders of magnitude in NA07. While in the Superior craton, the smallest possible change is one order of magnitude, again in model NA07. A strong, pronounced variation in viscosities between the lithosphere and asthenosphere is in agreement with the idea that the lithosphere is a relatively rigid plate riding on top of a relatively low- viscosity, convecting mantle asthenosphere. The implication of these findings suggests that the rheologically defined lithosphere-asthenosphere boundary is being imaged seismically. 6. Conclusions We calculated viscosity profiles for the Slave and Superior cratons using isotropic shear wave velocity models from several global-, continental- and regional- scale models. The ! 114! resulting viscosity models show a clear contrast in viscosity between the lithosphere and asthenosphere, consistent with the rheological definition of the lithosphere-asthenosphere boundary. In the lithosphere, a maximum viscosity of 1024 Pa*s is observed. Asthenospheric viscosities on the order of 1021 Pa*s agree with independent estimates based on post-glacial rebound, but are two orders of magnitude larger than estimates of asthenospheric viscosity for the western U.S. We also find that significant misfit exists between the velocity-derived and xenolith- constrained temperature profiles. A portion of this misfit can be explained by invoking a depleted (Mg# 94) composition, consistent with xenolith data, and by increasing garnet content in the mantle lithosphere. ! 115! References Artemieva, I.M., Global 1× 1 thermal model TC1 for the continental lithosphere: implications for lithosphere secular evolution. Tectonophysics 416, 245-277 (2006). Artemieva, Irina M., and Walter D. Mooney. "Thermal thickness and evolution of Precambrian lithosphere- A global study." Journal of Geophysical Research106.B8 (2001): 16387-16414 Bedle, H., and S. van der Lee, S velocity variations beneath North America. Journal of Geophysical Research 114, B7308 (2009). Bernstein, S., P.B. Kelemen, and K. Hanghøj, Consistent olivine Mg# in cratonic mantle reflects Archean mantle melting to the exhaustion of orthopyroxene. Geology 35, 459-462 (2007). Cammarano, F., A. Deuss, S. Goes, and D. Giardini, One-dimensional physical reference models for the upper mantle and transition zone: combining seismic and mineral physics constraints. Journal of Geophysical Research, 110, B01306 (2005). Cammarano, F., and B. Romanowicz, Insights into the nature of the transition zone from physically constrained inversion of long-period seismic data. Proceedings of the National Academy of Sciences 104, 9139-9144 (2007). Chapman, D. S. "Thermal gradients in the continental crust." Geological Society, London, Special Publications 24.1 (1986): 63-70. Chen, C., et al., New constraints on the upper mantle structure of the Slave craton from Rayleigh wave inversion. Geophysical Research Letters, L10301 (2007). Dalton, C. A., G. Ekström, and A.M. Dziewonski, Global seismological shear velocity and attenuation: A comparison with experimental observations. Earth and Planetary Science Letters 284, 65-75 (2009). Dalton, C. A., G. Ekström, and A.M. Dziewonski, The global attenuation structure of the upper mantle. Journal of Geophysical Research 113, B09303 (2008). Dalton, C.A. and U.H. Faul, The oceanic and cratonic upper mantle: Clues from joint interpretation of global velocity and attenuation models. Lithos 120, 160-172 (2010). Dziewonski, A.M. and D.L. Anderson. Preliminary reference Earth model. Physics of the Earth and Planetary Interiors 25, 297-356 (1981). Evans, R.L., et al., Electrical lithosphere beneath the Kaapvaal craton, southern Africa. Journal of Geophysical Research 116, B04105 (2011). Fischer, K.M., et al., The lithosphere-asthenosphere boundary. Annual Review of Earth and Planetary Sciences 38, 551-575 (2010). Ford, H.A., et al., The lithosphere–asthenosphere boundary and cratonic lithospheric layering beneath Australia from Sp wave imaging. Earth and Planetary Science Letters 300, 299-310 (2010). Forsyth, D.W., and A. Li, Array-analysis of Two-dimensional Variations in Surface Wave Velocity and Azimuthal Anisotropy in the Presence of Multipathing Interference, in Seismic Earth: Array Analysis of Broadband Seismograms, (A. Levander and G. Nolet, ed.), AGU Geophysical Monograph 157, 81-97. French, S.W., V. Lekic, and B.A. Romanowicz, Toward global waveform tomography with the SEM: Improving upper-mantle images. AGU Fall Meeting Abstracts. Vol. 1. 2011. ! 116! Freed, A.M., G. Hirth, and M.D. Behn, Using short-term postseismic displacements to infer the ambient deformation conditions of the upper mantle. Journal of Geophysical Research 117.B01409 (2012). Gaherty, J.B., M. Kato, and T.H. Jordan, Seismological structure of the upper mantle: a regional comparison of seismic layering. Physics of the Earth and Planetary Interiors 110, 21-41 (1999). Griffin, W. L., et al., Lithosphere mapping beneath the North American plate. Lithos 77.1 (2004): 873-922. Hacker, B.R., and G.A. Abers, Subduction Factory 3: An Excel worksheet and macro for calculating the densities, seismic wave speeds, and H2O contents of minerals and rocks at pressure and temperature. Geochemistry Geophysics Geosystems 5, Q01005 (2004). Hashin, Z.S., and S. Shtrikman, On some variational principles in anisotropic and nonhomogeneous elasticity. Journal of the Mechanics and Physics of Solids 10, 335- 342 (1962). Hasterok, D., and D.S. Chapman, Heat production and geotherms for the continental lithosphere. Earth and Planetary Science Letters 307, 59-70 (2011). Hirth, G., and D.L. Kohlstedt, Water in the oceanic upper mantle: implications for rheology, melt extraction and the evolution of the lithosphere. Earth and Planetary Science Letters 144, 93-108 (1996). Hirth, G., and D.L. Kohlstedt, Rheology of the upper mantle and the mantle wedge: A view from the experimentalists. Geophysical Monograph- American Geophysical Union 138, 83-106 (2003). Hirth, G., R.L. Evans, and A.D. Chave, Comparison of continental and oceanic mantle electrical conductivity: Is the Archean lithosphere dry. Geochemistry Geophysics Geosystem 1, 1030 (2000). Jacobsen, S.D., et al., Effects of hydration on the elastic properties of olivine. Geophysical Research Letters 35, L14303 (2008). Jordan, Thomas H., Structure and formation of the continental tectosphere. Journal of Petrology 1 (1988): 11-37. Karato, S., and H. Jung, Water, partial melting and the origin of the seismic low velocity and high attenuation zone in the upper mantle. Earth and Planetary Science Letters 157, 193-207 (1998). Karato, S., and H.A. Spetzler, Defecr microsystems in minerals and solid-state mechanisms of seismic wave attenuation. Reviews of Geophysics 28, 399-421 (1990). Kennett, B.L.N., E.R. Engdahl, and R. Buland, Constraints on seismic velocities in the Earth from traveltimes. Geophysical Journal International 122, 108-124 (1995). Kustowski, B., G. Ekström, and A.M. Dziewoński, Anisotropic shear-wave velocity structure of the Earth's mantle: a global model. Journal of Geophysical Research 113, B06306 (2008). Lebedev, S., and R.D. Van Der Hilst, Global upper‐mantle tomography with the automated multimode inversion of surface and S‐wave forms. Geophysical Journal International 173, 505-518 (2008). Lee, C.-T.A., Compositional variation of density and seismic velocities in natural peridotites at STP conditions: implications for seismic imaging of compositional ! 117! heterogeneities in the upper mantle. Journal of Geophysical Research 108, 2441 (2003). Lee, C.-T.A., P. Luffi, P., E. Chin, Building and destroying continental mantle, Annual Reviews of Earth and Planetary Sciences 39,59-90 (2011). Lekić, V., and B. Romanowicz, Inferring upper‐mantle structure by full waveform tomography with the spectral element method. Geophysical Journal International 185, 799-831 (2011). Mooney, Walter D., and John E. Vidale. "Thermal and chemical variations in subcrustal cratonic lithosphere: evidence from crustal isostasy." Lithos 71.2 (2003): 185-193. Nettles, M., and A.M. Dziewoński, Radially anisotropic shear velocity structure of the upper mantle globally and beneath North America. Journal of Geophysical Research 113, B02303 (2008). Pollack, H.N., and D.S. Chapman, On the regional variation of heat flow, geotherms, and lithospheric thickness. Tectonophysics 38, 279-296 (1977). Resovsky, J., J. Trampert, and R.D. van der Hilst, Error bars for the global seismic Q profile, Earth and Planetary Science Letters 230, 413–423 (2005). Ritsema, J., et al., S40RTS: a degree‐40 shear‐velocity model for the mantle from new Rayleigh wave dispersion, teleseismic traveltime and normal‐mode splitting function measurements. Geophysical Journal International 184, 1223-1236 (2011). Romanowicz, B., The thickness of tectonic plates. Science 324, 474-476(2009). Romanowicz, B., A global tomographic model of shear attenuation in the upper mantle. Jounral of Geophysical Research 100, (1995). Rudnick, Roberta L., and Andrew A. Nyblade. "The thickness and heat production of Archean lithosphere: constraints from xenolith thermobarometry and surface heat flow." Mantle petrology: Field observations and high pressure experimentation: A tribute to Francis R.(Joe) Boyd 6 (1999): 3-12. Sella, G.F., et al., Observation of glacial isostatic adjustment in stable North America with GPS. Geophysical Research Letters 34, 2306 (2007). Stixrude, L., and C. Lithgow-Bertelloni, Mineralogy and elasticity of the oceanic upper mantle: Origin of the low-velocity zone. Journal of Geophysical Research 110, B03204 (2005). Van der Lee, S., and G. Nolet, Upper mantle S velocity structure of North America. Journal of Geophysical Research 102 (1997). Williams, M.L., K.M. Fischer, J.T. Freymueller, B. Tikoff, A.M. Trehu, and others, Unlocking the Secrets of the North American Continent: An EarthScope Science Plan for 2010-2020, National Science Foundation, February, 2010, 78 pp. Yang, Y., D.W. Forsyth, and D.S. Weeraratne, Seismic attenuation near the East Pacific Rise and the origin of the low-velocity zone, Earth and Planetary Science Letters 258, 260-268 (2007). Yuan, H., et al., 3‐D shear wave radially and azimuthally anisotropic velocity model of the North American upper mantle. Geophysical Journal International 184, 1237-1260 (2011). ! 118! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! Figure' 1.' ' (a)! Map! overview! of! major! tectonic! provinces! (after! Canil,! 2008)! in! Canadian!North!America.!!Red!and!blue!filled!circles!indicate!locations!sampled!by! shear! wave! velocity! models! NA07,! ND08,! SAWum_NA2,! and! SEMum2.! ! Averaged! profiles! were! calculated! for! each! model! within/near! the! Slave! craton! (red! circles)! and! the! Superior! craton! (blue! circles).! ! (b)! Example! of! shear! velocity! profiles! through!the!Slave!and!Superior!cratons!for!models!SAWum_NA2!and!SEMum.!!Thick! lines!in!each!panel!correspond!to!the!average!of!the!thinned!lines!within!the!same! panel.!!Thin!lines!correspond!to!the!individual!points!in!part!(a).' ! 119! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! Figure'2.''Absolute!shear!wave!velocity!profiles!for!the!Slave!and!Superior!cratons,! before!(thick!solid!lines)!and!after!(thin!solid!lines)!taking!into!account!attenuation! effects.!!(a)!Profiles!from!ND08,!NA07,!SAWum_NA2!and!SEMum2!are!averaged!from! points! shown! in! Figure! 1.! ! The! shear! velocity! profile! from! Chen! et! al.! (2007)! is! calculated! as! a! 1D! profile! for! the! entire! Slave! craton.! ! (b)! Same! as! (a),! but! for! the! Superior!craton,!without!the!Chen!et!al.!(2007)!model.!' ! 120! a. Slave Craton b. Superior Craton ! 50 50 1300 ºC 1300 ºC 100 100 H & H C & Depth (km) 150 150 C 200 200 Chen ND08 Nettles 250 NA07 250 NA07 SAWumNA2 SAWumNA2 SEMum2 SEMum 0 500 1000 1500 0 500 1000 1500 Temperature (°C) Temperature (°C) c. Slave Craton d. Superior Craton ! 50 50 1300 ºC 1300 ºC 100 H 100 & H C & C Depth (km) Depth (km) 150 35% 150 48% added Garnet added Garnet 200 200 Chen ND08 Nettles 250 NA07 250 NA07 SAWumNA2 SAWumNA2 SEMum2 SEMum 0 500 1000 1500 0 500 1000 1500 Temperature (C) Temperature (C) ! ! ! ! ! ! ! ! ! Figure'3.''(a)!&!(b)!temperature!profiles!calculated!from!the!attenuationScorrected,! shear! velocity! profiles! shown! in! Figure! 2,! using! the! relationship! between! velocity! (Vs)!and!temperature!(T)!derived!by!Stixrude!and!LithgowSBertelloni!(2005)!for!the! Slave!craton!(a)!and!the!Superior!craton!(b).!!(c)!and!(d)!show!coarsely!discretized! thermal!models!computed!using!a!depleted!garnetSperidotite!composition!(Table!1)! and!the!Hacker!and!Abers!(2004)!Excel!macro.!!The!arrows!in!(c)!and!(d)!indicate! the!amount!of!additional!garnet!needed,!on!top!of!the!original!~2.5%,!to!remove!the! temperature! residual! at! 150! km.! ! The! solid! black! line! labeled! “H&C”! in! (a)S(d)! corresponds! to! the! bestSfitting! geotherm! calculated! for! each! craton! (Hasterok! and! Chapman,!2011).!!The!dashed!lines!represent!the!95%!confidence!intervals!for!the! bestSfitting!geotherm.!!The!1300°C!adiabat!is!also!shown.' ! 121! ! ! ! ! ! ! ! ! ! ! ! ! ! ! Figure' 4.' ' (a)! and! (b).! ! Residuals! in! velocity! for! each! craton! and! shear! velocity! model.!!Residuals!were!calculated!by!translating!the!predicted!geotherm!(Hasterok! and!Chapman,!2011)!into!velocity!using!the!relationship!from!Stixrude!and!LithgowS Bertelloni!(2005)!and!comparing!to!the!respective!velocity!models!for!each!craton.!! (c)! and! (d)! The! increase! in! Mg#! (=100! x! Mg/(Mg! +! Fe))! required! to! eliminate! residual!velocities!(a!and!b),!determined!using!the!relationship!dVs/dMg#,!from!Lee! (2003).! ! (e)! and! (f)! The! increase! in! garnet! (volume! %)! needed! to! eliminate! the! velocity!residual,!determined!after!increasing!Mg#!by!5!at!lithospheric!depths.' ! 122! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! Figure'5.''(a)!Plot!of!shear!wave!velocity!as!a!function!of!added!garnet!(volume!%)! using!GarnetSPeridotite!1”!from!Lee!(2003)!and!listed!in!Table!1.!!Colors!indicate!for! velocities! calculated! at! different! points! in! pressureStemperature! (PST)! space.! ! PSTs! correspond! to! the! bestSfitting! geotherm! calculated! for! the! Slave! craton! from! Hasterok!and!Chapman!(2011).!!!(b)!The!change!in!velocity!as!a!function!of!changing! garnet! content! (dVs/dGt%)! plotted! as! a! function! of! increasing! depth.! ! The! dVs/dGt%! calculated! at! each! depth! is! an! average! of! the! dVs/dGt%! calculated! for! each!of!the!three!garnetSperidotite!samples!listed!in!Table!1!(from!Lee,!2003).' ! 123! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! Figure' 6.' ' Temperature! profiles! calculated! from! the! compositionally! corrected! shear!velocity!profiles!using!the!relationship!between!velocity!(Vs)!and!temperature! (T)!derived!by!Stixrude!and!LithgowSBertelloni!(2005)!for!the!Slave!craton!(a)!and! the!Superior!craton!(b).!!Thick,!solid!lines!are!for!the!composition!corrected!profiles,! while!the!thin!lines!are!profiles!corrected!for!attenuation!only,!and!are!the!same!as! those!in!Fig.!3.!!Geotherms!and!adiabats!are!the!same!as!in!Figure!3.' ! 124! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! Figure' 7.' ' Viscosity!profiles!for!the!Slave!craton,!calculated!using!a!wet!dislocation! creep!flow!law!(Hirth!and!Kohlstedt,!2003).!!(a)!Profiles!calculated!using!a!range!of! activation!energies!(E*).!!The!resulting!range!of!viscosities!for!each!profile!is!shown! with!light!shading.!!The!thick,!solid!lines!correspond!to!a!viscosity!calculated!using! E*! =! 480! kJ/mol.! ! (b)! Profiles! calculated! using! a! range! of! activation! volumes! (V*).!! Range!of!V*!is!shown!with!light!shading,!with!a!V*!=!11!m3/mol!used!to!calculate!the! thick,!solid!lines.' ! 125! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! Figure 8. Viscosity!profiles!calculated!for!the!Slave!(a,c)!and!Superior!(b,d)!cratons.!! (a)!and!(b)!assume!a!strain!rate!of!10S15!1/s,!while!(c)!and!(d)!assume!a!strain!rate!of!! 10S14!1/s.!!Viscosities!in!all!panels!are!calculated!using!a!wet!dislocation!creep!flow! law!(Hirth!and!Kohlstedt,!2003),!and!E*!=!440S520!kJ/mol!and!V*!=!7S15!m3/mol. ! 126! Garnet-Peridotite 1 % GT added 0% 1% 3% 5% 10% 20% 30.3% OPX 10.1 9.8 9.1 8.4 6.8 3.4 0.0 Mg2Si2O6 9.2 8.9 8.3 7.7 6.2 3.1 0.0 Fe2Si2O6 0.9 0.8 0.8 0.7 0.6 0.3 0.0 CPX 14.3 14.0 13.3 12.7 11.0 7.7 4.2 CaMgSi2O6 12.9 12.6 12.0 11.4 9.9 6.9 3.8 CaFeSi2O6 1.4 1.4 1.3 1.3 1.1 0.8 0.4 GT 14.1 15.1 17.1 19.1 24.1 34.1 44.4 Fe3Al2Si3O12 2.1 2.3 2.6 2.9 3.6 5.1 6.7 Ca3Al2Si3O12 1.9 2.1 2.4 2.6 3.3 4.7 6.1 Mg3Al2Si3O12 10.1 10.8 12.2 13.6 17.2 24.3 31.6 OL 61.4 61.1 60.4 59.8 58.1 54.8 51.3 Mg2SiO4 55.7 55.4 54.8 54.2 52.7 49.7 46.6 Fe2SiO4 5.7 5.7 5.6 5.6 5.4 5.1 4.8 Garnet-Peridotite 2 % GT added 0% 1% 1.1% 5% 10% 20% 30% OPX 20.5 20.1 20.1 18.0 15.5 10.5 5.5 Mg2Si2O6 19.0 18.7 18.7 16.7 14.4 9.7 5.1 Fe2Si2O6 1.4 1.4 1.4 1.3 1.1 0.7 0.4 CPX 0.4 0.0 0.0 CaMgSi2O6 0.3 0.0 0.0 CaFeSi2O6 0.0 0.0 0.0 GT 4.6 5.6 5.6 9.6 14.6 24.6 34.6 Fe3Al2Si3O12 0.7 0.8 0.8 1.4 2.2 3.7 5.2 Ca3Al2Si3O12 0.6 0.8 0.8 1.3 2.0 3.4 4.8 Mg3Al2Si3O12 3.2 4.0 4.0 6.8 10.4 17.5 24.6 OL 74.6 74.3 74.3 72.1 69.6 64.6 59.6 Mg2SiO4 68.7 68.4 68.4 66.4 64.1 59.5 54.9 Fe2SiO4 5.9 5.9 5.9 5.7 5.5 5.1 4.7 Table 1. Continued on next page ! 127! Garnet-Peridotite 3 % GT added 0% 1% 5.4% 8% 10% 20% 30% OPX 24.5 24.1 22.7 20.5 19.5 14.5 9.5 Mg2Si2O6 23.1 22.8 21.4 19.3 18.4 13.7 8.9 Fe2Si2O6 1.4 1.4 1.3 1.2 1.1 0.8 0.5 CPX 1.8 1.5 0.0 CaMgSi2O6 1.7 1.4 0.0 CaFeSi2O6 0.1 0.1 0.0 GT 2.4 3.4 7.8 10.4 12.4 22.4 32.4 Fe3Al2Si3O12 0.4 0.5 1.2 1.6 1.9 3.4 4.9 Ca3Al2Si3O12 0.3 0.5 1.1 1.4 1.7 3.1 4.5 Mg3Al2Si3O12 1.7 2.4 5.5 7.4 8.9 16.0 23.1 OL 71.3 71.0 69.5 67.3 66.3 61.3 56.3 Mg2SiO4 66.7 66.4 65.0 62.9 62.0 57.3 52.7 Fe2SiO4 4.6 4.6 4.5 4.4 4.3 4.0 3.6 Table 1. Table of volume percent of end-member compositions for three garnet- periodite samples, in which additional garnet is added incrementally. Original garnet- peridotite (0%) based on weight percent and Mg# of garnet-peridotites listed in Lee (2003). ! 128! Garnet-Peridotite 1: Vs (m/s) dVs/dG T P 0% 1% 3% 5% 10% 20% 30.3% % 425 1.1 4708.12 4713.92 4725.61 4737.32 4766.76 4826.15 4888.33 6.01 573 2 4686.83 4692.80 4704.81 4716.85 4747.10 4808.15 4872.04 6.10 729 3 4660.10 4666.61 4679.25 4690.04 4722.07 4784.91 4850.69 6.16 889 4 4630.45 4637.13 4650.19 4661.35 4694.40 4759.27 4827.19 6.55 1044 5 4601.70 4608.54 4621.99 4633.51 4667.52 4734.34 4804.29 6.55 1200 6 4571.86 4578.85 4592.71 4604.58 4639.57 4708.35 4780.34 6.95 1487 8 4520.58 4527.84 4542.39 4554.85 4591.52 4663.66 4739.17 7.22 Garnet-Peridotite 2: Vs (m/s) dVs/dG T P 0% 1% 1.1% 5% 10% 20% 30% % 425 1.1 4755.68 4761.54 4762.07 4785.29 4808.67 4855.85 4903.60 4.82 573 2 4724.51 4730.55 4731.09 4755.37 4780.08 4829.91 4880.35 5.21 729 3 4691.91 4698.19 4698.42 4723.43 4749.91 4803.43 4856.01 5.53 889 4 4656.14 4662.65 4662.87 4689.10 4717.05 4773.53 4829.12 5.74 1044 5 4620.60 4627.33 4627.54 4654.95 4684.34 4743.72 4802.25 5.92 1200 6 4584.76 4591.69 4591.90 4620.48 4651.29 4713.54 4774.99 6.26 1487 8 4521.58 4528.88 4529.07 4559.75 4593.10 4660.50 4727.19 6.91 Table 2. Continued on next page ! 129! Garnet-Peridotite 3: Vs (m/s) dVs/dG T P 0% 1% 5.4% 8% 10% 20% 30% % 425 1.1 4755.68 4761.73 4789.26 4827.52 4840.76 4886.95 4933.14 5.86 573 2 4724.51 4731.73 4759.26 4797.52 4810.76 4856.95 4913.14 6.28 729 3 4691.91 4701.73 4729.26 4767.52 4780.76 4836.95 4883.14 6.42 889 4 4656.14 4661.73 4689.26 4727.52 4740.76 4796.95 4863.14 6.82 1044 5 4620.60 4631.73 4659.26 4697.52 4710.76 4766.95 4833.14 6.97 1200 6 4584.76 4591.73 4619.26 4657.52 4670.76 4736.95 4803.14 7.38 1487 8 4521.58 4521.73 4559.26 4597.52 4610.76 4676.95 4753.14 7.79 Table 2. Shear!velocities!calculated!for!garnet8peridotite!samples!(listed!in!Table!1)!for!varying!temperatures!and!pressures.!! The!best!fitting!slope!for!each!point!in!temperature8pressure!space!is!listed!on!the!far!right!(dVs/dGt%).!!A!corresponding!plot! for! garnet8peridotite! 1! is! shown! in! Figure! 5.! ! The! average! of! the! three! sets! of! values! is! used! when! computing! the! garnet! correction!needed!(Figure!4).! ! 130! Dislocation Creep Parameters "best fit" range A 90 3.1-2528 n 3.5 COH 50 H/106Si r 1.2 E* 480 440-520 kJ/mol V* 11 4-15 m3/mol Table&3.&&Parameters!used!in!viscosity!calculations.!!The!“best!fit”!corresponds!to! the!thick!solid!lines!in!Figures!7!and!8,!while!the!“range”!is!used!to!compute!the! error!bars!in!Figure!8.& ! 131! Chapter 4 Localized shear in the deep lithosphere beneath the San Andreas Fault system By: Heather A. Ford1, Karen M. Fischer1 and Vedran Lekic2 1 Department of Geological Sciences, Brown University 2 Department of Geology, University of Maryland 132 Abstract The extent and geometry of strike-slip plate boundaries in the deep mantle lithosphere is an important, yet unresolved, aspect of plate tectonics. Models range from localized shear zones that are deep extensions of individual crustal faults to broad zones of diffuse, distributed shear with widths of hundreds of kilometers1,2,3,4. Here we use seismic waves that convert from shear to compressional motion (Sp) to image the base of the lithosphere across the San Andreas Fault system. Amplitudes of Sp conversions indicate a systematically smaller shear-wave velocity gradient with depth across the lithosphere- asthenosphere boundary on the western side of the plate boundary than to its east. The change in shear-wave velocity gradient typically occurs over a horizontal length scale of less than 50 km. This result is best explained by the juxtaposition of mantle lithospheres with different properties across the plate boundary. The spatial correlation between the surface expression of the plate boundary and the laterally abrupt change in LAB properties points to the accommodation of relative plate motion on a narrow shear zone that extends throughout the entire thickness of the lithosphere. 133 The San Andreas fault system, comprising the San Andreas fault and adjacent (sub)parallel strike-slip faults (Fig. 1), accommodates approximately 75% of the relative motion between the North American and Pacific plates5 and is dominated by 23-37 mm/yr slip rates on the San Andreas fault5,6. At mantle depths, however, the width of the plate boundary zone and the distribution of shear across it have remained unclear. Deformed xenoliths beneath the Calaveras fault1 and a step in crust-mantle boundary (Moho) depth across the southern San Andreas7 are consistent with shear due to plate motion at the top of the mantle lithosphere, and shear-wave splitting in SKS phases2,8,9 has been used to argue for a ~130 km wide shear zone in the lithosphere beneath central California1,2. However, these studies do not directly constrain the distribution of shear in the deep mantle lithosphere. In contrast, converted seismic waves can image the lithosphere-asthenosphere boundary (LAB), including possible variations in LAB structure due to shear zones that extend through the entire thickness of the lithosphere. Prior converted phase studies10,11,12,13,14 have estimated lithospheric thicknesses in California, but did not note systematic changes in LAB properties across the San Andreas fault system. To explore changes in LAB structure across the plate boundary we stacked over 135,000 Sp receiver functions from 730 seismic stations (Fig. 1) using a common conversion point (CCP) approach (Methods Summary). An Sp phase that represents a velocity decrease with depth is observed throughout the study region at depths of 45-100 km, with an average depth of 70 km (Figs. 2 and 3). Because the depth of this phase falls within the transition from high velocity lithosphere to underlying low velocity asthenosphere 134 imaged by surface wave tomography15,16, and because phase amplitudes correspond to a significant drop in velocity (Supplementary Information) we interpret the phase as a conversion across the seismologically-defined LAB. A striking feature of our CCP Sp model is the pronounced change in LAB phase amplitude across the plate boundary (Figs. 2 and 3a). To the east of the San Andreas fault (e.g. cross-section 1-1’ in Fig. 2) the LAB phase is strong and spatially coherent, while to the west of the San Andreas (2-2’ in Fig. 2) the LAB phase is either very low amplitude or absent. In central California, the change in LAB phase amplitude occurs beneath the San Andreas fault (A-A’ in Fig. 2 and Fig. 3a), and in northern California the shift in amplitude occurs further eastward (for example at the Calaveras fault, B-B’ in Fig. 2). In southern California, LAB phase amplitudes are lower on the western side of the San Andreas fault, with the exception of a high amplitude patch beneath the Inner Borderlands (Fig. 3a) where LAB depth is also anomalously shallow12 (Fig. 3b). The observed patterns in LAB phase amplitude cannot be explained by spatial variations in Sp path density or ray parameter (Supplementary Information). When our study region is divided into crustal blocks defined by faults17,18 (Fig. 3a), the average amplitude of the LAB Sp phase beneath blocks to the east of the San Andreas that move with velocities close to North American motion is 70% higher than phase amplitudes beneath blocks to the west of the San Andreas with velocities close to Pacific plate motion (0.0680 ± 0.016 versus 0.0403 ± 0.0128) (Fig. 4a). Beneath two blocks in 135 northern California that have intermediate relative motions (ECCR and WCCR) phase amplitudes are low and comparable to blocks on the Pacific side of the San Andreas. Unlike amplitudes, LAB phase depths do not vary systematically across the San Andreas fault system, and LAB phase depths and amplitudes are not correlated throughout the region (Fig. 3a). However, local variations in LAB phase depth and/or amplitude are observed in the vicinity of the Walker Lane fault system to the east of the Sierras, the Isabella Anomaly19 in the Great Valley, and elsewhere (Figs. 1 and 3) (Supplementary Information). LAB phase depths are less reliable along the western edge of the CCP stack where phase amplitudes are very weak, and near its eastern edge where energy in the LAB phase depth range is complex. To assess vertical gradients in absolute shear velocity (Vs) at the LAB, we modeled the ratio of LAB phase amplitudes for two blocks that are located on adjacent sides of the San Andreas fault in central California, the Salinas (SALI) block to the west and the Great Valley Thrust Belt (GVTB) block to the east (Supplementary Information). If the decrease in Vs from lithosphere to asthenosphere occurs instantaneously in depth beneath both sides of the San Andreas, a 1% drop at the base of the SALI lithosphere would imply a 3.7% drop at the GVTB LAB. For an LAB Vs gradient over 30 km in depth on both sides, a 1% drop beneath the SALI block would indicate a 7.9% drop beneath the GVTB block. Alternatively, the LAB Vs gradient beneath SALI could occur over a much larger depth range than the gradient beneath GVTB. This modeling confirms that LAB 136 velocity gradients differ significantly between the east and west sides of the San Andreas fault. The Sp CCP stack thus shows that the portions of Californian lithosphere moving with velocities close to Pacific plate motion or with intermediate plate velocities have an LAB that is characterized by a smaller and/or more gradual decrease in shear velocity, in comparison to lithosphere that translates with velocities close to that of the North American plate. In addition, this transition in LAB properties occurs over short horizontal distances, typically 50 km or less (Figs. 2 and 3). At the periods used in this study and for depths of less than ~70 km, CCP-stacked Sp phases can distinguish lateral variations in discontinuity structure that occur over 50 km versus larger distances, but they cannot resolve lateral variations of less than 50 km12. The laterally abrupt change in LAB amplitude indicates a contrast in the properties of the mantle lithosphere across the plate boundary. Sharp lateral gradients correlated with the plate boundary would be hard to maintain in the low viscosity20 asthenosphere, which is being sheared at an angle nearly normal to the plate boundary by the absolute motion of the North American plate21. Assuming an asthenosphere of roughly constant Vs across the region, shear velocity in the lithosphere moving with the Pacific plate would be smaller than on the North American side. While this difference is seen in some regional tomographic models22,23, it is not consistently present in all models15,16,24, suggesting that Vs may not be reduced throughout the entire mantle lithosphere; indeed, smaller Sp amplitudes from the LAB on the Pacific side could be matched simply by a reduced Vs 137 contrast or a more gradual gradient at the base of the plate (Fig. 4b), which could be undetected by tomography. Variations in Vs anisotropy at the LAB are not a likely explanation for the differences in LAB phase amplitude to the east and west of the plate boundary (Supplementary Information). What processes could have produced a reduced Vs contrast or a more gradual gradient at the base of the lithosphere on the Pacific side of the San Andreas fault in central and southern California and the Calaveras/Bartlett Springs faults (the eastern edge of the ECCR block) in northern California? One possibility is that the observed difference in LAB properties across the fault system results from strike-slip plate boundary processes (within the last 30 Ma). For example, the lithosphere on the Pacific side of the plate boundary has an absolute plate motion in multiple reference frames that is two to three times higher than the absolute plate velocity on the North American side21, possibly producing greater heat flow due to viscous heating in the shearing asthenosphere27. The trouble with this hypothesis is that viscosities of the asthenosphere are estimated to be relatively low in the region20, making it difficult to produce significant amounts of shear heating (Kincaid and Silver, 1996). However, if the temperature of the mantle is close to its solidus, then even a slightly higher temperature due to greater rates of shear in the Pacific-side asthenosphere could produce a small amount of melt. If the buoyant melt percolated upward, over time a layer of concentrated melt at the top of the asthenosphere would then act to reduce the width of the effective shear zone, increasing shear heating (Turcotte and Schubert). If Vs in the lithosphere is higher because the lithosphere is more dehydrated, compositionally depleted and colder than the asthenosphere28,29, then greater 138 heat flow entering the base of the Pacific-side lithosphere would reduce the Vs contrast by increasing temperature and counteracting the contrast in composition. Alternatively, the small fractions of melt produced by shear heating could permeate the base of the lithospheric plate, making the LAB Vs gradient more gradual in depth on the Pacific side of the plate boundary, thus explaining the lower Sp amplitudes. Another scenario is that in the subduction zone that existed prior to the onset of right-lateral strike-slip motion between the Pacific and North American plates at ~30 Ma, all of the lithosphere now to the west of San Andreas and Calaveras/Bartlett Springs faults acquired distinct LAB properties relative to the lithosphere now to the east of the faults. Although crust from similar subduction zone terranes is present on both sides of the San Andreas fault25,26, there is also evidence for the existence of subducted lithospheric material extending to the edge of the plate boundary beneath much of northern and central California31, often inferred to be stalled microplates32. These microplates are thought to have originated from the fragmentation of the Farallon plate, as portions of the ridge system at its western edge began to subduct32. In such a case, the stalled microplates would be both young and hot, and may not have had the time to fully produce a dry and depleted lithospheric layer33, resulting in a LAB Vs gradient that is more gradual. Regardless of the underlying process, the primary conclusions of this study are that mantle lithospheres with different basal properties are juxtaposed across the San Andreas fault in central and southern California and the Calaveras/Bartlett Springs faults in northern California and that the lateral transition between the different lithospheres typically occurs over 50 km or less at the LAB (Fig. 4b). These results imply that the 139 plate boundary expressed at the surface by the San Andreas fault system extends through the entire body of the lithosphere in a relatively narrow shear zone. 140 Methods Summary Sp phases are produced by the conversion of an incident S wave to a P wave across an interface in seismic velocity, and Sp receiver functions are obtained by deconvolving the parent (SV) phase from the daughter (P) phase. For this study, we obtained waveforms for events at epicentral distances of 55-85° and depths of less than 300 km from 19 broadband networks, including the Northern and Southern California Seismic Networks and the NSF EarthScope Transportable and Flexible Arrays. Waveforms were bandpass- filtered between 0.03 - 0.5 Hz, windowed around the S phase, rotated into P and SV components using a free surface transfer matrix optimized for free surface velocities at each individual station11, and deconvolved12. (See Supplementary Information for discussion of filtering and other methods.) The polarity of the Sp receiver functions was reversed to match the typical Ps convention in which a negative phase corresponds to a velocity drop with depth. In order to construct a 3D model of mantle structure, individual Sp receiver functions were migrated and stacked at their conversion point locations (common conversion point or CCP stacking). Receiver function time series were migrated into space (depth, latitude, and longitude) using ray tracing in a 1D model that represents crust30 and mantle16 velocity structure beneath the station at which the waveform was recorded. Uncertainties in lithosphere-asthenosphere boundary depth due to uncertainties in the migration model are likely to be less than 5 km12. The migrated Sp receiver functions were stacked in a 3D model that was discretized at 0.1° increments in latitude and longitude, and 0.5 km in depth. At each node sampled by more than 300 paths, a 141 weighted average of individual receiver functions defines the CCP stack, with weights given by cubic spline functions that approximate the Sp phase Fresnel zone12. Acknowledgements Thanks to Don Forsyth, Greg Hirth, and Marc Parmentier for helpful discussion. Waveforms were obtained from the IRIS Data Management Center and the Northern and Southern California Earthquake Data Centers. This work was supported by the National Science Foundation EarthScope Program (EAR-0641772) 142 References 1. Titus, S. J., Medaris, L. G., Wang, H. F., & Tikoff, B. Continuation of the San Andreas fault system into the upper mantle: evidence from spinel peridotite xenoliths in the Coyote Lake basalt, central California Tectonophysics 429, 1-20 (2007).# 2. Bonnin, M., Barruol, G., & Bokelmann, G. H. Upper mantle deformation beneath the North American–Pacific plate boundary in California from SKS splitting. J. Geophys. Res. 115, B04306 (2010).# 3. Platt, J. P., & Behr, W. M. Lithospheric shear zones as constant stress experiments. Geology, 39, 127-130 (2011).# 4. Vauchez, A., Tommasi, A., & Mainprice, D. Faults (shear zones) in the Earth's mantle. Tectonophysics 558-559, 1-27 (2012).# 5. Molnar, P., & Dayem, K. E. Major intracontinental strike-slip faults and contrasts in lithospheric strength. Geosphere, 6, 444-467 (2010).# 6. Rolandone, F., et al. Aseismic slip and fault-normal strain along the central creeping section of the San Andreas fault. Geophys. Res. Lett. 35, L14305 (2008).# 7. Zhu, L. Crustal structure across the San Andreas Fault, southern California from teleseismic converted waves, Earth Planet. Sci. Lett. 179, 183-190 (2000).# 8. Özaleybey, S., & Savage, M. K. Shear-wave Splitting beneath Western United States in Relation to Plate Tectonics, J. Geophys. Res. 100, 18135-18149 (1995)# 9. Hartog, R., & Schwartz, S. Y. Depth-dependent mantle anisotropy below the San Andreas fault system: apparent splitting parameters and waveforms. J. Geophys. Res. 106, 4155-4167 (2001).# 10. Li, X., Yuan, X., & Kind, R. The lithosphere–asthenosphere boundary beneath the western United States. Geophys. J. Int. 170, 700-710 (2007).# 11. Abt, D. L., et al. North American lithospheric discontinuity structure imaged by Ps and Sp receiver functions. J. Geophys. Res. 115, B09301 (2010).# 12. Lekic, V., French, S. W., & Fischer, K. M. Lithospheric thinning beneath rifted regions of southern California. Science 334, 783-787 (2011).# 13. Kumar, P., Kind, R., Yuan, X., & Mechie, J. USArray Receiver Function Images of the Lithosphere-Asthenosphere Boundary. Seismol. Res. Lett. 83, 486-491 (2012).# 14. Levander, A., & Miller, M. S. Evolutionary aspects of lithosphere discontinuity structure in the western US. Geochem. Geophys. Geosyst. 13, Q0AK07 (2012).# 15. Rau, C. J., & Forsyth, D. W. Melt in the mantle beneath the amagmatic zone, southern Nevada. Geology 39, 975-978 (2011).# 16. Obrebski, M., Allen, R. M., Pollitz, F., & Hung, S. H. Lithosphere–asthenosphere interaction beneath the western United States from the joint inversion of body‐wave traveltimes and surface‐wave phase velocities. Geophys. J. Int. 185, 1003-1021 (2011).# 17. McCaffrey, R. Block kinematics of the Pacific–North America plate boundary in the southwestern United States from inversion of GPS, seismological, and geologic data. J. Geophys. Res. 110, B07401 (2005).# 18. McCaffrey, R., et al. Fault locking, block rotation and crustal deformation in the Pacific Northwest. Geophys. J. Int. 169, 1315-1340 (2007). 143 19. Zandt, G., et al. Active foundering of a continental arc root beneath the southern Sierra Nevada in California. Nature 431, 41-46 (2004). 20. Freed, A. M., Hirth, G., & Behn, M. D. Using short-term postseismic displacements to infer the ambient deformation conditions of the upper mantle. J. Geophys. Res. 117, B01409 (2012).# 21. Schellart, W. P., Stegman, D. R., & Freeman, J. Global trench migration velocities and slab migration induced upper mantle volume fluxes: Constraints to find an Earth reference frame based on minimizing viscous dissipation. Earth Planet. Sci. Lett, 88, 118-144 (2008).# 22. Lin, F. C., Ritzwoller, M. H., Yang, Y., Moschetti, M. P., & Fouch, M. J. Complex and variable crustal and uppermost mantle seismic anisotropy in the western United States. Nature Geosci. 4, 55-61 (2010).# 23. Lin, F. C., Schmandt, B., & Tsai, V. C. Joint inversion of Rayleigh wave phase velocity and ellipticity using USArray: Constraining velocity and density structure in the upper crust. Geophys. Res. Lett. 39, L12303 (2012). 24. Buehler, J. S., & Shearer, P. M. Localized imaging of the uppermost mantle with USArray Pn data. J. Geophys. Res. 117, B09305 (2012). 25. Barbeau, D. L. Jr., et al. U-Pb detrital zircon geochronology of northern Salinian basement and cover rocks, GSA Bulletin 117, 466–481 (2005). 26. Jacobson, C. E., et al. Late Cretaceous–early Cenozoic tectonic evolution of the southern California margin inferred from provenance of trench and forearc sediments, GSA Bulletin 123, 485–506 (2011). 27. Turcotte, D. L., & Schubert, G. Geodynamics: Applications of Continuum Mechanics to Geological Problems (Cambridge University Press, New York, 1982). 28. Hirth, G., & Kohlstedt, D. L. Water in the oceanic upper mantle: implications for rheology, melt extraction and the evolution of the lithosphere. Earth Planet. Sci. Lett. 144, 93-108 (1996). 29. Rychert, C. A., Rondenay, S., & Fischer, K. M. P-to-S and S-to-P imaging of a sharp lithosphere-asthenosphere boundary beneath eastern North America. J. Geophys. Res. 112, B08314 (2007). 30. Lowry A. R., & Pérez-Gussinyé M. The role of crustal quartz in controlling Cordilleran deformation. Nature 471, 353–357 (2011). 31. Brocher, Thomas M., et al. Synthesis of crustal seismic structure and implications for the concept of a slab gap beneath coastal California. International geology review 41, 263-274 (1999). 32. Nicholson, Craig, et al. Microplate capture, rotation of the western Transverse Ranges, and initiation of the San Andreas transform as a low-angle fault system. Geology 12, 491-495 (1994). 33. Hirth, Greg, and David L. Kohlstedt. Water in the oceanic upper mantle: implications for rheology, melt extraction and the evolution of the lithosphere. Earth and Planetary Science Letters 144, 93-108 (1996). 144 1 2 39 SA F CF HF GV SN 37 B A SA F IA 35 WTR IB 33 Figure 1. Map of study region. Red triangles are locations of seismic stations used in the study. White dotted lines delineate cross sections shown in Figure 2. Thin black lines are faults, with the San Andreas fault (SAF) shown as a thick yellow line. Solid grey lines are state boundaries. Other features include: Hayward fault (HF), Calavaras fault (CF), Great Valley (GV), Sierra Nevada (SN), Isabella Anomaly (IA), Western Transverse Range (WTR) and the Inner Borderlands (IB). 145 1 A B A 2 SA 1 Depth (km) Depth (km) Distance Along Profile (km) Distance Along Profile (km) 2 B SA H C B A 2 1 RF amplitude Depth (km) Depth (km) Distance Along Profile (km) Distance Along Profile (km) Figure 2. Vertical cross-sections through the Sp CCP stack. Negative phases (blue) indicate a velocity decrease with depth. Depths inferred for the LAB phase are shown with white squares (see Supplemental Information). A positive amplitude phase (red) at depths of 20-40 km is interpreted to be the Moho, although conversions from intra-crustal discontinuities, such as the base of sedimentary basins, also occur. The dashed black line corresponds to Ps estimates of Moho depth30. Cross-sections intersections are shown with gray labeling; major faults are shown with thick black lines. Surface topography is exaggerated times ten. The amplitude of the inferred LAB phase is significantly weaker in 2-2’ than 1-1’, and west of plate boundary in cross-sections A-A’ and B-B’. 146 a. b. WeBR ECCR WCCR EBnR SNev SA PA L G I VT N A B IN Paci YO EMOJ MOJA SG VENT AB A W NZ SA SA BA A LT SMOJ C AT NA J A | Amplitude of LAB phase | LAB Depth (km) Figure 3. a Map of LAB phase amplitude strength. Grey lines outline crustal blocks defined by faults and plate motion17,18. Amplitudes are significantly lower beneath crustal blocks translating close to Pacific Plate motion, relative to blocks translating close to North American Plate motion (see Fig. 4). b Map view of LAB depth. Dots show volcanism from the past 3 Ma (www.navdat.org). Amplitude and depth values were smoothed with a Gaussian filter (Supplementary Information). 147 a. b. 0. 1 0.09 0.08 EMOJ EBnR 0.07 WeBR GVTB PANA MOJA SNev 0.06 INYO CATA 0.05 SGAB SAF GVT WBAJ HoF ECCR LIT 0.04 SALI VENT Paci HO ? SP WCCR HE ? 0.03 RE AS TH 0.02 EN OS PH 0.01 ER 0 10 20 30 40 50 60 E Figure 4. a Mean LAB phase amplitude in each crustal block17,18 as a function of block velocity relative to North America and b Schematic view of deformation in the lithosphere beneath the San Andreas fault system. a Red: blocks east of the San Andreas Fault with velocities close to North American Plate motion. WCCR and ECCR (black) are also east of the San Andreas but have intermediate plate motion. Blue: blocks west of the San Andreas, with light blue for CATA which includes anomalously thin lithosphere and a high amplitude phase in the Inner Borderlands that may be related to recent lithospheric thinning12. Blocks were included only if CCP stack nodes were > 35 and covered more than 50% of the block area. b A shear zone extends to the base of the lithosphere beneath the San Andreas fault and at the LAB it is less than ~50 km wide. Faults in foreground are the San Andreas fault (SAF), the Hosgri Fault (HoF) and the Great Valley Thrust (GVT). 148 Appendix A Supplementary materials for: Localized shear in the deep lithosphere beneath the San Andreas Fault system By: Heather A. Ford1, Karen M. Fischer1 and Vedran Lekic2 1 Department of Geological Sciences, Brown University 2 Department of Geology, University of Maryland 149 1. Additional information on receiver function and CCP stack methods 1.1. Optimal receiver function filtering Band-pass filtering of waveforms has a significant impact on the phases observed in individual receiver functions and in CCP stacks. Filters should include a broad enough spectrum of frequencies to ensure stable deconvolution and accurately image robust, small-scale variations in structure. To illustrate the consequences of eliminating too much short period energy from Sp waveforms, we calculated synthetic receiver functions using a propagator matrix method8, filtered them with a variety of bandpasses, and calculated receiver functions. The velocity model consists of a 30 km thick crust and an LAB of variable depth (50-100 km) with a constant velocity gradient (5.2% drop in shear velocity over 0 km). The filters used to calculate the synthetics all had a low frequency corner at 0.03 Hz (~33 s), with an upper corner ranging from 0.5 Hz to 0.125 Hz (2 s to 8 s). For the model with an LAB depth of 50 km, the only filters that produce an LAB phase within ±5 km of the correct depth are the two filters with upper-corners periods of 2.5 s and 2 s (Fig. 1). For the lower-pass filters (4 s to 8 s), the synthetic LAB phase has a peak between 10 and 35 km deeper than the 50 km depth in the model (Fig. 1). To a certain point, the offset between input depth and that predicted by the synthetics decreases for each of these low-pass filters as the input depth increases. However, the synthetics for each of these three filters eventually reach a cross-over point where the depth of the synthetic LAB phase begins to underestimate the input LAB depth. Overall, the depth of the LAB phase cannot reliably be estimated from the long-period filters, and the offsets between input and predicted LAB depths are typically related to interference between the Moho phase and the LAB phase. 150 The bandpass filter applied to the data (33 s to 2 s) in the CCP stacks presented in the main text is one that accurately retrieves LAB depth in the synthetic receiver function modeling. If LAB interface depths in the real earth are in fact at depths of less than ~90 km, as the CCP stacks obtained with the 33 s to 2 s filter suggest is true over most of the study region (Fig. 3 of the main text), the synthetic modeling indicates that applying a 33 s to 8 s filter should result in deeper apparent LAB phase depths. Overall, this effect is observed in the CCP stacks with the 33 s to 8 s filter. With the 33 s to 8 s filter, the width of the Moho and LAB phases are significantly broader (Fig. 2 versus Fig. 2 in the main text) and apparent LAB depth is typically deeper and less variable (Fig. 3 versus Fig. 3 in the main text), although some exceptions exist. This comparison supports the decision to interpret the CCP stacks obtained with the 33 s to 2 s filter. However, we note that the decrease in Sp phase amplitude across the plate boundary towards the Pacific is observed in models for both filters. 1.2 LAB phase selection When identifying an LAB phase at a given location within the CCP stacked receiver function images we make two assumptions. The first assumption is that the impedance contrast between the seismologically-defined lithosphere and asthenosphere is large and/or sharp enough in depth to produce an observable converted phase. The second is that the transition from lithosphere to asthenosphere occurs within a depth range (40-130 km) defined by surface wave tomography studies9, 10. 151 The depth and amplitude maps shown in the main body of text, as well as the LAB phase picks on the cross-sections in the main text, are the result of semi-automated user picks that were smoothed with a rotationally symmetric, lowpass Gaussian filter to eliminate high frequency variations. However, a depth and amplitude map was also constructed using a completely automated picking algorithm. In the user-guided approach, once a given peak in the CCP stack was chosen, a Gaussian fit was applied to guarantee peak amplitude selection. The largest negative peak in the 40 km to 130 km depth range was typically selected, but where two negative peaks of comparable amplitude exist at closely spaced depths, the mean depth of the two peaks was chosen. In the automated approach, the LAB phase depth and amplitude were calculated from the mean of the ten most negative points between 40 and 130 km. Other approaches to estimating phase amplitude (for example the amplitude at the single point closest to the mean depth of the ten most negative points, or the mean amplitude of all points within 5 km of that depth) yielded very similar amplitude values. To be considered for an LAB depth estimate, a given geographic location in the CCP stack had to have a peak negative amplitude greater than 0.01. For both semi-automated and automated picks, an LAB phase was selected at a given geographic location only if the model node for that location at a depth of 70 km had more than 300 receiver functions with non-zero weights in the spline function approximation of the Sp Fresnel zone (see Methods section in main text) and if amplitude information is present at 20 km. These cut-offs are a very conservative threshold, meant to ensure that LAB depths and amplitudes were measured only where the CCP model space is densely 152 sampled by Sp receiver functions. Other than the western and northeastern edges of the model and a few zones to the east, the entire study region is very well-sampled (Fig. 4). In particular, the pattern of Sp path density (Fig. 4) is not correlated with Sp phase amplitude (Fig. 3 of the main text), and the low Sp amplitudes on the Pacific side of the plate boundary cannot be attributed to sparse data. This lack of correlation is also illustrated by a plot of user-picked LAB phase amplitudes versus the number of receiver functions with non-zero weight for the same stack at 70 km depth (Fig. 4 inset). Overall, the depths of the LAB Sp phases obtained with the two picking approaches agree well. At 89.3% of model nodes, the depths agree to within 10 km, and 73.2% agree to within 5 km. Where depths are different by more than ~20 km (Fig. 5), the two methods are typically identifying different phases as the LAB arrival, and these points largely lie to the west of the plate boundary (beneath crustal blocks with plate velocities close to Pacific plate motion) where phase amplitudes are small overall, or close to the eastern edge of the model where multiple strong peaks that are closely spaced in depth occur. The user picks typically identify the largest negative amplitude within the 40 km to 130 km depth range as the LAB Sp phase, but where two negative peaks of comparable amplitude exist and the LAB Sp phase is assigned the mean depth of the two peaks, which explains some of the discrepancies between the methods. Other differences between the methods relate to a relatively small number of cases where the typical automated-pick criteria would create large oscillations in LAB depth over short horizontal length scales, and lateral continuity of the phase is favored in the user-picks. 153 Agreement in LAB Sp amplitudes from the two picking methods is very good (Fig. 5). At the largest amplitudes, user-picked values are slightly smaller than those picked automatically. This trend is largely due to zones near the eastern edge of the model where two negative peaks of comparable amplitude that are closely spaced in depth often exist in the LAB depth range. For these cases the automatic pick has a higher amplitude because it is the most negative of the two peaks, while the amplitude from the hand- picking represents their mean depth. 1.3 Influence of ray parameter Because the amplitudes of Sp converted phases depend on their incidence angles, we investigated whether variations in incident S phase ray parameters in our dataset across the study region could influence the observed variations in Sp phase amplitudes. Epicentral distances in our waveform data vary from 55 to 85 degrees, corresponding to ray parameters of 0.1206 and 0.0891 respectively, when the ak135 model11 is assumed. A plot of the mean ray parameter at each node in the model at 70 km depth (weighted by the spline function value for each phase) (map view in Fig. 6) shows that the ray parameter distribution is not in general well-correlated with the observed Sp amplitudes (Fig. 3 of the main text) although there is a slight tendency for the nodes with the smallest amplitudes to span a lower (but still broad) range of mean ray parameters. We therefore determined correction factors that scale Sp LAB phase receiver function amplitudes from different ray parameters to a constant amplitude, applied these to each observed receiver function, and calculated a corrected CCP stack. The result of applying this correction factor leads to only small differences in amplitude (less than 5%) between the corrected 154 stack and that shown in the main text (left-side of Fig. 7), and all of the significant amplitude features are retained (right-side of Fig. 7). 2. Modeling the relative velocity drop across the plate boundary The difference in Sp receiver function LAB phase amplitude between the Pacific and North American sides of the plate boundary can be attributed to changes in either the drop in velocity at the LAB or the depth range over which the velocity drop is distributed. To better understand the range of permissible LAB velocity gradients, we modeled the relative amplitudes of the Sp LAB phases between the SALI and GVTB blocks (Fig. 4 of main text) with synthetic Sp receiver functions computed for velocity drops of 1-10% over LAB gradient thicknesses (depth ranges) of 0-40 km. Velocity drops of greater than 10% were not considered as they exceed the upper limit of velocities inferred from surface wave tomography10. The amplitudes of the LAB phases computed from the synthetics (scaled so that the largest LAB phase has an amplitude of 1) are shown in a plot of amplitude versus gradient thickness for different values of velocity drop (Fig. 8). For a given gradient thickness and assuming a velocity drop of 1% for the SALI block, we calculated the velocity drop beneath the GVTB block that would be required to match the ratio of the GVTB and SALI LAB phase amplitudes. Similarly, for a given gradient thickness and assuming a velocity drop of 10% for GVTB, we calculated the velocity drop beneath SALI that would be required to match their LAB phase amplitude ratio. These calculations define the range of differences in LAB velocity drop between the two blocks (ΔVsGVTB-ΔVsSALI) as a function of gradient thickness (bold bars in Fig. 8). We also repeated this process for the minimum and maximum ratios of the SALI and GVTB 155 phase amplitudes to define uncertainties (thin bars in Fig. 8). Although there is significant overlap in the error bars, the calculated drop in velocity needed to explain the difference in amplitudes increases with increasing LAB gradient thickness. Assuming the LAB thickness is the same on both sides, the best-fitting required differences in velocity (ΔVsGVTB-ΔVsSALI) range from 2.7% to 8.3%. 3. Other variations in LAB phase depth and amplitude Although they are not the focus of this study, a number of other interesting variations in LAB phase depth and amplitude are observed. Throughout most of the three-dimensional Sp CCP stack, the depth of the LAB phase varies gradually, with a few exceptions where the depth of the LAB phase changes rapidly over short (25-50 km) lateral distances (Fig. 3 of main text). For example, the Inner Borderlands are underlain by an anomalously shallow LAB (~50 km) (Fig. 9), thought to be the result of lithospheric thinning associated with the clockwise rotation of the Western Transverse Range block starting at ~19 Ma7. The westward decrease in LAB depth from Los Angeles and other onshore areas to the east into the Inner Borderlands is laterally abrupt (Fig. 3 of the main text), as is the north to south transition from the thicker lithosphere beneath the Western Transverse Range (~90 km) to the north. 156 A reduction in LAB depth to ~60 km is also found beneath much of Walker Lane, the zone of deformation to the east of the Sierra Nevada that accommodates approximately 25% of the right-lateral motion between the North American and Pacific plates12, 13, and shallow LAB depths are particularly pronounced beneath the northern Walker Lane. Levander and Miller14 also found a zone with shallow LAB depths between the Sierra Nevada and the Basin and Range in their combined analysis of Ps and Sp receiver functions. Because strike-slip deformation is oriented parallel to the NW-SE trend of Walker Lane and has a relatively modest magnitude, for example ~30 km of cumulative displacement in the northern Walker Lane13, this deformation alone does not provide a clear explanation for the relatively thin lithosphere beneath Walker Lane. Basins bounded by normal faults are found within northern Walker Lane, but they appear to be accommodating the same overall right-lateral strike-slip motion found elsewhere in the Walker Lane zone12, 13. Rather, the thin lithosphere beneath Walker Lane may predate the onset of strike-slip motion, and in fact could represent a pre-existing zone of lithospheric weakness that allowed plate boundary deformation to migrate into this region. The Isabella anomaly is a high velocity mantle feature seen in tomography models beneath the southern Great Valley. It was originally hypothesized to represent delamination of the Sierra Nevada lithosphere15, 16, 17 . More recently, some have suggested that it is evidence for a remnant fossil slab18. In the vicinity of the anomaly we observe an LAB phase that is shallower and higher amplitude than in the surrounding area. In Fig. 10 we show cross sections that intersect in the region containing the 157 anomalously shallow LAB phase. To ensure that these depth and amplitude anomalies are not artifacts of the thick sedimentary layer present in the Great Valley, we calculated synthetic Sp receiver functions for a velocity model containing an eight-layer crust, simplified from a seismic refraction experiment in the Great Valley19, 20, a lithospheric layer, and an asthenospheric halfspace. This modeling shows that an Sp phase from an LAB at 70 km underlying a complex crust containing a ~3-5 km thick sedimentary layer should still appear at 70 km in an Sp receiver function. Moreover, since a perturbed LAB Sp phase is not consistently observed in the CCP stack elsewhere beneath the Great Valley sediments, the anomalously shallow and high amplitude feature in the vicinity of the Isabella Anomaly appears to be structurally significant. The implications of this correlation warrant further study. Another area with an unusually shallow LAB exists between the Clear Lake Volcanic Field and Sutter Butte, two regions that have been volcanically active in the past 3 My (Fig. 10; Fig. 3 of main text). Regionally, the location of recent volcanism often (but not always) appears to be co-located with rapid lateral gradients in LAB depth (Fig. 3 of main text). A regional correlation between LAB depth variations, volcanism and seismicity was previously noted by Levander and Miller14. Another intriguing observation is the presence of a double positive phase in the northwestern portion of our CCP model, in the vicinity of the subducting Juan de Fuca slab. The first positive phase arrives at a depth of ~110-120 km and the second phase is observed at ~150 km (Fig. 11). Although there is a slight variation depth, these positive 158 phases do not obviously follow the dip of the subducting lithosphere. Further to the east, a zone with a particularly strong LAB phase can be found at the top of the mantle wedge above the subducting Juan de Fuca slab (Fig. 2 of the main text). This zone is located west of the magmatic arc, which itself is underlain by an average LAB phase of unremarkable amplitude. These results are consistent with models in which the upper plate lithosphere has been permeated with melt or otherwise effectively removed in the vicinity of the magmatic arc, while to the west melt is trapped at the base of the upper plate lithosphere, producing a particularly strong LAB Sp arrival. Other large offsets in LAB depth occur near the western and eastern edges of the Sp CCP stack (Fig. 3 of main text). However, low amplitude arrivals on the Pacific side of the plate boundary make identifying a clear Sp phase from the LAB difficult, and the existence of multiple arrivals that are closely spaced in depth beneath parts of the eastern Sierra Nevada and the Basin and Range also complicates selection of a single LAB phase. Thus, LAB phase depths in these two zones are somewhat suspect and should not be over-interpreted. 159 References 1. Kennett, B. L. N. The removal free surface interactions from three‐component seismograms, Geophys. J. Int. 104, 153–154 (1991). 2. Abt, D. L., et al. North American lithospheric discontinuity structure imaged by Ps and Sp receiver functions. J. Geophys. Res. 115, B09301 (2010). 3. Helffrich, G. Extended-time multitaper frequency domain cross-correlation receiver-function estimation. Bull. Seismol. Soc. Am. 96, 344-347 (2006). 4. Lowry, A. R., & Pérez-Gussinyé, M. The role of crustal quartz in controlling Cordilleran deformation. Nature 471, 353-357 (2011). 5. Obrebski, M., Allen, R. M., Xue, M., & Hung, S. H. Slab-plume interaction beneath the Pacific Northwest. Geophys. Res. Lett. 37, L14305 (2010). 6. Obrebski, M., Allen, R. M., Pollitz, F., & Hung, S. H. Lithosphere–asthenosphere interaction beneath the western United States from the joint inversion of body‐wave traveltimes and surface‐wave phase velocities. Geophys. J. Int. 185, 1003-1021 (2011). 7. Lekic, V., French, S. W., & Fischer, K. M. Lithospheric thinning beneath rifted regions of southern California. Science 334, 783-787 (2011). 8. Keith, C. M., & Crampin, S. Seismic body waves in anisotropic media: synthetic seismograms. Geophys. J. Roy. Astron. Soc. 49, 225-243 (1977). 9. Rau, C. J., & Forsyth, D. W., Melt in the mantle beneath the amagmatic zone, southern Nevada. Geology 39, 975-978 (2011). 10. Obrebski, M., Allen, R. M., Pollitz, F., & Hung, S. H. Lithosphere–asthenosphere interaction beneath the western United States from the joint inversion of body‐wave traveltimes and surface‐wave phase velocities. Geophys. J. Int. 185, 1003-1021 (2011). 11. Kennett, B. L. N., Engdahl, E. R., and Buland, R. Constraints on seismic velocities in the Earth from traveltimes. Geophys. J. Int., 122, 108-124 (1995). 12. Wesnousky, S. G. The San Andreas and Walker Lane fault systems, western North America: Transpression, transtension, cumulative slip and the structural evolution of a major transform plate boundary. J. Struc. Geol. 27, 1505-1512 (2005). 13. Wesnousky, S. G., Bormann, J. M., Kreemer, C., Hammond, W. C., & Brune, J. N. Neotectonics, geodesy, and seismic hazard in the Northern Walker Lane of Western North America: Thirty kilometers of crustal shear and no strike- slip?. Earth Planet. Sci. Lett. 329, 133-140 (2012). 14. Levander, A., & Miller, M. S. Evolutionary aspects of lithosphere discontinuity structure in the western US. Geochem. Geophys. Geosyst. 13, Q0AK07 (2012). 15. Zandt, G., et al. Active foundering of a continental arc root beneath the southern Sierra Nevada in California. Nature 431, 41-46 (2004). 16. Boyd, O. S., Jones, C. H., & Sheehan, A. F. (2004). Foundering lithosphere imaged beneath the southern Sierra Nevada, California, USA. Science 305, 660- 662. 17. Schmandt, B., & Humphreys, E. Complex subduction and small-scale convection revealed by body-wave tomography of the western United States upper mantle. Earth Planet. Sci. Lett. 297, 435-445 (2010). 160 18. Pikser, J. E., Forsyth, D. W., & Hirth, G. Along-strike translation of a fossil slab. Earth Planet. Sci. Lett. 331, 315-321 (2012). 19. Holbrook, W. S., & Mooney, W. D. The crustal structure of the axis of the Great Valley, California, from seismic refraction measurements. Tectonophysics 140, 49-63 (1987). 20. Colburn, R. H., & Mooney, W. D. Two-dimensional velocity structure along the synclinal axis of the Great Valley, California. Bull. Seismol. Soc. Am. 76, 1305- 1322 (1986). 161 1 0 0 60 80 100 160 180 10 5 0 50 55 60 65 70 75 80 85 90 95 100 Figure 1. (upper) Plot of Sp receiver functions calculated with different bandpass filters. All receiver functions normalized to maximum amplitude. (lower) Differences in depth between LAB discontinuities in the input velocity model and the predicted LAB phase in the synthetic Sp receiver functions. 162 1 1’ 0 40 Depth (km) 80 120 160 0 100 200 300 400 500 600 700 800 900 1000 Distance Along Profile (km) 2 2’ 0 40 Depth (km) 80 120 160 0 100 200 300 400 500 600 700 800 900 Distance Along Profile (km) 0 RF amplitude Figure 2 Figure 2. Cross-sections 1-1’ and 2-2’ striking parallel to the plate boundary. Cross sections are the same as those in Fig. 2 in the main body of text, but were calculated using a filter of 0.03-0.175 Hz. Compared to results obtained from bandpass filters with a higher frequency upper corner (Fig. 2) the cross sections contain phases that are significantly longer period and in general deeper. A lower LAB phase amplitude is still observed to the west of the plate boundary. 163 a. b. | LAB Amplitude | LAB Depth (km) Figure 3 Figure 3. Automatically picked LAB Sp phase amplitude (A) and depth (B) maps for CCP stacked Sp receiver function model calculated using a filter of 0.03-0.175 Hz (8-33 sec). (A) Similar to the 0.03-0.5 Hz filter map shown in Fig. 4 of the main text, lower amplitudes are observed west of the plate boundary. (B) With the 0.03-0.175 Hz filter, LAB phase depths are deeper and less variable than those with the 0.03-0.5 Hz filter shown in the main text. 164 39 37 35 33 Figure 4 Figure 4. Number of receiver functions contributing information to the spline function at each node point in the model at 70 km depth. White line outlines the spatial extent of model points containing more than 300 receiver functions, which is used as a lower limit for LAB phase selection. Inset is a plot of the number of receiver functions with non- zero weights in the spline function versus the amplitude of the user-picked LAB phase. Inset illustrates the lack of correlation between amplitude and number of receiver functions at a given model point. Red line lies at the 300 receiver function cut-off. 165 a. b. 130 0.2 120 0.18 0.16 110 0.14 100 0.12 90 0.1 80 0.08 70 0.06 60 0.04 50 0.02 40 0 40 50 60 70 80 90 100 110 120 130 0 0.05 0.1 0.15 0.2 Depth of automatically picked LAB phase (km) Amplitude of automatically picked LAB phase (km) Figure 5. Relationship between LAB depth (A) and LAB amplitude (B) determined using the automatic picking and hand picking methods describe in the Supplemental Information. (A) 89% of the depth picks for both methods fall within 10 km of each other, and 73% fall within 5 km. (B) User-picks versus automated-picks illustrate a tendency for user-picks to have slightly lower amplitudes. 166 Figure 6. Weighted mean ray parameter at 70 km for each point in the CCP stack. A subtle shift from smaller to larger mean ray parameter, which would result in increased LAB amplitude, occurs across the study region. Inset demonstrates the weak correlation between ray parameter and amplitude. 167 a. b. Percent Difference in Amplitude Amplitude of LAB phase Figure 7 Figure 7. (A) Percent difference in LAB Sp phase amplitude between the original CCP stack and the ray parameter corrected stack. (B) LAB Sp phase amplitude in the ray parameter corrected model. 168 1 0.9 Normalized Amplitude Relaitve to Maximum 0.8 0.7 Decreasing LAB 0.6 Velocity Drop- 10% to 1% 0.5 0.4 0.3 0.2 0 5 10 15 20 25 30 35 40 Gradient Thickness 9 Difference in LAB Velocity Drop (%) Across the SAF 8 7 6 5 4 3 2 1 0 0 5 10 15 20 25 30 35 40 Gradient Thickness Figure 8 Figure 8. (upper) LAB phase amplitude (with the maximum value normalized to 1) as a function of LAB gradient thickness for varying velocity contrasts, calculated from synthetic Sp receiver functions. Each line corresponds to a different value of LAB velocity drop. (lower) Difference in percent LAB velocity drop (ΔVsGVTB-ΔVsSALI) required to match the observed ratio in Sp phase amplitudes between the GVTB and SALI crustal blocks for different gradient thicknesses. Thick bars indicate best fitting range of values, and thin bars indicate maximum range of values (one standard deviation). 169 C’ C C C’ 0 40 Depth (km) 80 120 160 0 100 200 300 400 500 600 Distance Along Profile (km) 0 RF amplitude Figure 9 Figure 9. Cross-section C-C’ through the Sp CCP stack (0.03-0.5 Hz filter) that runs perpendicular to the plate boundary and intersects the Inner Borderlands. At the eastern edge of the Inner Borderlands (~150 km distance) an offset in LAB depth occurs. Observation matches results of Lekic et al. (2011) and is interpreted to be the result of lithospheric thinning due to the rotation of the Western Transverse Range lithosphere. 170 V X’ X X’ V 0 Depth (km) 40 80 X 120 V’ 160 0 100 200 300 400 500 Distance Along Profile (km) V V’ X 0 40 Depth (km) 80 120 160 0 100 200 300 400 500 600 700 800 900 1000 Distance Along Profile (km) 0 RF amplitude Figure 10 Figure 10. Cross-sections X-X’ and V-V’ through the Sp CCP stack (0.03-0.5 Hz filter). The shallow LAB imaged at the intersection of the cross-sections (“V” in X-X’ and “X” in V-V’) is located in the vicinity of the Isabella Anomaly imaged in surface and body wave tomography. 171 W Y’ Y Y Y’ W 0 Depth (km) 40 80 120 W’ 160 0 100 200 Distance Along Profile (km) W W’ Y 0 40 Depth (km) 80 120 160 0 100 200 300 400 500 600 700 800 900 1000 Distance Along Profile (km) 0 RF amplitude Figure 11 Figure 11. Cross-sections Y-Y’ and W-W’ through the Sp CCP stack (0.03-0.5 Hz filter). Two positive phases at different depths are imaged in the vicinity of the subducting Juan de Fuca slab (~0-200 km distance on W-W’ and ~74-200 km distance on Y-Y’). 172