Seismic Wavefield Imaging of the Earth: the Regional, the Local, and the Remote Citation Castellanos, Jorge Alberto Castillo (2022) Seismic Wavefield Imaging of the Earth: the Regional, the Local, and the Remote. Dissertation (Ph.D.), California Institute of Technology. doi:10.7907/fm32-1t74. https://resolver.caltech.edu/CaltechTHESIS:07192021-182024547 Abstract In this thesis, I use seismic wavefield methods to illuminate the interior structure and the dynamics of the Earth across different scales. First, I image the large-scale lithospheric structure at the eastern sector of the Trans-Mexican Volcanic Belt to constrain on the transition from flat to steeper subduction in central Mexico. Then, I move to a regional scale and image the dynamics of the Wallowa Mountain block in northeastern Oregon, where mantle-based stresses appear to have played an essential role in shaping the crustal structure. With the findings of this investigation, I was able to illuminate a deformation mechanism of mantle origin, which I also use here to explain other near-surface processes in different parts of the North America continent. After, I move to a local scale, where I use dense oil-industry instrumentation to image the sub-kilometer crustal structure of Long Beach, California. In the first part of this investigation, I use noise-derived surface waves to create a high-resolution shear wave velocity model of the first kilometer of the crust, which I use to numerically determine the variability in the expected ground shaking intensity of the area. In the second part, I move past the traditional surface wave analysis and use the body wave portion of the noise-derived Green's functions to create a high-resolution compressional wave velocity model beneath one of the surveys. Finally, I present a waveform-based method of analysis that shows great promise as a new way of investigating the seismic behavior and the physical conditions of isolated marine environments. Item Type: Thesis (Dissertation (Ph.D.)) Subject Keywords: Seismology, Ambient noise, Seismic Imaging Degree Grantor: California Institute of Technology Division: Geological and Planetary Sciences Major Option: Geophysics Thesis Availability: Public (worldwide access) Research Advisor(s): Clayton, Robert W. Thesis Committee: Jackson, Jennifer M. (chair) Zhan, Zhongwen Gurnis, Michael C. Ross, Zachary E. Clayton, Robert W. Defense Date: 8 July 2021 Non-Caltech Author Email: jcastillo (AT) google.com Record Number: CaltechTHESIS:07192021-182024547 Persistent URL: https://resolver.caltech.edu/CaltechTHESIS:07192021-182024547 DOI: 10.7907/fm32-1t74 Related URLs: URL URL Type Description https://doi.org/10.1029/2018JB015783 DOI Article adapted for Chapter 2. https://doi.org/10.1126/sciadv.abb0476 DOI Article adapted for Chapter 3. https://doi.org/10.1029/2019JB018855 DOI Article adapted for Chapter 6. https://doi.org/10.1029/2019JB018941 DOI Article adapted for Chapter 7. ORCID: Author ORCID Castellanos, Jorge Alberto Castillo 0000-0002-0103-6430 Default Usage Policy: No commercial reproduction, distribution, display or performance rights in this work are provided. ID Code: 14308 Collection: CaltechTHESIS Deposited By: Jorge Castillo Castellanos Deposited On: 30 Jul 2021 19:02 Last Modified: 04 Aug 2022 21:08 Thesis Files PDF - Final Version See Usage Policy. 145MB Repository Staff Only: item control page

SEISMIC WAVEFIELD IMAGING OF THE EARTH: THE REGIONAL, THE LOCAL, AND THE REMOTE Thesis by Jorge Alberto Castillo Castellanos In Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy in Geophysics CALIFORNIA INSTITUTE OF TECHNOLOGY Pasadena, California 2022 Defended July 8th , 2021 ii © 2022 Jorge Alberto Castillo Castellanos ORCID: 0000-0002-0103-6430 All rights reserved iii To my beloved grandmother, Irma iv ACKNOWLEDGEMENTS First and foremost, I would like to express my deepest appreciation to my advisors, Robert Clayton and Zhongwen Zhan. It is because of their incommensurable guidance and support that I am able to write this thesis and continue my academic journey with more enthusiasm than ever. Yet, above all else, I am grateful for their friendship, empathy and patience, all of which have shaped my experience as a graduate student and propelled my professional growth forward. I am extremely grateful to the members of my committee: Jennifer Jackson, Michael Gurnis, and Zachary Ross. I thank them for taking the time to know me both as a researcher and as a person, and for providing me with direction when I found myself floundering. I am also very much indebted to several Seismo Lab professors: Victor Tsai, Mark Simons, Joann Stock, Jean-Philippe, Pablo Ampuero, and Don Helmberger. It was a true privilege to sit through their lectures and learn from some of the most inspirational thinkers I have ever known. I would of course not have been able to get to this point of my PhD without the company and support of my friends and colleges at the Seismolab. To the Seismo Fam, thank you for all the good AGU memories and impromptu lunches. To the current and former members of the SMS-362 office, Semechah Lui, Yiran Ma, Kangchen Bai, Voonhui-Lai, Zhichao Shen, Yuan-Kai Liu and, for a brief yet enjoyable moment, Jack Muir, thank you for all the stimulating conversations and for being a constant source of encouragement. Last, but certainly not least, to the always-loving SM-252 office, Kim, Donna, Priscilla, Sarah, and Rosemary, thank you for making me feel like I never left home. I must also thank my former advisor and dear friend Xyoli Pérez-Campos for introducing to the world of seismology and for putting me in this path. Still to this day, knocking at her door has been the greatest decision I have ever made. This section would not be complete without acknowledging the people who put up with me outside the office. To the group of misfits at Caltech, Marcel, Amanda, Abbas, Bin, Akshita, Magnus, and Gautham, thank you for sharing your PhD with me and reminding me really matters when work occupied v too much of our lives. To the Mexican crew, Porfi, Ellen, Alejandro, Orika, Manuel, Daniela, Enrique, Manolo, and David, thank you for all those Ernie’s lunches and late-night mezcales–they will always be part of my most treasured memories at Caltech. To my fiancée Helen Han, I cannot express how fortunate I was for your arrival at this time of my life. Your support, warmth, and encouragement have been the most wonderful part of this journey and the biggest drive for me to move forward. Finally, I would like to give my most heartfelt thanks to my family. To Yanina and Francisco, thank you for teaching me the value of hard work and for always encouraging me to pursue my passions. To Javier and José, thank you for being my rock and for keeping my feet on the ground. To Irma, Claudia, Daniel, Daniela, and Fernanda, thank you for believing in me and being present in every step of the way. Everything I am is because of every one of you. vi ABSTRACT In this thesis, I use seismic wavefield methods to illuminate the interior structure and the dynamics of the Earth across different scales. First, I image the large-scale lithospheric structure at the eastern sector of the Trans-Mexican Volcanic Belt to constrain on the transition from flat to steeper subduction in central Mexico. Then, I move to a regional scale and image the dynamics of the Wallowa Mountain block in northeastern Oregon, where mantle-based stresses appear to have played an essential role in shaping the crustal structure. With the findings of this investigation, I was able to illuminate a deformation mechanism of mantle origin, which I also use here to explain other near-surface processes in different parts of the North America continent. After, I move to a local scale, where I use dense oil-industry instrumentation to image the sub-kilometer crustal structure of Long Beach, California. In the first part of this investigation, I use noise-derived surface waves to create a high-resolution shear wave velocity model of the first kilometer of the crust, which I use to numerically determine the variability in the expected ground shaking intensity of the area. In the second part, I move past the traditional surface wave analysis and use the body wave portion of the noise-derived Green’s functions to create a high-resolution compressional wave velocity model beneath one of the surveys. Finally, I present a waveform-based method of analysis that shows great promise as a new way of investigating the seismic behavior and the physical conditions of isolated marine environments. vii PUBLISHED CONTENT AND CONTRIBUTIONS J. C. Castellanos, R. W. Clayton, and X. Pérez-Campos (2018). “Imaging the eastern Trans-Mexican Volcanic Belt with ambient seismic noise: Evidence for a slab tear.” In: Journal of Geophysical Research 123.9, pp. 7741-7759. doi: 10.1029/2018JB015783. (Top 10% most downloaded JGR papers of 2018). J.C.C. processed the seismic data, performed the analysis, and wrote the manuscript. R.W.C. and X.P.-C. conceptualized the project. J. C. Castellanos, R. W. Clayton, and A. Juarez (2020). “Using a timebased subarray method to extract and invert noise-derived body waves at Long Beach, California.” In: Journal of Geophysical Research 125.5. doi: 10.1029/2019JB018855. (Paper selected for Eos Editor’s highlight). J.C.C. processed the seismic data, performed the analysis, and wrote the manuscript. A.J. help in designing the methodology. R.W.C. conceptualized the project. J. C. Castellanos, Z. Zhan, and W. Wu (2020). “Absolute centroid location of submarine earthquakes from 3D waveform modeling of water reverberations.” In: Journal of Geophysical Research 125.5. doi: 10.1029/2019JB018941. J.C.C. processed the seismic data, performed the analysis, and wrote the manuscript. W.W. helped in setting up the wavefield simulations. Z.Z. conceptualized the project. J. C. Castellanos, J. Perry-Houts, Robert W. Clayton, Y. Kim, A. C. Stanciu, B. Niday, and E. Humphreys (2020). “Seismic anisotropy reveals crustal flow driven by mantle vertical loading in the Pacific NW.” In: Science Advances 6.28. doi: 10.1126/sciadv.abb0476. J.C.C. processed the seismic data, wrote the manuscript. J.P.-H. performed the geodynamic modeling. Y.K., A.C.S., and B.N. helped in the interpretation of the data. R.W.C. and E.H. conceptualized the project. J. C. Castellanos, and R. W. Clayton (In review). “The fine-scale structure of Long Beach, California, and its impact on ground motion acceleration.” In: Journal of Geophysical Research. J.C.C. processed the seismic data, performed the analysis, and wrote the manuscript. R.W.C. conceptualized the project. J. C. Castellanos, E. Humphreys, and R. W. Clayton (In review). “Evidence for mantle-based deformation across the western US.” In: Geophysical Research Letters. J.C.C. processed the seismic data, performed the analysis, and co-wrote the viii manuscript. E.H. conceptualized the project and co-wrote the manuscript. R.W.C. conceptualized the project and co-wrote the manuscript. ix TABLE OF CONTENTS Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Published Content and Contributions . . . . . . . . . . . . . . . . . . . Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . List of Illustrations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter I: Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter II: The Enigmatic Trans-Mexican Volcanic Belt . . . . . . . . 2.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Data and Method . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . 2.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter III: The Dynamic Interaction Between the Crust and the Mantle: A Case Study in the Pacific NW. . . . . . . . . . . . . . . . . . 3.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Discussion and Conclusion . . . . . . . . . . . . . . . . . . . . Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter IV: The Dynamic Interaction Between the Crust and the Mantle: In Search for the Upside-Down Water Bed Model. . . . . . . . . 4.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Data and Methods . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Discussion and Conclusions . . . . . . . . . . . . . . . . . . . . Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter V: The Fine-Scale Structure at Long Beach, California: Surface Waves. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Data and Methods . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . 5.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter VI: The Fine-Scale Structure at Long Beach, California: Body Waves. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv vi vii ix xi 1 4 4 5 8 17 29 30 39 39 39 43 47 50 54 54 54 56 59 63 64 68 68 68 70 81 90 92 99 99 x 6.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Ambient Noise Correlation and Body Waves . . . . . . . . . . 6.4 Double Beamforming and Body Wave Extraction . . . . . . . . 6.5 Traveltime Measurements and Body Wave Tomography . . . . 6.6 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . 6.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter VII: A Novel Framework for Investigating Mid-Oceanic Environments. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.3 Waveform Modeling of Water-Reverberations . . . . . . . . . . 7.4 Earthquake Relocation by Waveform Modeling of Water Reverberations: A Case Study in the Gofar Transform Fault . . . . . 7.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter VIII: Conclusions and Future Direction . . . . . . . . . . . . . 99 102 104 110 117 120 121 126 126 127 129 132 142 148 148 155 xi LIST OF ILLUSTRATIONS Number Page 2.1 Tectonic setting of central and southern Mexico. The brown box indicates the main focus area of this study. The black contour lines depict the depth of the subducted slab compiled using the results of receiver functions, tomography studies, and hypocentral re-localization (from Ferrari et al., 2012). The dashed part of these contour lines indicates the transition between the flat and normal dipping portion of the South Cocos plate. The inverted white triangles show the location of the broadband seismic stations that have operated within this region. Regional seismicity at different depths, as reported by the Servicio Sismológico Nacional (SSN) since 1998, is shown as color-coded circles. The volcanic provinces: the Trans-Mexican Volcanic Belt (TMVB), the Anegada High (AH), the Los Tuxtlas Volcanic Field (LTVF), and the Modern Chiapanecan Volcanic Arc (MCVA) are delimited by the orange areas, and the location of the main stratovolcanoes and calderas are marked by the red triangles and red circles, respectively. The light purple region delimits the Veracruz basin (VB), and the light magenta region delimits the Jalisco Block. The green dashed line denotes the Rivera-Cocos plate boundary, and the purple dashed line denotes the projected path of the Orozco Fracture Zone (OFZ) beneath North America (Dougherty et al., 2012). The blue dashed line indicates the location where the slab exhibits a significant change in dip and where Dougherty and Clayton (2014) propose the existence of a trench-perpendicular tear. The ages of the plates (Ma) and the convergence rates (cm/yr) along the Middle America Trench (MAT) are shown in blue and red numbers, respectively. The inset map shows this work study area and the distribution of every station used in this study. . . . . . . . . . . . . . . . . . 6 xii 2.2 Elevation map of the eastern section of the Trans-Mexican Volcanic Belt showing the abrupt termination of the volcanic arc. Note how altitude drops more than 5,000 m from Pico de Orizaba to the coastal Veracruz basin over a horizontal distance of 120 km. This region corresponds to the one encompassed by the brown box in Figure 2.1. . . . . . . . . . . . . . . . . . . . . . 2.3 ZZ, RR, ZR, and TT cross-correlations as a function of interstation distance for the period band of 5-100 s. The traces are stacked by bins of 2 km and 0.5 s. As expected, the extracted Green’s functions show clear Rayleigh waves in the ZZ, RR, and ZR components and clear Love waves in the TT component; although several body wave arrivals emerge as well. Note the presence of the anomalous cS-cP core phase in the ZZ and ZR cross-correlations (Pham et al., 2018). . . . . . . . . . . . . . . 2.4 Ambient noise phase and group velocity maps (left), checkerboard resolution maps (center), and model error maps (right) at 34 s period. The volcanic provinces and the Veracruz basin are delimited by the solid and dashed contour lines, respectively. Note how the velocity distribution varies for all four velocity types due to their difference in sensitivity. . . . . . . . . . . . . 2.5 Cross-sections of the Rayleigh wave group (A) and phase (B) velocity models and the Love wave group (C) and phase (D) velocity models along the MASE line (F-F’ in Figure 2.1). Topography is shown above each profile. . . . . . . . . . . . . . . 8 10 13 14 xiii 2.6 Vertical cross-sections of the Voigt-averaged VS (A and B) and radial anisotropy models (C and D) along the MASE and VEOX lines (F-F’ and G-G’ in Figure 2.1, respectively). Topography is shown above each profile. The hypocenters from the SSN catalog within ±10 km are projected to each cross-section. The slab models from Pérez-Campos et al. (2008) and Melgar and PérezCampos (2011) are plotted in black lines assuming a constant plate thickness of 40 km. The thin black dashed line in profiles A, B, and D represent the Moho interface, defined by where the shear wave velocity first exceeds 4.1 km s−1 , and the thick green and red dashed lines in profiles C and D represent regions of elevated positive (≥VSH ) and negative radial anisotropy (>VSV ), respectively. Cross-sections E and F show the average misfit of the 1-D inversion along the same profiles. . . . . . . . . . . . . 2.7 Sensitivity kernels for the fundamental mode Rayleigh-wave phase velocity calculated from the modified Tectonic North America Model (A; Stubailo et al., 2012). Two-dimensional histograms over the backazimuth-phase velocity space for the 3-20 seconds period band (B) and the 20-50 seconds period band (C). The orange bars on top of each panel represent the number of crosscorrelations beamformed in each backazimuth bin. The dashed red line indicates the best fit of Equation 2.3 and the gray diamonds mark the beam velocity in each bin (with size as a function of beam power). The best fitting parameters are given at the bottom right of each panel. . . . . . . . . . . . . . . . . . . 2.8 Vertical cross-sections of the Voigt-averaged VS model along the A-A’ and H-H’ lines in Figure 2.1. Both profiles show extreme crustal thinning in the Veracruz basin. The dashed red line in panel (A) indicates the location where Dougherty and Clayton (2014) propose the existence of a plate tear. Topography is shown above each profile. Other symbols are as in Figure 2.6ab. 15 18 19 xiv 2.9 Vertical cross-sections of the Voigt-averaged VS model in the Veracruz basin (from the A-A’ to the E-E’ lines in Figure 2.1). The color scale is saturated to only show the velocity anomalies present in the upper mantle. The solid black line indicates the Moho interface. The thick red line indicates the location where Dougherty and Clayton (2014) propose the existence of a plate tear, and the dashed red line delimits the 75-km-wide zone where they observe a sharp change in intraslab seismicity. Cross-sections are separated by a constant offset of 20 km. The eye marker in Figure 2.1 indicates the viewpoint of the observer. 2.10 Vertical cross-sections of the radial anisotropy model in the Veracruz basin (same profiles as in Figure 2.9). . . . . . . . . . . 2.11 Map of the azimuthal anisotropy results for the lower crust and upper mantle. The orientation of the vectors gives the seismically fast direction, φ, and the length of the lines is proportional to the amplitude of the anisotropy, A. The thick green arrows depict the possible trajectories of mantle flows inferred from our measurements. The red oval delimits the eastern TMVB volcanic chain. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 22 24 xv 2.12 3-D schematic illustration of the proposed tectonic setting and inferred mantle flow (green arrows). Asthenospheric mantle materials are flowing almost trench-parallel beneath the slab in the forearc region until they are redirected to the mantle wedge through a tear that separates Central and South Cocos. The accelerated rollback rate of South Cocos relative to Central Cocos introduces a suction force that further displaces the mantle materials laterally toward the south of Mexico. The shear stress exerted by the toroidal flow around the slab tear may be transporting the melts and fluids that are feeding the volcanoes in the LTVF from central Mexico, thus allowing their existence without any slab material at an appropriate depth directly underneath them. Some of the active volcanoes composing the eastern TMVB volcanic chain however might still be fed by some slabedge melting mechanism. Anisotropy in the back-arc of central Mexico (i.e. Central Cocos) is primarily controlled by a trenchperpendicular 2-D corner flow that is induced by the down-dip motion of the slab. The flow of mantle materials through a slab tear, with the accompanying slab-edge melting, may explain the exceptional spatial distribution of stratovolcanoes in the eastern sector of the TMVB, as well as the abrupt change in the source of melts observed in young rocks in this segment of the MAT. The eye in Figure 2.1 marks the viewpoint of the observer. . . . 3.1 Regional and location maps for northeastern (NE) Oregon. (A) Regional map showing the broadband seismic stations used in this study (black inverted triangles). The dashed blue line depicts the Snake River Plain (SRP). WA, Washington; OR, Oregon; ID, Idaho. (B) Global map centered in NE Oregon. The thick black line encloses the region shown in (A). (C) Elevation map of the topographic bullseye region (red area in A). The dashed larger ellipse is the outer limit of the bullseye, whereas the inner ellipse locates the Wallowa batholith. . . . . . . . . . 27 41 xvi 3.2 Example of beamformer outputs and final azimuthal anisotropy model. (A) Two-dimensional histograms over the azimuth-velocity space for the 3- to 17-s period band with the best-fitting anisotropy model (red dashed lines). The green bars on top of each panel indicate the number of noise cross-correlations available for each azimuth. (B) Azimuthal anisotropy model for the crust of the Pacific NW. Bar orientation gives the fast direction of azimuthal anisotropy, and bar length is proportional to anisotropy amplitude. The background color represents the intersection density of the anisotropy vectors assuming that they are of infinite length (i.e., projected to the bounds of the study region). The green and blue dots indicate the location of the two stations beamformed in (A). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.3 Comparison of seismic and geodynamic results with crustal anisotropy. (A) Azimuthal anisotropy model for the crust of the Pacific NW overlying a depth slice through the Vp tomography model at 250 km (Schmandt and Humphreys, 2010). The red dashed lines depict the Wallowa anomaly and the Siletzia slab curtain. (B) Modeled Moho stress and mid-lower crustal flow velocity for the Pacific NW. The colored contours represent the vertical stress at the Moho based on a global geodynamic model driven by density anomalies derived from the P-wave velocity structure. The black arrows denote the predicted mid-lower crustal flow velocity that results from the application of the modeled Moho stress to a viscously heterogeneous crust. The red bars represent the anisotropy measurements derived from this study. . . . . . . . . 44 xvii 3.4 Station averaged shear wave splitting measurements for the Pacific NW. (A) SKS splitting measurements for the entire western U.S. (Becker et al., 2012). The red arrow depicts the relative motion between the North American plate (NA) and the hotspot reference frame (HS) (Gripp and Gordon, 2002). (B) SKS splitting measurements for our study region (red area in A). The thick blue vectors depict the measurements of Niday and Humphreys (2020), and the black vectors are from the database of Becker et al. (2012). The orientation of the vectors gives the angle of the fast polarization, and the length of the bars is proportional to the magnitude of the shear wave splitting. The white trajectories through the anisotropy field lines in (B) are used to represent the streamlines of the mantle flow assuming an east-oriented flow (Silver and Holt, 2002). The red region in the background marks the location of the Wallowa anomaly. Note how mantle materials appear to flow smoothly around the lateral boundaries of the Wallowa anomaly. . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.5 Schematic representation of the upside-down water-bed model. The load of the mantle lithosphere is a force creating vertical stresses (σzz ) on the Moho. The lithospheric load pulls down on the crust, which creates a lateral pressure gradient that drives Poiseuille flow in the ductile mid-lower crust. The asthenosphere flows independently (as evidenced by its independent anisotropy field; Figure 3.4), creating a local Couette flow that is decoupled from the mid-lower crust by the mantle lithosphere. . . . . . . 48 3.6 Crustal anisotropy and upper mantle velocity structure of California. The black vectors depict the surface wave anisotropy measurements at 12 s period (Lin et al., 2011). Bar orientation gives the fast direction of azimuthal anisotropy, and bar length is proportional to anisotropy amplitude. The background color corresponds to a depth slice through the Vp tomography model at 195 km (Schmandt and Humphreys, 2010). The red dashed lines denote the two seismically fast and likely dense mantle anomalies. 49 xviii 4.1 Upper mantle velocity structure of the western US and geography of the target regions. The main map shows a depth slice through the Vp tomography model at 195 km (Schmandt and Humphreys, 2010). The dashed black rectangle delimits the region of study of Castellanos et al. (2020), whereas the continuous black rectangles delimit the regions that are focused in this investigation. The zoomed maps show the distribution of broadband seismic stations used to resolve the crustal anisotropy beneath each region (red inverted triangles). . . . . . . . . . . . . . . . 4.2 Crustal properties of the western US crust. (A-B) Show the mid-crust and lower crust shear wave velocity structure as derived from ambient seismic noise and earthquake tomography, receiver functions, and Rayleigh wave ellipticity (H/V) measurements (Shen and Ritzwoller, 2016). (C) Shows the crustal thickness as derived from teleseismic P-to-S receiver functions, and Rayleigh wave phase velocities and ellipticity from noise interferometry and earthquakes (Schmandt et al., 2015). (D) Shows the surface topography (Ryan et al., 2009). (E-D) Shows the shortest Rayleigh wave period that is sensitive to the middle crust and the longest Rayleigh period that is sensitive to the lower crust, respectively. The average upper and lower period for the two target regions is shown in the bottom of each of the geographic borders. . . . . . . . . . . . . . . . . . . . . . . . . 57 58 xix 4.3 Lower crust and mantle anisotropy for the two target regions. (A) Lower crust azimuthal anisotropy around the Rocky Mountains. (B) Station-averaged shear-wave splitting measurements around the Rocky Mountains (Becker et al., 2012). (C) Lower crust azimuthal anisotropy around California. (D) Station-averaged shear-wave splitting measurements around California (Becker et al., 2012). The bar orientation gives the fast direction of anisotropy, and the bar length is proportional to anisotropy amplitude. The background color represents the mantle P-wave structure at 195 km depth (Schmandt and Humphreys, 2010). In panels (A) and (B) the green arrows depict our preferred interpretation of crustal flow. The blue arrows in (C) depict the relative motions between the Pacific (PA) and North American plates (NA), the Juan de Fuca plate (JdF) and NA. The purple lines in (C) mark the plate boundaries. ISA stands for the Isabella mantle anomaly. . . . . 60 xx 4.4 Proposed mechanisms that give rise to the observed crustal anisotropy. (A) Negatively buoyant mantle anomalies. The load of the mantle lithosphere creates vertical stresses on the Moho and pull the crust down, creating a lateral pressure gradient that drives Poiseuille flow in the ductile mid-lower crust towards the mantle anomaly. The green arrow shows the sign of the lithospheric load. (B) Positively buoyant mantle anomalies. The hot and buoyant mantle rises and pushes the crust upward, creating a lateral pressure gradient that drives Poiseuille flow in the ductile mid-lower crust outwards from the mantle anomaly. The green arrow shows the sign of the lithospheric load. (C) Plate interactions. For the case of a transform fault, the deformation in the crust and mantle lithosphere is dominated by strain in a narrow shear zone, whereas in the asthenosphere deformation is in a zone whose width increases with depth. The green arrows show velocity relative to the fault, so that velocity is zero beneath the fault. Motion is in opposite directions on either side of the fault, and far from the fault the velocity is at plate rate. (D) No crustal anisotropy. There are no stresses from plate boundaries and no mantle-buoyancy anomalies. The mantle lithospheric strength isolates the crust from the sub-horizontal asthenospheric flow (i.e. mantle flow does not drive Couette flow in the mid-lower crust). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 5.1 Regional map of southern California and the Continental Borderland. The black rectangle marks the area where the three petroleum industry surveys were deployed and the inset map shows the distribution of the instruments (with the Long Beach array in red, the Extended Long Beach in green, and the Seal Beach array in blue). The stations of the Southern California Seismic Network (SCSN) are shown as inverted triangles and the regional mapped faults are delineated with thick gray lines. Seismic station LAF is shown as an inverted red triangle. . . . 71 xxi 5.2 Wavefield snapshots of the ambient noise correlation functions from a single virtual source of the Long Beach array (left), Extended Long Beach array (middle), and the Seal Beach array (right). Clear ballistic waves propagating away from all three of the sources can be observed. The number of instruments and time of operation of each survey is presented on top of each panel. 72 5.3 Stacked gather of the ambient noise cross-correlations and inversion for an average 1-D shear-wave velocity model. (a) Shows the first 30 seconds of the stacked cross-correlations from all three surveys. The bin-size of the spatial stacking is of 100 meters. The waveforms are band-passed filtered between 0.5 and 4 Hz. Clear fundamental (FM) and first overtone (FO) surface waves traveling at velocities between 0.3 and 1.5 km/s can be observed for the entire distance range. (b) Shows the average FM and FO dispersion curves that are obtained through the slant stack analysis of each survey’s dataset. The continuous red lines represent the predicted dispersion curve that results from our inversion process. (c-d) Shows the sensitivity kernels for the FM and FO, respectively. (e) Shows a comparison between shear wave velocity profile that is obtained from jointly inverting dispersion curves in (b) (black line) and a profile that is extracted from the CVM-S4 model (Shaw et al., 2015) at the geographic center of the nodal stations (red line). The black dashed line on top of panels (c-e) represents the free surface (FS). . . . . . . . . . . . 73 5.4 Schematic representation of the neighborhood-based cross-correlation method for phase arrival picking. Receivers that are less than one wavelength apart (<1λ) from a reference station (R1 ) are grouped together to form a subarray. The differential delay times, ∆t, between the reference waveform and all the other waveforms in the subarray are measured through waveform crosscorrelation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 xxii 5.5 Example of relative phase delay measurements for a virtual source of the Long Beach array. (a) Shows the L-curve analysis that was performed to determine the strength of the regularization term, S, in the inversion of the relative traveltimes. (b) Shows how the resulting relative traveltimes of the 1-s period waveforms vary as a function of distance from the virtual source. The dashed red line on top of the measurements has a corresponding slope of 0.70 km/s. (c) Shows the spatial distribution of the relative traveltimes that are in (b). The contours are drawn in time increments of 1 second. The red star denotes the location of the virtual source and the thick dashed circle around the source marks the three-wavelength minimum distance criterion that was imposed in our analysis. (d) Shows the process by which we extend the differential delay measurements to higher frequencies. The waveforms on the top are two ambient noise cross-correlations functions that are filtered at 1-s period and windowed around the expected time of arrival of the ballistic surface waves. The source-receiver raypath of these two waveforms is shown in (c). The bottom panel shows the correlation function of the two traces as a function of period. To determine the differential delay time of these two waveforms across the entire frequency range, we use the 1-s measurement to track down the phase delay to higher frequencies (red dashed line). . . . . . 77 xxiii 5.6 Example of phase velocity measurements for a virtual source of the Long Beach array at 1-s period. (a) Shows the source and receiver configuration that is used to make a single phase velocity measurement at a particular site. The yellow star marks the reference station that is used to form the receiver subarray and the green circles correspond to all the stations that were grouped together using the half-a-wavelength maximum distance criterion. The aperture of the subarray is approximately one-wavelengthlong. (b) Shows the relative traveltime of the receiver subarray in (a) as a function of their distance from the virtual source. Each marker is color-coded by the amount of times that a differential measurement of that station was used to construct the relative traveltime surface. The red circles depict the predicted relative arrival times that result from our time-based beamforming. (c-d) Show the distribution of our bootstrapped backazimuth and velocity measurements for the data shown in (b). (e-g) Shows the spatial distribution of all the velocity, backazimuth, and error measurements that were obtained for the same virtual source. . 5.7 Schematic representation of the construction process of isotropic phase velocity maps. (a-c) Show the phase velocities, backazimuths, and error measurements for different virtual sources of the Long Beach array at 1-s period. The errors in (c) are used to weight the averaging of the velocities in (a) to ultimately construct a frequency-dependent isotropic velocity map. The backazimuth measurements in (b) are reserved and used at a later step to characterize the anisotropy across the array. (d) Shows, from left to right, the isotropic phase velocity maps for the 1-, 0.8-, 0.6-, and 0.40-s Rayleigh waves. Each of these maps are generated in the same fashion, with the exception that the relative phase delay measurements of the higher frequencies were derived from the 1-s period waveforms. For periods larger than 1-s, such measurements were done independently. . . . . . . . . 79 80 xxiv 5.8 Comparison of our 1-s isotropic phase velocity map (left) with a reflective time slice of the active part of the survey at 0.40 s (right). Despite their difference in nature, the two different quantities in these images appear to illuminate similar tectonic features. The reflection image is courtesy of 3D Seismic Solutions. 81 5.9 Inversion process of the regionalized dispersion curves. The maps on the left show how the phase velocities vary as a function of period. The black dashed black lines on the right panels correspond to the dispersion curves that are extracted at locations A and B of the maps. The red continuous lines correspond to the dispersion’s best-fit. The number of updates that were required for the inversion to converge in these two examples is shown on the top-right of each panel. . . . . . . . . . . . . . . . . . . . . 82 5.10 Depth slices and cross-sections of our velocity model. (a) Shows a depth slice our velocity model at 50-, 100-, and 600-m depths. (b) Shows four cross-sections of our velocity model plotted on top of the migrated reflection images that were generated by the petroleum company that carried out the Seal Beach survey. To highlight the amount of lateral variation in each of our profiles, we removed the average velocity at each depth. The black dashed lines are used to mark prominent discontinuities in the migrated images. The location of each cross-section is shown in the leftmost map in (a). The migrated reflection profiles are courtesy of 3D Seismic Solutions. . . . . . . . . . . . . . . . . . . . . . . 83 5.11 Example of the Rayleigh wave azimuthal dependence at 1-s period for two different locations at the Long Beach array. The color of the point cloud represents the density of the velocity measurements and the black curve marks the predicted directionality using the first three terms in Equation 5.6. The best fitting coefficients, along with their 95% confidence intervals, are shown in the bottom right of each panel. . . . . . . . . . . . . . 85 xxv 5.12 Anisotropy measurements for the 1-s period Rayleigh wave for the Long Beach array. The bars point to the fast direction of anisotropy and their length is proportional to the magnitude. The colored markers behind each anisotropy vector denote the isotropic velocity (C0 term in equation 5.6). The continuous red line on the right map marks the mapped surface trace of the regional faults, whereas the dashed red lines delineate sharp changes in the anisotropy measurements that we suggest might be indicative of the presence of faults or other geologic boundaries. We associate these structures to the following: (a) NIF zone; (b) Silverado aquifer; (c) Compton-Los Alamitos Fault; (d) Garden Grove Fault; (e) Prograding clinoforms (see reflection image in Figure 5.8). The yellow star marks the location of Signal Hill. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.13 Wavefield snapshots of the ambient noise correlation functions using station LAF of the SCSN as a virtual source. Clear fundamental mode (FM) and first overtone (FO) surface waves can be observed. The dashed lines are the same as in Figure 5.12 and delineate the zones where the anisotropy measurements exhibit a drastic change in direction. The wavefield was band-passed filtered between 0.5 and 1.5 Hz. The time-stamp of each frame is shown in the upper right side of the panels. . . . . . . . . . . 5.14 Wavefield snapshots of the ambient noise correlation functions using station LAF of the SCSN as a virtual source. Clear fundamental mode (FM) and first overtone (FO) surface waves can be observed. The dashed lines are the same as in Figure 5.12 and delineate the zones where the anisotropy measurements exhibit a drastic change in direction. The wavefield was band-passed filtered between 0.5 and 1.5 Hz. The time-stamp of each frame is shown in the upper right side of the panels. . . . . . . . . . . 86 87 88 xxvi 5.15 Map of relative amplification for the source configuration shown in Figure 5.14. These values describe the expected PGA relative to the one that is obtained from a 1-D model. The red star marks the location of the source. The right panels show the synthetic waveforms that are recorded at locations A and B of the map for both our 3-D model (black) and a 1-D model (red). Note how the structure of the NIF generates elongated zones of intense amplification that result in ground motion changes that vary for several factors on the sub-kilometer scale. The waveforms are low-passed filtered at 4 Hz, which is the maximum frequency that is resolvable in our simulations. . . . . . . . . . . . . . . . 5.16 Experiment set-up for the estimation of the amplification factors for the case of a vertically propagating shear wave. The source wavelet, and its frequency content, are shown in the upper right panel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.17 Map of relative amplification for the case of a vertically propagating shear wave. These values describe the expected PGA relative to the one that is obtained from a 1-D model. . . . . . 5.18 Comparison between the Vs30 and the velocity measurements of this investigation. (a) Shows the Vs30 values across the surveys. (b) Shows our shear wave velocity estimates at 15-m depth. (c) Shows the ratio of (a) and (b). (d) Shows the sensitivity kernels of our measurements for the frequency range of analysis. The shallowest kernel peaks at around 30-m depth (red circle). . . . 6.1 Regional map of southern California and the Continental Borderland. The black rectangle marks the location of the Long Beach array and the inset map shows the distribution of the seismic sensors (as black circles). The mapped faults are from Jennings and George, 1994. The main trace of the Newport Inglewood fault is labeled as NIF. The major strands of the Newport Inglewood fault are denoted by red lines in the magnified map. To provide some context on the level of density of the Long Beach array, the stations of the Southern California Seismic Network are plotted as inverted red triangles. . . . . . . . . . . . . . . . 89 90 91 92 101 xxvii 6.2 Stacked record section of the Long Beach virtual shot gathers. The waveforms are band-passed filtered between 2 to 8 Hz. The bin size of the spatial stacking is of 50 m. Clear fundamental mode Rayleigh waves, some combination of higher-mode Rayleigh waves and S-waves, and P-waves are visible at almost all offset ranges. The dashed black lines mark the different time windows in which each of these phases are expected to arrive (from Lin et al., 2013). . . . . . . . . . . . . . . . . . . . . . . 103 6.3 Observed and synthetic wavefields windowed around the expected P-wave time of arrival. (A) Record section of the first 5 seconds of the stacked cross-correlations. The thick black lines mark 3 different reference velocities that approximate the variations in the apparent wave speed of the first-arrival P-wave as a function of offset. (B) Synthetic wavefield computed by 2D finite difference modeling using the inverted 1D velocity model (black profile in Figure 6.4B). The frequency range of the waveforms is the same as Figure 6.2. (C) Snapshots of the 2D velocity wavefield at different times (1.5 s, 2.5 s, and 3.5 s). The explosion at the surface marks the location of the source. Note that attenuation was not considered in the wavefield simulation. . . 104 6.4 T au-p representation of the first 5 seconds of the stacked crosscorrelated waveforms (A) and inverted velocity model (B). As a comparison, the model that is obtained from the inversion of the T au-p curve (black line in A) is shown together with a 1D profile of the CVM-S4 model (Shaw et al., 2015) at the center of the Long Beach array. . . . . . . . . . . . . . . . . . . . . . . . . . 105 xxviii 6.5 Examples of traveltime picking of individual traces. (A) Map of the Long Beach array configuration. The red lines mark the shortest path between the virtual source (station S) and the receivers (stations R1 and R2) involved in the computation of the correlograms shown in (B) and (C). (B) Cross-correlated waveforms between stations S and R1 together with its time-frequency spectra computed via the S-transform. The dashed red line on top of the waveforms marks the predicted arrival time of the Pwave. The star markers in the lower panel depict all of the energy peaks identified in the time-frequency spectra and the red box delimits the time window used to identify the potential P-wave arrivals. The energy peaks in the spectrogram are determined using a watershed algorithm (Meyer, 1994). (C) Same as (B) but for the correlogram between stations S and R2. . . . . . . . 107 xxix 6.6 Double-beamforming process for a single source-receiver array configuration. (A) Map showing the entire Long Beach survey together with a source array (red circles) and a receiver array (gray circles). The red radius around the source array marks the 1-km distance threshold imposed in our stacking scheme. For this particular source-receiver configuration, the source radius is of 0.32 km, the receiver radius is of 0.42 km, and the distance between the arrays’ geographic centers is of 2.12 km. (B) Entire source and receiver side velocity space of the DBF. The yellow star indicates the best-fitting velocity pair obtained in the leastsquares inversion. The black dots correspond to the best-fitting velocity pairs obtained in the bootstrapping process and the red box marks their 95% confidence limits. For this case, the area of the 95% confidence region is 0.09 km2 /s2 . (C) Empirically derived probability density function of the source-side velocity (SSV) and the receiver-side velocity estimations (RSV) calculated from the bootstrapped results. (D) Beamed trace (black waveforms) plotted behind one of the 900 original correlograms used to build the beam (red waveforms). Note how the P-wave is essentially unidentifiable in the original correlogram and how its signal-to-noise ratio improves dramatically after applying DBF. The retrieved phase corresponds to a diving P-wave propagating between the center of the source and receiver arrays. . . . . . . 6.7 Wavefield emitted by a virtual source. The location of the virtual source is marked by the yellow star. The lag-times are shown on the upper left corner of each panel. The top panels show snapshots of the 2-8 Hz band-passed wavefield observed at each station. The middle panels show snapshots of the beamed wavefield observed at each station. Strong body waves are now visible propagating away from the virtual source. The lower panels shows snapshot of a synthetic wavefield using the inverted velocity model (black profile in Figure 6.4B). Note how the observed wavefield is not completely spherical as in the synthetic case. . 109 111 xxx 6.8 Examples of P-wave refractions from a N-S line (A) and an E-W line (B). The profiles are formed by stacking all available beams along corridors that are 0.5 km wide (green areas in C) in a 50-m offset bin. The location of the virtual sources are marked by the yellow stars. Note how there exists a clear jump in the N-S profile that is coincident with the Newport Inglewood fault (thick red line in all panels), and shorter jump in the far end of the E-W profile that might be associated with a smaller fault (thick blue line in panel B). . . . . . . . . . . . . . . . . . . . . 6.9 Stacked record section of beams (A) and waveform decomposition process of a single reference trace (B-D). (A) Stacked recordsection showing how the DBF has successfully isolated the body wave energy and suppressed any other arrival that is not arriving with a slowness of a refracted P-wave. The dashed colored lines mark the time window where the higher-mode surface waves were present before the DBF. (B) Reference trace (black waveform) and reconstructed trace (red waveform) that is obtained after linearly adding all of the decomposed waveforms (all traces in D). The reference trace corresponds to the stacked beam in the 7 km offset bin (red arrow in A). (C) Time-frequency spectra of the reference trace computed via the S-transform. The star markers indicate the arrivals identified in the power spectrum and the numbers indicate their order in time. (D) Decomposed traces sorted by time. A total of 27 apparent phases are identified in the reference trace spectra, with the diving P-wave being the most prominent amongst all. The green trace alone is what is used as a template to pick the refracted P-wave arrival in individual beams that have a 7-km offset. . . . . . . . . . . . . . 112 113 xxxi 6.10 Traveltime measurements for two virtual sources localized on opposite sides of the survey. The left panels show the distribution of the traveltime differences between the picked time and the theoretical arrival. The middle panels show the normalized correlation coefficient between each picked phase with their respective waveform template. The right panels show the signal-to-noise ratio of the beams. Here, we define the SNR by the ratio between the peak amplitude within a window containing the refracted P-wave arrival to the root-mean-square of the amplitudes before this window. The locations of the virtual sources are marked by the yellow stars. The red lines depict the major strands of the Newport Inglewood Fault and the black dashed lines delimit the regions where traveltime measurements appear to be random. . 6.11 Horizontal (A) and vertical (B) slices of the inverted velocity model. Velocities are presented as perturbation (in percent) from the horizontally averaged inverted model. To emphasize the lateral velocity variations, each horizontal slice is plotted with its independent scale bar. The depth of the horizontal slices is given at the upper right corner of each panel (0.14 km, 0.49 km and 0.98 km). The blue lines in the leftmost panel in A show the location of the vertical slices in B (P1, P2, P3 and P4). For this set of plots, the local coordinates have been shifted to the center of the survey. The red lines in every panel mark the major strands of the Newport Inglewood fault. . . . . . . . . . . . . . 6.12 Depth slices of the velocity uncertainties obtained after bootstrapping the virtual sources 50 times. The depth of the horizontal slices is given at the upper right corner of each panel (0.14 km, 0.49 km, and 0.98 km). . . . . . . . . . . . . . . . . . . . . 115 116 117 xxxii 6.13 Fault-perpendicular cross-sections of the inverted velocity model. (A) Map showing the Long Beach survey with the location of the profiles. (B) Cross-sections illustrating the general structural pattern across different segments of the Newport-Inglewood fault zone (modified from Wright, 1991). The thick back lines in A show the location of the two profiles (X-X’ and Y-Y’). (C-G) Vertical cross-sections of the inverted velocity model (colored profiles in A). The red arrows on top of each profile indicate the location of the surface trace of the Newport Inglewood fault. The black dashed lines in each profile mark our interpretation of the fault. The inferred fault geometry of Wright (1991) across profile X-X’ is overlain on profile E. . . . . . . . . . . . . . . . 119 6.14 Wavefield emitted by the Carson Mb2.5 earthquake at 4 seconds after its origin time (A). Note how the earthquake surface waves clearly break up along the Newport Inglewood fault. The dashed black box in the zoomed panel is the same as in Figure 6.11A and delimits a region where the wavefield is being most advanced and stretched by the high velocities of the fault. The velocity perturbation at 0.98 km depth of this small area is shown in (B). 120 7.1 Schematic diagram showing the ray paths of the direct P-wave, the depth phases (pP, sP) and the first four water multiples (pwn P, swn P) (left). Synthetic velocity waveforms for a 20 km deep strike-slip earthquake under a homogeneous water layer (right). For simplicity, the synthetic seismogram is shown without water phases (top) and with water phases (bottom). The water multiples alone are presented in the middle. The waveforms were calculated for a station at a distance of 84◦ using the Preliminary Reference Earth Model (PREM; Dziewonski and Anderson, 1981) with a 3-km-thick water layer on the source side. Note how the ray parameter of the water phases is the same as the one of the direct P-wave. . . . . . . . . . . . . . . . . . . . 130 xxxiii 7.2 Example illustrating the effect of seafloor roughness on teleseismic waveforms. The synthetic waveforms show a noticeable increase in amplitude and complexity as the ocean bottom roughness increases (A). The dashed green line on top of the elevation profiles depict the event’s epicenter. The high-resolution bathymetry profile used in the simulations is extracted from ETOPO1 global relief model along profile P-P’ in (B). The lowresolution bathymetry profiles are built from applying running averages of different lengths to the high-resolution one. Note that we only show the first 1000 km of each profile to accentuate the near-source structure but the models extend all the way to the inland station. Focal mechanism and takeoff angle of the direct P wave are shown in (C). . . . . . . . . . . . . . . . . . . 7.3 Waveform modeling results for the 2013 Mw6.1 Gofar earthquake using different structural models (see Figure 7.4 for the event location). Black traces are observed data and red traces are synthetics (0.01-0.2 Hz). Station network and name is shown in panel (A) only but is the same for all 4 panels. Synthetics in (A) are computed using a 1D Earth (Model 1). Synthetics in (B) are computed using a 2D Earth with a flat ocean bottom on the source side (Model 2). Synthetics in (C) are computed using a 2D Earth with realistic topography (Model 3). Synthetics in (D) are computed using the hybrid method with 3D source side structures (Model 4). The average waveform misfit of each set of synthetics is shown in the upper part of the panels. Examples of the velocity models (for the 1D and 2D simulations) and the bathymetry mesh used in the 3D simulations are shown on the right. The elevation data for Model 4 is extracted from the Global Multi-Resolution Topography synthesis (GMRT; Ryan et al., 2009). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 133 xxxiv 7.4 Topographic map of the Gofar transform fault system. G1, G2, G3 mark the different segments of the fault that are separated by intratransform spreading centers. Beachballs show the focal mechanisms and centroid locations from the Harvard GCMT catalogue. The circled numbers represent candidate locations for the 2002 Mw5.1, 2007 Mw6.1, and 2016 Mw5.7 earthquakes. The average intersource spacing between neighboring candidate locations is 5 km. The small inset map is a zoom of the 2016 Mw 5.7 event epicentral region. The dashed purple line marks the extension of a potential sub-parallel fault (from surface bathymetry). The red circled numbers mark the centroid locations determined by this study. Elevation data is extracted from the Global MultiResolution Topography synthesis (GMRT; Ryan et al., 2009). . 7.5 Waveform modeling results for the 2002 Mw5.1 (A), 2007 Mw6.1 (B), 2008 Mw6.0 (C), and 2016 Mw6.1 (D) Gofar earthquakes using the reported GCMT solutions and a 5-km source depth. Black traces are observed data and red traces are synthetics (0.01-0.2 Hz). Station name and network are shown above each trace pair. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.6 Record section of 3D synthetic waveforms generated for one teleseismic station using source locations 7-27 in Figure 7.4 (A). All waveforms have been band-passed between 0.01 to 0.2 Hz. The blue shaded region delimits the traces for the fault-perpendicular sources whereas the yellow shaded region delimits the traces for the fault-parallel sources. The mean azimuth and distance to the receiver is of 0.70◦ and 51.76◦ , respectively. Beamforming output for the 3 time windows that are enclosed by the red dashed lines in the record section (B). Slow and prominent arrivals can be observed being generated from multiple azimuths. Collection of all maximums in every beamforming output for the full station synthetic array (C). The size of the markers represents the beam power whereas the color represents the apparent wave speed. The presence of slow arrivals concentrated in the azimuth along the fault’s strike suggests that the fault itself is behaving as a strong scatterer. . . . . . . . . . . . . . . . . . . . . . . . . 135 136 137 xxxv 7.7 Waveform modeling results for the 2002 Mw5.1 using source locations 1-6 in Figure 7.4. Black traces are observed data and red traces are synthetics (0.01-0.2 Hz). Station name and network are shown above each trace pair in the first panel only. The green contoured panel shows the waveform comparison of the observed and synthetic traces for the location that results in the best fit (location 1) whereas the the blue contoured panel shows the waveform comparison for the location that is closest to the GCMT solution (location 5). . . . . . . . . . . . . . . . . . . . 7.8 Illustration exemplifying the beamforming process applied to the 2007 Mw6.1 recordings (A). The map shows the stations used to produce a beamed trace at station CI.ARV (red inverted triangle). The blue inverted triangles represent the 19 closest stations used to form the subarray. The inset plot shows the beamforming output for a time window around P-wave arrival of the 20 traces. The subarray and beam information is shown in the bottom white box. The record section in (B) shows the raw waveforms of the 2007 Mw6.1 event at all available US Array stations (gray inverted triangles in A). The record section in (C) shows the waveforms at the same stations after applying the spatiotemporal beamforming. Both sets of records have been aligned with the P-wave onset time and band-pass filtered between 0.01 and 0.2 Hz. The red dashed lines in (B) and (C) mark the Pwave and PP-wave arrival time. The green arrows in (C) mark the arrival time of two prominent water phase arrivals. . . . . . 138 139 xxxvi 7.9 Waveform modeling results for the 2007 Mw6.1 earthquake. Black traces are observed data and red traces are synthetics (0.01-0.2 Hz). The leftmost panel shows the synthetic waveforms of earthquake sources 7-17 in Figure 7.4 for station CI.SWS. The black waveform plotted on top corresponds to the observed beamed trace around the same station. The top green arrows mark the arrival time of the two prominent water phases arrivals. Note how source locations 13-16 are the only ones capable of producing the first prominent water phase arrival (blue-circled numbers). The three rightmost panels correspond to the waveform simulation results for a strike-slip source placed at the GCMT reported location at 5 km depth, at source location 15 at 5 km depth, and at source location 15 but at 15 km. Station name and network are shown above each trace pair in the first panel only. The green contour marks the waveform fit of our preferred solution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 7.10 Search space for the 2016 Mw5.7 event optimal location. The circles depict the earthquake candidate locations (see magnified map in Figure 7.4) and the traces show the waveform fit between observed (black) and synthetics (red) at station TA.L02F for 4 different locations (0.01-0.2 Hz). The color of the markers represent the waveform misfit for the same station. The waveform modeling results for this earthquake best-fitting location are shown on the right. Station name and network are shown next to each trace pair. The dashed purple line marks the extension of a potential sub-parallel fault (from surface bathymetry). . . 143 xxxvii 7.11 Heatmap showing the waveform sensitivity to bathymetry resolution as a function of frequency. The thick lines depict the P-wave wavelength in the water (λ=1 equals to one wavelength assuming a constant water velocity of 1.5 km/s). The white dashed lines mark the 1-arc minute base resolution of current global bathymetry models, the resolution limit of the resampled SRTM30 global bathymetry model (Becker et al., 2009), and the the resolution limit of the resampled SRTM15+ global bathymetry model (Tozer et al., 2019). The thick black arrow marks the bathymetry resolution for the Gofar transform fault (61 m; Ryan et al., 2009). The topographic profile used in these calculations is cut from the EPR to the center of the United states (from -4.75◦ /-105.40◦ to 41.38◦ /-119.17◦ lat/lon). The misfit is quantified via the shape-fitting criteria (Equation 7.1) in which the record generated in the highest resolution run is taken as the "true" waveform. Synthetics are calculated for a station at a distance of 48◦ using a 15-km-deep strike-slip source. 145 7.12 Waveform comparison of two similar-sized pairs of earthquakes occurring in a fast spreading oceanic environment (top) and a slow spreading environment (bottom). Both pairs of events are separated by the same distances, as determined by the time shifts of their fundamental-model Rayleigh waves for a full azimuth range. Event names and focal mechanisms are shown on the left size of each panel. Station name, distance, and azimuth are shown above each pair of waveforms. Distances and azimuths are approximate and calculated from the center of the fault. Fault name and slip rate are shown in the top part of each panel together with the estimated intersource distance. The normalized cross-correlation coefficient of each pair of waveforms is shown in blue next to the traces. Note that, despite having the same intersource distances, the water phases for the events occurring in the fast spreading ridge are more similar than the ones occurring in the slow spreading ridge. All waveforms have been band-passed between 0.01 to 0.2 Hz. For each panel, the focal mechanism and waveforms are color-associated. . . . . . . . . . 146 1 Chapter 1 INTRODUCTION The understanding of the Earth’s structure and dynamics is one of the main scientific targets of geophysics. To date, some of the most compelling contributions in this area have come from the observation and study of the seismic wavefield. In principle, when a mechanical disturbance (either man-made or natural) occurs, seismic waves are generated and start to propagate away from the source in every possible direction. These waves are recorded by seismographic instruments that are deployed all over the world, and the resulting time series (i.e. seismograms) are analyzed to extract information about Earth’s structure and the nature of the source. In this thesis, I center on the characterization of seismic wavefields to gain a deeper understanding of the elastic properties of the Earth, and use this information to solve geologic puzzles. In detail, I present five case studies, where I correlate the Earth’s background vibrations to extract coherent waves traveling between pairs of stations, and use the differences in their propagation properties to construct high-resolution images of the structure that they sampled. Although the principle in all five studies is similar, the focus and, ultimately, contribution of each investigation is unique. Below, I provide a brief introduction to each case study and summarize their motivations. In Chapter 2, I image the lithospheric structure beneath the Trans-Mexican Volcanic Belt (TMVB). The TMVB is a prominent and enigmatic feature of the subduction system in Mexico. Its volcanic style diversity and oblique orientation to the trench are explained by the large along-strike variations in the subduction parameters of the Rivera and Cocos plates. However, the abrupt termination of the TMVB on its eastern end with the Pico de Orizaba volcano is puzzling as the current slab model suggests that the transition of the Cocos flat-slab geometry to normal subduction is smooth through this region. There is evidence that suggests that a tear in the slab might be developing, but it is unclear how this feature can support the unusually large topographic gradient that connects the volcanic high peaks with the Veracruz basin just south of the volcanic front. To provide further insight into the transition anatomy of this portion of the slab, I construct a detailed and unified model 2 of the structure of the crust and uppermost mantle of south-central Mexico. In Chapter 3, I image the dynamics of the Wallowa Mountain block in northeastern Oregon. It is a generally understood that buoyancy anomalies within the Earth’s mantle create large convective currents that control the evolution of the lithosphere. Although there is observational evidence for this relation, the scale and mechanism by which mantle processes influence near-surface processes remain obscure. To shed light on this matter, I place constraints on the seismic properties of the lithosphere beneath the Wallowa mountains, where there is an over-thickened crust sitting on top of a major high-velocity anomaly that appears to be dangling in the upper mantle. In Chapter 4, I investigate the conditions under which mantle-based stresses affect the dynamics of the lithosphere. Inspired by the observations of Chapter 3, I use short-period ambient noise surface waves to resolve the lateral variations of crustal anisotropy in regions of the US where near-surface mantle-based deformation might be have occurred or is presently occurring. In particular, I focus around the California and Rocky Mountain areas since high-resolution tomographic images reveal the presence of dense upper mantle structures similar to the ones that I suggest to be controlling the evolution of the Pacific northwest lithosphere. In Chapter 5, I image the shear wave velocity structure of Long Beach, California. With the development of portable seismic instrumentation, the deployment of dense geophone surveys has become a common practice in the field of crustal geophysics. However, because of their financial cost, only a few regions around the world have hosted this type of experiments. This limitation has resulted in a lack of knowledge on both the scale of structural variations and the spatial variability in shaking intensity that we can expect on different types of tectonic settings. To aid in this issue, I use ambient noise data recorded at three dense petroleum industry surveys deployed at Long Beach, California, to produce a high-resolution model of the top kilometer of the crust. In Chapter 6, I image the compressional wave velocity structure at Long Beach, California. As it will become evident in Chapters 2-5, the analysis of ambient noise surface waves is a well-established routine in the seismic imaging community. However, the reconstruction of body waves from the cross-correlation of random wavefields has only recently emerged as an alternate and promising way to probe the elastic structure of the Earth. To test the feasibility of 3 their use, I extract body waves from ambient noise data recorded at a dense petroleum industry survey and generate a high-resolution model of the top two kilometer of the crust. In addition to the five case studies listed above, this thesis includes one last chapter (Chapter 7) in which I do not investigate any seismic properties of a particular region. Rather, I present a new method of analysis that shows great promise as a new way of characterizing the earthquake behavior of oceanic transform faults (OTFs). This issue is of great interests to the seismological community as OTFs are amongst the simplest tectonic environments on Earth and, as such, are one of the most attractive locations to investigate how slip is accommodated within the crust. However, because they tend to grow in the deep ocean, far away from where local seismic observations are available, the retrieval of accurate OTF source parameters represents an immensely complicated task. This limitation has, to date, hindered our ability to characterize the elastic structure of mid-oceanic environments and their tectonic behavior on a global scale. 4 Chapter 2 THE ENIGMATIC TRANS-MEXICAN VOLCANIC BELT This chapter was adapted from: J. C. Castellanos, R. W. Clayton, and X. Pérez-Campos (2018). “Imaging the eastern Trans-Mexican Volcanic Belt with ambient seismic noise: Evidence for a slab tear.” In: Journal of Geophysical Research 123.9, pp. 7741-7759. doi: 10.1029/2018JB015783. 2.1 Abstract The eastern sector of the Trans-Mexican Volcanic Belt (TMVB) is an enigmatic narrow zone that lies just above where the Cocos plate displays a sharp transition in dipping angle in central Mexico. Current plate models indicate that the transition from flat to steeper subduction is continuous through this region, but the abrupt end of the TMVB suggests that the difference in subduction styles is more likely to be accommodated by a slab tear. Based on a high-resolution shear wave velocity and radial anisotropy model of the region, we argue that a slab tear within South Cocos can explain the abrupt end of the TMVB. We also quantify the azimuthal anisotropy beneath each seismic station and present a well-defined flow pattern that shows how mantle material is being displaced from beneath the slab to the mantle wedge through the tear in the subducted Cocos plate. We suggest that the toroidal mantle flow formed around the slab edges is responsible for the existence of the volcanic gap in central Mexico. Moreover, we propose that the temperature increase caused by the influx of hot, less-dense mantle material flowing through the tear to the Veracruz area may have significant implications for the thermomechanical state of the subducted slab, and explain why the intermediate-depth seismicity ends suddenly at the southern boundary of the Veracruz basin. The composite mantle flow formed by the movement of mantle material through the slab tears in western and southern Mexico may be allowing the Cocos plate to rollback in segments. 5 2.2 Introduction The Trans-Mexican Volcanic Belt (TMVB) is one of the largest volcanic arcs on the North America plate. This Neogene continental arc grows over the central Mexican margin of the North-American plate, as a result of the subduction of the Rivera and Cocos plates along the Middle-America Trench (MAT), and is comprised of nearly 8,000 igneous structures that extend from the coast of Jalisco to the Gulf of Mexico in Veracruz [1] (Figure 2.1). The remarkable compositional variability of the TMVB and its unusual lack of parallelism to the trench are unique aspects of the Central America tectonics that have been extensively debated over the years. Today, however, geophysical evidence indicates that the Rivera and Cocos oceanic plates are being subducted with highly variable dip angles that explain the obliquity of the volcanic arc [2]. In the present-day configuration, the western sector of the TMVB overrides the Rivera micro-plate, whereas the central and eastern sectors are underlain by the Cocos plate. Slab contours determined by precise hypocentral localization show that the Rivera plate dips at a steep and constant angle of ∼50◦ beneath the Jalisco block and that earthquakes extend to a depth of ∼120 km [3]. Conversely, in central Mexico, the Cocos plate exhibits marked changes in its geometry, and seismicity does not exceed ∼80 km depth [4]. Beneath the Michoacán area (from 103◦ to 102◦ W), the Cocos plate displays a shallow slab dip that gradually decreases eastward under Guerrero (from 102◦ to 98 ◦ W), where the slab is sub-horizontal for about 250 km [3]. From 98◦ W, the slab dip increases eastward until reaching an uniform angle of ∼50-60◦ in central America (90◦ W) beneath Guatemala, El Salvador, Nicaragua, and Costa Rica [5–7]. More recent studies have refined the morphology of the Rivera and Cocos plates beneath the continent and investigated their relation with the particular space distribution of volcanoes along the convergent regime. Waveform modeling of moderate-size intra-slab earthquakes recorded by the MARS (Mapping the Rivera Subduction Zone) experiment [9] and seismic anisotropy measurements [11] suggest that the Cocos slab could be currently fragmenting along the landward projection of the Orozco Fracture Zone (OFZ) by a similar process to that which occurred when the Rivera plate separated from the proto-Cocos plate [12]. This tear is proposed to indicate the ongoing fragmentation of the Cocos slab into North Cocos and South Cocos plates and would explain the 6 Figure 2.1: Tectonic setting of central and southern Mexico. The brown box indicates the main focus area of this study. The black contour lines depict the depth of the subducted slab compiled using the results of receiver functions, tomography studies, and hypocentral re-localization [from 8]. The dashed part of these contour lines indicates the transition between the flat and normal dipping portion of the South Cocos plate. The inverted white triangles show the location of the broadband seismic stations that have operated within this region. Regional seismicity at different depths, as reported by the Servicio Sismológico Nacional (SSN) since 1998, is shown as color-coded circles. The volcanic provinces: the Trans-Mexican Volcanic Belt (TMVB), the Anegada High (AH), the Los Tuxtlas Volcanic Field (LTVF), and the Modern Chiapanecan Volcanic Arc (MCVA) are delimited by the orange areas, and the location of the main stratovolcanoes and calderas are marked by the red triangles and red circles, respectively. The light purple region delimits the Veracruz basin (VB), and the light magenta region delimits the Jalisco Block. The green dashed line denotes the Rivera-Cocos plate boundary, and the purple dashed line denotes the projected path of the Orozco Fracture Zone (OFZ) beneath North America [9]. The blue dashed line indicates the location where the slab exhibits a significant change in dip and where Dougherty and Clayton [10] propose the existence of a trench-perpendicular tear. The ages of the plates (Ma) and the convergence rates (cm/yr) along the Middle America Trench (MAT) are shown in blue and red numbers, respectively. The inset map shows this work study area and the distribution of every station used in this study. 7 observed offset of the volcanic arc in the Michoacán area [9]. East of this region, receiver functions and seismic velocity tomography along the MASE (MesoAmerican Subduction Experiment) line confirmed that the Cocos plate is sub-horizontal in the Guerrero region, and showed that the slab then plunges steeply into the mantle at a dip of ∼75◦ where it is finally truncated at 500 km [13–15]. This atypical subduction geometry suggests that the slab is rolling back and forcing the volcanic arc to retreat [13], as confirmed by the age progression of volcanism migrating trenchward [16]. Further to the south, detailed analysis of receiver functions along the VEOX (Veracruz-Oaxaca) line indicates that the Cocos plate dips at a constant angle of ∼26◦ to a depth of 150 km, where it appears to be truncated by an anomalous south-west dipping slab that extends to a depth of 250 km [17, 18]. In between the MASE and VEOX experiments, near the transition from shallow to normal subduction, the TMVB ends abruptly with Pico de Orizaba. This stratovolcano is the highest point in Mexico and sits at the front of an active volcano chain that is oriented almost perpendicular to the trench. Coupled with this feature is an extremely sharp topographic gradient that connects the volcanic arc high peaks with the Veracruz basin (altitude drops ∼5,000 m over a horizontal distance of just 120 km) (Figure 2.2). Crustal thickness measurements derived from gravity data [19, 20] and receiver functions [21] indicate that a simple isostasy compensation model is insufficient to explain the thickness difference between these two regions. Furthermore, the presence of the nearby late Miocene Anegada High submarine volcanic complex [22] and the active Los Tuxtlas volcanic field (LTVF) [23] marks an interruption of arc volcanism that is thought to be associated with the steepening and rollback motion of the slab. However, the mechanisms of their origin remain unclear. Altogether, these features indicate that the transition from shallow to normal subduction in central Mexico is more complicated than it was originally suggested in the Pardo and Suarez [3] plate model. To account for the abrupt end of the TMVB, and the discontinuity of arc volcanism in central-southern Mexico, Dougherty and Clayton [10] propose the existence of a possible tear located within the downdip portion of the South Cocos plate. Conversely, patterns of local seismicity, tectonic tremor, and slow-slip events support the idea that a sharp flexure of the slab in both the downdip and along-strike directions accommodates the transition from flat to steeper subduction [24]. In this study, 8 Figure 2.2: Elevation map of the eastern section of the Trans-Mexican Volcanic Belt showing the abrupt termination of the volcanic arc. Note how altitude drops more than 5,000 m from Pico de Orizaba to the coastal Veracruz basin over a horizontal distance of 120 km. This region corresponds to the one encompassed by the brown box in Figure 2.1. we determine the shear wave velocity structure and radial anisotropy of the upper crust and lithosphere in central Mexico by inverting Rayleigh and Love dispersion curves. In addition, we quantify the azimuthal anisotropy present beneath each seismic station using an array analysis to resolve variations in the direction of mantle flow in the region. The combined understanding of the velocity distribution and seismic anisotropy allows us to place constraints on the transition structure of the subducted slab and its relation with the abrupt end of the TMVB at its eastern limit. 2.3 Data and Method The data used in this study consists of surface wave signals obtained from the three-component cross-correlation of background noise recorded at over 2,000 broadband stations. This dataset results from combining every available seismic network that operated within Mexico and its surroundings (from 5◦ to 40◦ N and -125◦ to -60◦ E) from January 2006 to December 2016. The reader is referred to Pérez-Campos et al. [25] and Córdoba-Montiel et al. [26] for a summary of the permanent stations in this area. Ambient Noise Cross-Correlations Both theoretical and experimental studies have demonstrated that by crosscorrelating the ambient noise recorded at two stations over a sufficiently-long period of time, the Green’s function between the two stations can be retrieved 9 [e.g., 27–31]. Here, we use a technique very similar to the one described by Bensen et al. [32] to compute the three-component cross-correlations of continuous recordings between all synchronous station pairs. The single-station data preparation consists of (i) down-sampling the records of all three-components to 1 sample per second and dividing them into 1-day time windows, (ii) removing the mean and trend value in each time window, (iii) band-pass filtering between the 3-100 s period band, (iv) whitening the spectra, and (v) normalizing in the time domain. Once the preprocessing is complete, each time window is cross-correlated, normalized to unit peak amplitude, and averaged over time. The cross-correlation traces are then rotated from the east-northvertical (ENZ) frame into the radial-transverse-vertical (RTZ) frame between all station pairs. To simultaneously determine all nine components of the Green’s tensor, we follow Muir and Tsai [33] and build the rotation matrix as M = M1 ⊗ M2 , where M1,2 are the three-component rotation matrices from the ENZ to the RTZ frame for the individual stations, and ⊗ is the Kronecker product. Finally, to enhance the signal-to-noise ratio (SNR), the causal and anticausal parts of the cross-correlations are stacked to obtain the so-called symmetric cross-correlations. This process gave rise to more than 136,000 noise correlation functions for each of the nine calculated components. As evidence of the effectiveness of the process described above, Figure 2.3 shows the ZZ, RR, ZR, and TT cross-correlations as a function of interstation distance for the period band of 5-100 s. Prominent surface wave signals and various body wave arrivals are visible in all four component combinations. Dispersion Measurements For the 3-70 s period band, we apply an automated image transformation technique [35] to retrieve the Rayleigh and Love wave phase and group velocity dispersions from the ZZ, RR, ZR, and TT cross-correlations. Three selection criteria are imposed before accepting a measurement at a certain period. First, the minimum inter-station spacing is set to one wavelength [e.g., 36]. Second, the maximum phase and group velocity measurement deviation from the average dispersion curve is limited to 0.5 km/s. Third, the SNR threshold is empirically set to 5. Here we define the SNR by the ratio of the peak amplitude within a window containing the surface wave signal to the mean of the noise trailing the direct arrival window. Generally, the SNR decreases with interstation distance due to scattering and attenuation; however, the re- 10 Figure 2.3: ZZ, RR, ZR, and TT cross-correlations as a function of interstation distance for the period band of 5-100 s. The traces are stacked by bins of 2 km and 0.5 s. As expected, the extracted Green’s functions show clear Rayleigh waves in the ZZ, RR, and ZR components and clear Love waves in the TT component; although several body wave arrivals emerge as well. Note the presence of the anomalous cS-cP core phase in the ZZ and ZR cross-correlations [34]. laxed restriction of one wavelength for the shortest path still allows sufficient coverage to extend our measurements up to 70 s. Finally, to increase the robustness of our measurements and ensure a larger period coverage, we stack the ZZ, RR, and ZR dispersion curves for identical paths [e.g., 37]. For this last step, we introduce an additional selection criterion that requires that the standard deviation at every common period be smaller than than 0.1 km/s. 11 Tomographic Inversion The Rayleigh and Love interstation dispersion measurements are used to invert for phase and group velocity maps using the method of Barmin et al. [38]. But, before formulating the tomographic problem, it is necessary to consider that the elongated shape of our study region limits the number of long period trench-perpendicular raypaths to the MASE and VEOX lines only. To mitigate the uneven ray coverage at long periods, we choose to include distant stations into the inversion and project all crossing seismic rays to a relatively wide rectangular area (from 14◦ to 23◦ N and -107◦ to -90◦ E) so as to provide an appropriate azimuthal coverage. A fundamental disadvantage of this approach however is the assumption that all velocity anomalies are contained within the study region. To that end, all measurements are down-weighted throughout the inversion according to the percentage of the ray that lies outside the target area. The size and extension of this area is chosen by trial-and-error until the inverted models show the minimal amount of smearing while preserving all first-order features when compared to velocity maps that are produced with wider boxes and coarser grids. To carry out the actual inversion, homogeneous maps are constructed on a 0.4◦ x 0.4◦ regular grid across the study region and are defined relative to the average slowness observed at each period of interest. The optimal grid size is chosen empirically based on a resolvability test for the periods with the lowest path density using Voronoi diagrams [39]. We then perform two full tomographic inversions. In the first inversion, no additional weights are added to the data and we apply heavy damping, which results in a highly smoothed model. Synthetic traveltimes are then computed and used to identify and discard highly anomalous measurements in the observations. A 3σ residual threshold is set for this criterion [e.g., 40]. For the second inversion, the remaining data are further down-weighted according to their misfits in the first inversion and the optimum smoothing factor is determined from the misfit and model smoothness trade-off curve. The resulting slowness maps are then converted to phase and group velocity maps. To assess the capability of the different ray path geometries to resolve for contrasting slowness distributions, we follow Ma and Clayton [41] and use the resolution matrix R = (G 0 C −1 G + Q) −1 G 0 C −1 G (where G is the generalized inverse or forward operator on the slowness model, C is the data covariance matrix, and Q is the regularization matrix) from the tomographic inversion to generate standard checkerboard resolution maps. For this test, we set up 12 input models containing ±1 km/s perturbations and evaluate how accurately the tomography is able to retrieve the anomaly distribution. To characterize the error of the model, we use the diagonal elements of the model covariance matrix Cmm = (G 0 C −1 G + Q) −1 ; which reflects the variance of the model subject to the variance of the data [41]. As an example, Figure 2.4 shows the velocity maps, checkerboard resolution maps and model error maps at 34 s period for both Rayleigh and Love surface waves. Perhaps the most important feature of these maps is that the velocity distribution varies for all four velocity types even though the measurements are taken at the same period. This can be explained by their differences in depth sensitivity and is the main reason why their joint analysis provides better constraints on the radial structure of the crust and lithosphere [e.g., 37]. Nonetheless, all models exhibit a similar pattern and reveal low velocities beneath the TMVB (Figure 2.5). Based on the resolution maps, we find that the checkerboard structure is satisfactorily reproduced in all inversions except in the southeast sector of the study area. As expected, error along the coast is high, but it decreases to values smaller than 0.1 km/s as we approach central Mexico, where the path coverage in denser. Inversion for Shear Wave Velocity and Radial Anisotropy Once the velocity maps are constructed for each frequency, we extract the velocity dispersion curves at each (x,y)-point in the grid and use a linearized inversion algorithm [42] to simultaneously map the phase and group velocities as a function of period to shear wave velocity as a function of depth. This process is done independently for Rayleigh and Love waves to obtain a VSV and a VSH model, respectively. For each inversion, we use a linearly increasing 1-D model discretized into 2km layers at the top 50 km, then 5-km layers to 100 km, and finally 10 km layers to 140-km depth, as a starting model to solve for a smooth structure without a Moho discontinuity. The implications of not imposing a Moho in the inversion are addressed in Ma and Clayton [41]. We then assume a constant Vp /Vs ratio of 1.73 for the whole structure and determine the density from the compressional wave velocity [43]. Throughout the inversion, both the Vp /Vs ratio and density remain fixed and only the shear wave velocity is perturbed. The damping factor is chosen from the misfit and model smoothness tradeoff 13 Figure 2.4: Ambient noise phase and group velocity maps (left), checkerboard resolution maps (center), and model error maps (right) at 34 s period. The volcanic provinces and the Veracruz basin are delimited by the solid and dashed contour lines, respectively. Note how the velocity distribution varies for all four velocity types due to their difference in sensitivity. curve. We also tested the case of reducing the damping factor with depth to account for the large heterogeneity of the uppermost crust but found no significant differences in the results. The final model is obtained by iteratively perturbing the initial model until a good fit to the two dispersion curves is achieved. Finally, we use regional average structures (1◦ x 1◦ ) as initial models and perform one last inversion on a finer grid. For this step, each part of the dispersion curves is weighted inversely proportional to the model error in the tomographic inversions. The final 1-D velocity profiles are then combined to form two orthogonally polarized three-dimensional VS models of the whole study region. In a simple isotropic medium, both the VSV and VSH models would be identical, as shear waves travel at the same speed regardless of their polarization. However, in the presence of variable tectonic stresses, complicated structures 14 Figure 2.5: Cross-sections of the Rayleigh wave group (A) and phase (B) velocity models and the Love wave group (C) and phase (D) velocity models along the MASE line (F-F’ in Figure 2.1). Topography is shown above each profile. begin to emerge and reorient so that the shear wave velocity of the medium also depends on direction of propagation and polarization. The velocity dependence of the latter is termed radial anisotropy and is estimated here as: γ= VSH − VSV , VS (2.1) where VS is the isotropic or effective shear wave velocity and is directly computed from the VSV and VSH models via a Voigt’s average: s VS = 2 +V2 2VSV SH . 3 (2.2) The velocity distribution derived from Equation (2.2) is then use to construct our final velocity model since it better reflects the apparent variations in elastic properties [44, 45]. Figure 2.6 shows vertical cross-sections of our derived Voigt-averaged VS and radial anisotropy model along the MASE and VEOX seismic lines (F-F’ and G-G’ in Figure 2.1, respectively) with their associated average misfits. An obvious feature in the VS profiles is that superficial slow velocities correspond well with the TMVB and even more so for the LTVF, where the recently active San Martín Tuxtla Volcano is located. The radial anisotropy distribution, on the other hand, appears to be more contrasting and related to the subduction geometry. The main mechanisms responsible for 15 causing radial anisotropy in a subduction environment, and how its presence can be interpreted, are addressed in the "Results and Discussion" section. Figure 2.6: Vertical cross-sections of the Voigt-averaged VS (A and B) and radial anisotropy models (C and D) along the MASE and VEOX lines (FF’ and G-G’ in Figure 2.1, respectively). Topography is shown above each profile. The hypocenters from the SSN catalog within ±10 km are projected to each cross-section. The slab models from Pérez-Campos et al. [13] and Melgar and Pérez-Campos [17] are plotted in black lines assuming a constant plate thickness of 40 km. The thin black dashed line in profiles A, B, and D represent the Moho interface, defined by where the shear wave velocity first exceeds 4.1 km s−1 , and the thick green and red dashed lines in profiles C and D represent regions of elevated positive (≥VSH ) and negative radial anisotropy (>VSV ), respectively. Cross-sections E and F show the average misfit of the 1-D inversion along the same profiles. Inversion for Azimuthal Anisotropy To characterize the azimuthal anisotropy, or the dependence of wave-speed with azimuth of propagation, we adopt a beamforming approach and fit the first three parameters in Smith and Dahlen [46] anisotropy model for Rayleigh waves: 16 υ(T, θ) =a0 (T ) + a1 (T ) cos (2θ) + a2 (T ) sin(2θ)+ a3 (T ) cos (4θ) + a4 (T ) sin(4θ), (2.3) where υ is the surface wave phase velocity, T the period, θ the backazimuth, a0 the isotropic velocity, and a1−4 the azimuthal coefficients [47], to the ZZ crosscorrelations. Only the isotropic and 2θ coefficients are considered here because of the very small contribution that came from the rest of the parameters in our initial inversion. After characterizing the wavefield’s azimuthal dependence at every station, we calculate the amplitude of the anisotropy, A, and its seismically fast direction, φ, using: A= φ= q a1 2 + a2 2 , (2.4) a2 1 arctan . 2 a1 (2.5) In the traditional beamforming method, one inverts the phase information by finding the best fitting slowness and backazimuth of a plane wave, thus providing a detailed characterization of the seismic wavefield at a given location [e.g., 48]. Here, we isolate stations one at a time, and use the remainder as virtual sources to find the average phase velocity of Rayleigh waves traveling to the reference station from all available azimuths. To ensure the robustness of our measurements, we only beamform cross-correlations with a broadband SNR higher than 10 and an inter-station distance larger than one wavelength of the lowest period of the bandpass filters. We also assume that the wavefield’s full azimuthal dependance can only be characterized if the azimuth range of 180◦ is sampled by at least 3 paths in a 5-bin range [39]. In the actual beamforming process, we search for the maximum coherent output over velocities from 1-5 km/s and every 5◦ from 0-360◦ backazimuth with 70% overlap for the 3-20 and 20-50 s period bands. Such period bands are determined empirically based on the Rayleigh wave phase velocity sensitivity to perturbations in VS in an effort to characterize the upper crust anisotropy and the lower crust and upper mantle anisotropy independently (Figure 2.7a). The sensitivity kernels are computed using the modified Tectonic North America Model (mTNA) [11]. A visual inspection of Figure 2.7bc clearly shows that the surface wave phase velocity varies with backazimuth differently over the two frequency bands. To 17 find the best fitting coefficients to these variations, we perform a grid-search using a weighted L2 regularization in which each azimuth bin is weighted proportionally to the amount of energy that was beamformed. The resulting best fits for two different frequency bands are shown in Figure 2.7bc as dashed red lines with the best fitting parameters given at the bottom right. A possible source of error in our procedure may come from the inherent tradeoff between seismic heterogeneity and azimuthal anisotropy. However, given the broad frequency bands we make our measurements in, it is reasonable to assume that any effect caused by lateral heterogeneities will be minimal when compared to the one introduced by large-scale tectonic processes. To validate this assumption, we compared a stacked map of all the de-trended Rayleigh wave phase velocity maps for the 20-50 s period band with the results of our beamforming process in the same period band for a given station in central Mexico. The comparison showed clear systematic velocity variation as a function of backazimuth with no direct correlation with the isotropic velocity distribution. Lastly, to assess the uncertainty in our fit parameters as well as their statistical significance, we estimate the 95% bootstrap confidence limits using 100 resamples. 2.4 Results and Discussion Shear Wave Velocity and Radial Anisotropy The shear wave velocity structure in central and southern Mexico has been discussed in a number of seismic imaging studies [e.g., 40, 49, 50]. However, the enhanced ray coverage that is now provided by the GECO, OXNET, Veracruz, and expanded SSN networks allows us to resolve new features of the continental crust and uppermost mantle that were averaged out due to insufficient resolution, particularly at the eastern sector of the TMVB. Figures 2.8 and 2.9 shows the results of our three-dimensional VS model near the coastal Veracruz basin. Among the most striking features observed in the profiles of Figure 2.8 is the dramatic variation in crustal thickness that occurs in the transition of the TMVB and the LTVF, and how the crust thins out towards the Gulf of Mexico. Taking a 4.1 km/s limit in VS as the crust and mantle transition zone [e.g., 51], our results indicate that the Moho reaches a depth of about 45 km beneath the eastern TMVB and rises sharply to 20 km beneath the Veracruz basin just before lowering once more to an average depth of 40 km beneath the LTVF. Interestingly enough, this abrupt change in the Moho 18 Figure 2.7: Sensitivity kernels for the fundamental mode Rayleigh-wave phase velocity calculated from the modified Tectonic North America Model (A) [11]. Two-dimensional histograms over the backazimuth-phase velocity space for the 3–20 seconds period band (B) and the 20-50 seconds period band (C). The orange bars on top of each panel represent the number of cross-correlations beamformed in each backazimuth bin. The dashed red line indicates the best fit of Equation 2.3 and the gray diamonds mark the beam velocity in each bin (with size as a function of beam power). The best fitting parameters are given at the bottom right of each panel. depth occurs just above where Dougherty and Clayton [10] propose the existence of a plate tear, and it encompasses the same ∼50-75-km-wide zone where they observe a sharp decrease in intraslab seismicity (Figure 2.9). The crustal thickness measurements inferred from this study are, in general, consistent with previous experiments in the area [19–21]. However, using a fixed velocity contour to define the Moho can be problematic if other physical processes that are altering the crust’s chemistry are present. In the shallow crust, radial anisotropy is thought to result from shape-preferred orientation (SPO) in the form of fine layering, metamorphic foliations, and multi-scale fractures in the crystalline rocks exposed to regional stresses [e.g., 52, 53]. Crystallographic-preferred orientation (CPO) of intrinsically anisotropic minerals can also cause strong anisotropy at shallow depths, but this mechanism seems to only be dominant in the upper mantle, where olivine and 19 Figure 2.8: Vertical cross-sections of the Voigt-averaged VS model along the A-A’ and H-H’ lines in Figure 2.1. Both profiles show extreme crustal thinning in the Veracruz basin. The dashed red line in panel (A) indicates the location where Dougherty and Clayton [10] propose the existence of a plate tear. Topography is shown above each profile. Other symbols are as in Figure 2.6ab. Figure 2.9: Vertical cross-sections of the Voigt-averaged VS model in the Veracruz basin (from the A-A’ to the E-E’ lines in Figure 2.1). The color scale is saturated to only show the velocity anomalies present in the upper mantle. The solid black line indicates the Moho interface. The thick red line indicates the location where Dougherty and Clayton [10] propose the existence of a plate tear, and the dashed red line delimits the 75-km-wide zone where they observe a sharp change in intraslab seismicity. Cross-sections are separated by a constant offset of 20 km. The eye marker in Figure 2.1 indicates the viewpoint of the observer. orthopyroxene aggregates are deformed during oriented geodynamic processes [54]. Anisotropy in the lower part of the continental crust generally does not have a strong seismic signature [55]. However, in subduction zone environments, significant radial anisotropy can be found in the overall crust 20 and uppermost mantle arising from melt-filled cracks, lenses of partial melts, and other large-scale intrusive bodies such as dykes and sill complexes [e.g., 37, 56]. The analysis of the strength and spatial extent of radial anisotropy can therefore be used to identify heterogeneous compositions in a wide range of depths. However, due to the potential biases in the background VSV and VSH velocities that arise from the difference in the regularizations imposed in the tomographic inversions, we limit the interpretation of our results to the changes in the sign of the radial anisotropy rather than the amplitude. An additional source of error in our radial anisotropy measurements may come from the leakage of Rayleigh wave energy in the TT components, which can be particularly strong at long periods and large interstation distances [57]. As previously stated, Figure 2.6C shows a complicated radial anisotropy pattern along the MASE line. A clear plume-like structure with negative radial anisotropy is present in the mantle wedge, where the VS model predicts low velocities and the subducted plate starts diving into the continental mantle at a steep angle [13]. This feature is most likely related to the presence of partial melts and fluid upwelling since SV waves are traveling faster than SH waves in this zone. The ascent of this material appears to be confined to the base of the continental crust, at the point in our model where there is a sharp transition from negative to positive radial anisotropy. Considering that SV waves are being slowed in this part of the crust, we believe that this anisotropy contrast is due to the existence of molten material piled horizontally such as a large sill complex fed by the partial melts just beneath it. It is also worth noting a secondary positive radial anisotropic anomaly rooting-up from the flat portion of the Cocos plate. Although the cause of this anomaly is unclear, its location is coincident with where [58] revealed a low resistivity zone that is related to the dehydration of the subducted slab. The slab itself, on the other hand, has a predominantly negative radial anisotropy that might reflect vertically oriented faults or fluid-filled cracks that are functioning as the primary conduit for slab dehydration. Interestingly enough, the zones with the highest negative anisotropy within the slab are close to the two patches in which tremors appear to occur regularly (the so-called sweet spot and transient zone [59, 60]). This last observation is in agreement with previous receiver function results that revealed elevated shear wave splitting in these portions of the subducted oceanic crust along the MASE line [61]. 21 Our radial anisotropy model along the VEOX line shows a marked negative anomaly beneath the LTVF (Figure 2.6D). This feature is commonly observed in volcanic environments and is probably related to the presence of a set of vertically oriented conduits or dikes that is storing the magma beneath the stratovolcanoes, such as the one Spica et al. [37] imaged beneath the Colima volcano in Jalisco, southwestern Mexico. The source of the fluids and melts feeding the magma storage are unknown since the projected Cocos slab is far too deep at this point to have an immediate effect on the overriding plate. We also do not find any evidence in the surface wave images of the anomalous south-dipping structure seen in the Kim et al. [18] receiver function image. However, the isolated large positive radial anisotropy present in the 30-40 km depth range beneath the LTVF indicates the possibility that melts are migrating laterally rather than being generated directly beneath the volcanoes. This observation is more intriguing in view of the fact that a magnetotelluric profile collected at approximately 100 km west of the LTVF revealed an unresolved low resistivity zone beneath the Veracruz basin [58]. The radial anisotropy distribution along the VEOX line is a key element for establishing a relationship between the transition structure of the Cocos slab and the Central America volcanism and is discussed with more detail in the "The Eastern End of the TMVB" subsection of this chapter. Figure 2.10 shows the radial anisotropy along the same profiles as the ones presented in Figure 2.9. Although the anisotropy distribution appears to be smoother than the one observed along the MASE and VEOX experiments (where slab material is present at much shallower depths), there is a sharp positive radial anisotropy anomaly that delimits most of the crust-mantle transition zone. This observation is not too surprising since anisotropy in the uppermost mantle is generally governed by the systematic flow alignment of olivine crystals under dislocation creep, thus allowing SH waves to travel faster on average than SV waves [54, 62]. What is intriguing however is a subtle discontinuity in the anisotropic pattern at again, just above where Dougherty and Clayton [10] propose the existence of a plate tear. This interruption in what appears to be the horizontal creeping flow of the mantle is indicative of a sharp structural change and may be related to some vertical flow component [e.g., 63]. Note that the lack of stations in and near the Gulf of Mexico may limit the resolution of these images. 22 Figure 2.10: Vertical cross-sections of the radial anisotropy model in the Veracruz basin (same profiles as in Figure 2.9). Azimuthal Anisotropy The slab roll-back that has built the TMVB has proceeded a distance of 150 km over the past 20 Myr [64]. This process requires a substantial movement of mantle material from the back to the front of the slab that can only be achieved by mantle flow beneath or around the subducted slab. Regardless of its path and orientation, this ductile flow is likely to produce a strong olivine CPO and lead to a bulk seismic anisotropy that should be observable at various scales. Under dry mantle conditions, the seismically fast olivine a-axis generally aligns with the shear direction [65, 66]. However, recent experimental work has shown that the presence of water in the medium can change the olivine a-axis orientations to perpendicular to the mantle flow direction [67]. This configuration is referred to as olivine type-B, whereas the relationship of dry olivine is type-A. In typical subduction zones, the mantle wedge tip meets the conditions for the existence of type-B olivine whereas type-A olivine is found throughout the mantle wedge core [68]. However, due to the young age (∼14 Ma) of the Cocos plate and its high temperature (>900◦ C), we can expect that most of the azimuthal anisotropy in the mantle wedge in central and southern Mexico is dominated by type-A olivine CPO [3, 14, 61, 69, 70]. Consistent with a previous anisotropy study of the region [11], our measure- 23 ments show a preferential trench-parallel alignment of the seismically fast direction beneath the subducted slab in most of the forearc (Figure 2.11). However, our anisotropy map displays, for the first time, a well-defined flow pattern that supports the hypothesis of a slab tear in South Cocos. In the case of a slab tear or gap, mantle material could flow through the slab window and create the mantle fabric depicted by our anisotropic parameters. The variations in slab dip and rollback rate of the plate along the MAT would also induce a 3-D flow field, in which mantle materials move from beneath the slab to the mantle wedge underneath Veracruz. A similar type of flow in the Rivera-Cocos plate boundary has been proposed to control the asthenospheric anisotropy that is observed beneath the MARS and NARS seismic arrays [71, 72], which would explain the subtle north-northeast rotation in our fast axes in the Rivera segment of the MAT. Our results cannot provide definite evidence for the slab tear that has been proposed to exist along the projection of the OFZ in western Mexico [9] as the seismic coverage around this area is sparse. Anisotropy in the back-arc of central Mexico (i.e. the northern section of the MASE array), on the other hand, appears to be primarily controlled by a trench-perpendicular 2-D corner flow that is induced by the abrupt down-dip motion of the slab [70]. The Eastern End of the TMVB There has been significant progress in understanding the TMVB and the Central American subduction system. However, the nature of the pronounced change in arc volcanism in central and southern Mexico has remained ambiguous. On the basis of structural arguments, our study indicates that a tear in the South Cocos slab can explain most of the enigmatic features that characterize this segment of the MAT (Figure 2.12). For the remainder of this chapter, we will refer to the northern part of South Cocos as Central Cocos whereas the slab segment that is south of the possible tear will remain as South Cocos. The first piece of evidence for a possible slab tear separating Central and South Cocos comes from geomorphological features. The general extent of the NNE volcanic chain at the easternmost TMVB suggests that this composite structure is linked to a source of magma and fluids that reached the surface in a linearly distributed order. This narrowly localized source of melts may have been associated with an upwelling of the isotherms due to the asthenospheric 24 Figure 2.11: Map of the azimuthal anisotropy results for the lower crust and upper mantle. The orientation of the vectors gives the seismically fast direction, φ, and the length of the lines is proportional to the amplitude of the anisotropy, A. The thick green arrows depict the possible trajectories of mantle flows inferred from our measurements. The red oval delimits the eastern TMVB volcanic chain. mantle material flowing around the edge of the slab during tear propagation, a mechanism similar to the one that formed the cross-back-arc volcanic trail in the Ryukyu subduction zone in Japan [73]. There is also a systematic southward progression in age along the volcanic chain [8], which is characteristic of plate rollback and/or tear development [74]. However, Siebert and CarrascoNúñez [75] report the presence of young basaltic rocks in the northern part of the chain. The chronological similarity between the activity at the current volcanic front and in the back-arc region suggests that volcanism is not only controlled by regular subduction and rollback but also by a mechanism, such as an elongated slab window, that could allow melts to reach the back-arc at a relatively contemporary time. Moreover, Gómez-Tuena et al. [76] suggest that the dramatic change in the composition of volcanic rocks in the easternmost TMVB is associated with a gradual increase in the angle of subduction at the end of the Miocene, which would allow the partial melting of a relatively deeper mantle source. Nonetheless, the flow of mantle materials through a slab tear, with the accompanying slab-edge melting, may also explain the acute change in the source of melts and the adakitic signature observed in young rocks in 25 the eastern sector of the TMVB [e.g., 77–79]. We suggest that the building of the volcanic chain in the eastern TMVB represents an early stage of the slab tear development. The abrupt termination of the TMVB at its eastern end, and the resulting discontinuity of arc volcanism, introduces complexity into the slab tear hypothesis but does not abate it. To explain the absence of surficial volcanism directly in the coastal Veracruz basin, Dougherty and Clayton [10] suggest that the tear is a less developed or young feature located in the downdip portion of the slab, and that consequently any material flowing through it is not rising to a sufficiently shallow depth to have an effect on the overriding plate. We propose a slightly different scenario. Rather than the tear being underdeveloped, we hypothesize that the rapid rollback rate of South Cocos relative to Central Cocos, accelerated by the influx of less dense asthenosphere material into the mantle wedge through the slab window [71, 80], introduces a suction force that increases the strength of the toroidal flow through the tear and drags the mantle materials under Central Cocos laterally toward the south of Mexico. As a result, any material flowing through the tear is unable to reach the surface directly at the Veracruz basin. The main evidence for this scenario comes from the azimuthal anisotropy results (Figure 2.11), where the asthenospheric mantle appears to be flowing almost trench-parallel beneath the slab until it flows out into the mantle wedge through the possible tear and heads towards southeastern Mexico. The subtle interruption in the radial anisotropy pattern at the eastern end of the TMVB (Figure 2.10) also suggests the dominance of a horizontal mantle flow in the presence of a weak and localized vertical flow just where the tear is proposed to be. Furthermore, we suggest that some of the active volcanoes of the eastern TMVB volcanic chain are still being fed by some slab-edge melting mechanism, and that the southeastern drag of mantle materials below the Veracruz area is responsible for the isolated volcanic expressions that are present in central and southern Mexico. As described earlier, there is no deep source of melts feeding the LTVF in the same way that we observe for the TMVB in central Mexico (Figure 2.6). Instead, there is a strong positive radial anisotropy, indicative of accumulation of horizontal or sub-horizontal lava flows, that appears to be transporting the magma that is building the volcanoes at the LTVF. We suggest that the shear stress exerted by the toroidal flow around the slab tear is responsible for orienting the magma flow horizontally and that the melts are then able to reach the surface 26 through a series of vertically oriented conduits that are just below the LTVF. This scenario is similar to the process that created Mount Etna in Europe [81] and would explain the younger age of LTVF magmas (7 Ma) relative to the overall age of the TMVB [82]. The horizontal transport of mantle material in the uppermost mantle and lower crust underneath the Veracruz basin is further supported by geochemical data that advocates that the LTVF source is likely to reside in the lithosphere rather than the asthenosphere, and would also explain why the LTVF evolved magmas do not require a slab component in their genesis [83]. A similar scenario, in which the anomalous setting of the southern Mexico volcanism is related to a mantle flow readjustment, was proposed by Manea et al. [84]. However, different from this study, Manea et al. [84] argues that materials in the mantle wedge are being pushed laterally out of central Mexico due to the ongoing slab flattening process, and that the combination of this flow with a decrease in temperature around the flat slab area may be responsible for the discontinuity of arc volcanism. Even though the slab flattening in central Mexico can also provide a plausible explanation for the spatial variation of the volcanic arcs, the idea of a slab tear and a perturbed mantle wedge underneath the Veracruz basin is more appealing for two main reasons. First, the influx of hot, less-dense asthenospheric mantle material flowing through the tear to the Veracruz area may explain how the contrasting topographic relief between the eastern TMVB and coastal Veracruz basin is supported by the variations in the crustal thicknesses. Second, the temperature change introduced by the hot mantle flow may significantly influence the thermomechanical state of the subducting slab and explain the abrupt termination of the seismicity associated with the Cocos plate at the southern end of the Veracruz basin. Among the many hypotheses put forward to explain the occurrence of intermediate-depth earthquakes, dehydration embrittlement [85], which is the brittle failure associated with dehydration reactions of hydrous minerals in the slab and upper mantle, is considered to be the leading mechanism [86]. This hypothesis asserts that intermediate-depth earthquakes occur in subducting slabs where dehydration is expected but are absent from parts of the slab predicted to be anhydrous [87]. To this day, a direct link between lithospheric tearing and intermediate-depth seismicity has not been established. Meighan et al. [88] report a positive correlation between slab tears and intermediate-depth seismicity and suggest that most of the seismic sequences that occur in these environments are associated with the 27 Figure 2.12: 3-D schematic illustration of the proposed tectonic setting and inferred mantle flow (green arrows). Asthenospheric mantle materials are flowing almost trench-parallel beneath the slab in the forearc region until they are redirected to the mantle wedge through a tear that separates Central and South Cocos. The accelerated rollback rate of South Cocos relative to Central Cocos introduces a suction force that further displaces the mantle materials laterally toward the south of Mexico. The shear stress exerted by the toroidal flow around the slab tear may be transporting the melts and fluids that are feeding the volcanoes in the LTVF from central Mexico, thus allowing their existence without any slab material at an appropriate depth directly underneath them. Some of the active volcanoes composing the eastern TMVB volcanic chain however might still be fed by some slab-edge melting mechanism. Anisotropy in the back-arc of central Mexico (i.e. Central Cocos) is primarily controlled by a trench-perpendicular 2-D corner flow that is induced by the down-dip motion of the slab. The flow of mantle materials through a slab tear, with the accompanying slab-edge melting, may explain the exceptional spatial distribution of stratovolcanoes in the eastern sector of the TMVB, as well as the abrupt change in the source of melts observed in young rocks in this segment of the MAT. The eye in Figure 2.1 marks the viewpoint of the observer. fluid-related embrittlement of mantle rocks. However, we argue that, in the case of the subduction system in Mexico, the elevated rate of shear heating due to strong lateral mantle materials flowing through the tear has deprived northern South Cocos of most of the fluids stored in its crust. This, in turn, prevents intra-plate earthquakes from occurring due to the absence of dehy- 28 dration embrittlement at the northern South Cocos’ intermediate depth. An alternate scenario that may explain the absence of subduction-related seismicity in the Veracruz region is slab-detachment, which involves the detachment of a portion of the slab during ongoing subduction [89]. However, the lack of a strong topographic response directly in the Veracruz basin (i.e. uplift) or any sharp lateral contrasts in the trench topography [90] suggests that an abrasive process such as slab break-off might not be responsible for the absence for seismicity and surface volcanism in this area. The existence of a horizontal separation of the two slab segments with a relatively large gap in between due to slab bending [91] may also explain the absence of intermediate-depth seismicity underneath Veracruz. Yet, focal mechanisms of regional intraslab earthquakes do not show a trench-parallel trending T-axis orientation along the sharp transition in slab dip [10], suggesting that the tear, if existent, is more likely to be propagating vertically rather than horizontally. Note that both the down-dropped side of the slab (in the case of a vertically propagating tear) and a large horizontal gap between the slab segments could have created the topographic hole for the Veracruz basin, but only the vertically propagating tear scenario would explain the existence of the LTVF without a deep source of melts directly underneath the stratovolcanoes. Slab tearing has been observed in numerous subduction systems around the world and is generally believed to be triggered by local collisional events [92] or by the variations in age, temperature, geometry, and convergence rate of the subducting plate along the trench [e.g., 93, and references therein]. For the case of southern Mexico, Dougherty and Clayton [10] suggest that the possible tear separating South and Central Cocos may be related to the subduction of several parallel ridges or topographic heterogeneities off the coast of Oaxaca, or, alternatively, to the accommodation of strain due to the sharp change in slab geometry (or some combination of the two). Together with the slab window along the Rivera-Cocos plate boundary, and the tear in western Mexico, a tear in southern Mexico has important implications in the subduction dynamics of Central America. Aside from exerting strong influence on the surface volcanism, the composite mantle flow formed by the movement of asthenospheric mantle materials from the back to the front of the slab through the gap and tears may be allowing the separate segments of the Cocos plate to rollback independently and further promote the large variations in dip that characterize the Mexico subduction system Dougherty and Clayton [10]. Our 29 proposed scenario wherein the hot, less-dense asthenospheric mantle material is flowing through the tear to the Veracruz area is consistent with regional thermal models that suggest the existence of an upper mantle with relatively low density and high temperatures in southern Mexico [94]. However, further heat flow measurements as well as geodynamical modeling of the mantle flow in this setting are required in order to provide more quantitative insights into the role of the slab tear in the Mexico subduction system. 2.5 Conclusions We have used the discrepancy between Rayleigh and Love waves derived from ambient seismic noise cross-correlations to develop high-resolution anisotropic images for central Mexico and showed that the transition from flat to steeper subduction is more likely to be accommodated by a slab tear than a sharp flexure. A 2-D profile of the radial anisotropy model along the MASE line shows apparent melt migration pathways that go from the steeply dipping portion of the subducted slab to a sill complex in the lower crust beneath the TMVB. Moreover, the flat portion of the slab appears to be dominated by vertically oriented structures, such as fluid-filled cracks, that may be functioning as the primary conduit for slab dehydration. Along the VEOX line, our images suggest that volcanoes at the LTVF are fed by lava flows that are sub-horizontally transported in the lower crust. Our azimuthal anisotropy results indicate that mantle material flows trench-parallel underneath the slab until it is redirected to the mantle wedge through a possible tear that separates Central and South Cocos. The difference in the rollback rates of the plates introduces a suction force that further displaces the asthenospheric mantle material towards southern Mexico. We hypothesize that the shear stress exerted by this toroidal flow is transporting the lava flows from the slab window to the LTVF, and may also be responsible for the few other isolated volcanoes present in southern Mexico. The lack of subduction-related earthquakes underneath the Veracruz basin may therefore be explained by two plausible scenarios: A) The elevated rate of shear heating due to the rapid influx of hot, less-dense material flowing through the tear underneath Veracruz has depleted most of the fluids that were once stored in the slab’s oceanic crust. As a result, intra-plate earthquakes do not occur due to the absence of dehydration embrittlement at intermediate depths. B) The tear is accommodated by a lateral movement promoted by slab bending such that there exists a relatively large gap in between the two 30 slab segments underneath the Veracruz basin. Given that a large lateral slab gap beneath Veracruz would not explain the existence of the LTVF without a deep source of melts below the stratovolcanoes, and that there appears to be no obvious topographic response to this process, we favor the first scenario over the second one. References [1] Alain Demant. Características del eje neovolcánico transmexicano y sus problemas de interpretación. Revista mexicana de ciencias geológicas, 2 (2):172–187, 1978. [2] Arturo Gómez-Tuena, Ma Teresa Orozco-Esquivel, and Luca Ferrari. Petrogénesis ígnea de la faja volcánica transmexicana. Boletín de la Sociedad Geológica Mexicana, 57(3):227–283, 2006. [3] Mario Pardo and Gerardo Suarez. Shape of the subducted rivera and cocos plates in southern mexico: Seismic and tectonic implications. Journal of Geophysical Research: Solid Earth, 100(B7):12357–12373, 1995. [4] Gerardo Suárez and Osvaldo Sánchez. Shallow depth of seismogenic coupling in southern mexico: implications for the maximum size of earthquakes in the subduction zone. Physics of the earth and planetary interiors, 93(1-2):53–61, 1996. [5] Stephan Husen, R Quintero, E Kissling, and B Hacker. Subduction-zone structure and magmatic processes beneath costa rica constrained by local earthquake tomography and petrological modelling. Geophysical Journal International, 155(1):11–32, 2003. [6] Ivonne Gabriela Arroyo, Stephan Husen, Ernst R Flueh, J Gossler, E Kissling, and GE Alvarado. Three-dimensional p-wave velocity structure on the shallow part of the central costa rican pacific margin from local earthquake tomography using off-and onshore networks. Geophysical Journal International, 179(2):827–849, 2009. [7] Ellen M Syracuse, Geoffrey A Abers, Karen Fischer, Laura MacKenzie, Catherine Rychert, Marino Protti, Víctor González, and Wilfried Strauch. Seismic tomography and earthquake locations in the nicaraguan and costa rican upper mantle. Geochemistry, Geophysics, Geosystems, 9(7), 2008. [8] Luca Ferrari, Teresa Orozco-Esquivel, Vlad Manea, and Marina Manea. The dynamic history of the trans-mexican volcanic belt and the mexico subduction zone. Tectonophysics, 522:122–149, 2012. 31 [9] Sara L Dougherty, Robert W Clayton, and Don V Helmberger. Seismic structure in central mexico: Implications for fragmentation of the subducted cocos plate. Journal of Geophysical Research: Solid Earth, 117 (B9), 2012. [10] Sara L Dougherty and Robert W Clayton. Seismicity and structure in central mexico: Evidence for a possible slab tear in the south cocos plate. Journal of Geophysical Research: Solid Earth, 119(4):3424–3447, 2014. [11] Igor Stubailo, Caroline Beghein, and PM Davis. Structure and anisotropy of the mexico subduction zone based on rayleigh-wave analysis and implications for the geometry of the trans-mexican volcanic belt. Journal of Geophysical Research: Solid Earth, 117(B5), 2012. [12] WL Bandy, TWC Hilde, and CY Yan. The rivera-cocos plate boundary: Implications for rivera-cocos relative motion and plate fragmentation. Geological Society of America Special Papers, 334:1–28, 2000. [13] Xyoli Pérez-Campos, YoungHee Kim, Allen Husker, Paul M Davis, Robert W Clayton, Arturo Iglesias, Javier F Pacheco, Shri K Singh, Vlad Constantin Manea, and Michael Gurnis. Horizontal subduction and truncation of the cocos plate beneath central mexico. Geophysical Research Letters, 35(18), 2008. [14] Allen Husker and Paul M Davis. Tomography and thermal state of the cocos plate subduction beneath mexico city. Journal of Geophysical Research: Solid Earth, 114(B4), 2009. [15] Y Kim, RW Clayton, and JM Jackson. Geometry and seismic properties of the subducting cocos plate in central mexico. Journal of Geophysical Research: Solid Earth, 115(B6), 2010. [16] Luca Ferrari. Slab detachment control on mafic volcanic pulse and mantle heterogeneity in central mexico. Geology, 32(1):77–80, 2004. [17] Diego Melgar and Xyoli Pérez-Campos. Imaging the moho and subducted oceanic crust at the isthmus of tehuantepec, mexico, from receiver functions. Pure and applied geophysics, 168(8-9):1449–1460, 2011. [18] YoungHee Kim, Robert W Clayton, and Fraser Keppie. Evidence of a collision between the yucatán block and mexico in the miocene. Geophysical Journal International, 187(2):989–1000, 2011. [19] R Molina-Garza and J Urrutia-Fucugauchi. Deep crustal structure of central mexico derived from interpretation of bouguer gravity anomaly data. Journal of Geodynamics, 17(4):181–201, 1993. 32 [20] Jaime Urrutia-Fucugauchi and Jesus Hernan Flores-Ruiz. Bouguer gravity anomalies and regional crustal structure in central mexico. International Geology Review, 38(2):176–194, 1996. [21] Victor Hugo Espíndola, Luis Quintanar, and Juan Manuel Espíndola. Crustal structure beneath mexico from receiver functions. Bulletin of the Seismological Society of America, 107(5):2427–2442, 2017. [22] Luca Ferrari, Takahiro Tagami, Mugihiko Eguchi, Ma Teresa OrozcoEsquivel, Chiara M Petrone, Jorge Jacobo-Albarrán, and Margarita López-Martínez. Geology, geochronology and tectonic setting of late cenozoic volcanism along the southwestern gulf of mexico: The eastern alkaline province revisited. Journal of Volcanology and Geothermal Research, 146 (4):284–306, 2005. [23] Stephen A Nelson, Erika Gonzalez-Caver, and T Kurtis Kyser. Constraints on the origin of alkaline and calc-alkaline magmas from the tuxtla volcanic field, veracruz, mexico. Contributions to Mineralogy and Petrology, 122(1-2):191–211, 1995. [24] Shannon Fasola, Michael R Brudzinski, Noorulann Ghouse, Katharine Solada, Stefany Sit, Enrique Cabral-Cano, Alejandra Arciniega-Ceballos, Nicholas Kelly, and Kevin Jensen. New perspective on the transition from flat to steeper subduction in oaxaca, mexico, based on seismicity, nonvolcanic tremor, and slow slip. Journal of Geophysical Research: Solid Earth, 121(3):1835–1848, 2016. [25] Xyoli Pérez-Campos, Víctor Hugo Espíndola, Jesús Pérez, Jorge A Estrada, Caridad Cárdenas Monroy, Delia Bello, Adriana GonzálezLópez, Daniel González Ávila, Moisés Gerardo Contreras Ruiz Esparza, Rafael Maldonado, et al. The mexican national seismological service: An overview. Seismological Research Letters, 2018. [26] Francisco Córdoba-Montiel, Arturo Iglesias, Xyoli Pérez-Campos, K Sieron, Grupo de Trabajo del Servicio Sismológico Nacional, C ValdésGonzález, SK Singh, and JF Pacheco. The broadband seismological network of veracruz, mexico: Toward a regional seismotectonic interpretation. Seismological Research Letters, 89(2A):345–355, 2018. [27] Oleg I Lobkis and Richard L Weaver. On the emergence of the green’s function in the correlations of a diffuse field. The Journal of the Acoustical Society of America, 110(6):3011–3017, 2001. [28] Michel Campillo and Anne Paul. Long-range correlations in the diffuse seismic coda. Science, 299(5606):547–549, 2003. [29] Nikolai M Shapiro and Michel Campillo. Emergence of broadband rayleigh waves from correlations of the ambient seismic noise. Geophysical Research Letters, 31(7), 2004. 33 [30] Roel Snieder. Extracting the green’s function from the correlation of coda waves: A derivation based on stationary phase. Physical Review E, 69(4): 046610, 2004. [31] Korbinian Sager, Christian Boehm, Laura Ermert, Lion Krischer, and Andreas Fichtner. Sensitivity of seismic noise correlation functions to global noise sources. Journal of Geophysical Research: Solid Earth, 0(ja), 2018. doi: 10.1029/2018JB016042. URL https://agupubs.onlinelibrary. wiley.com/doi/abs/10.1029/2018JB016042. [32] GD Bensen, MH Ritzwoller, MP Barmin, AL Levshin, F Lin, MP Moschetti, NM Shapiro, and Yanyan Yang. Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements. Geophysical Journal International, 169(3):1239–1260, 2007. [33] Jack B Muir and Victor C Tsai. Rayleigh-wave h/v via noise cross correlation in southern california. Bulletin of the Seismological Society of America, 107(5):2021–2027, 2017. [34] Thanh-Son Pha.m, Hrvoje Tkalčić, Malcolm Sambridge, and Brian LN Kennett. Earth’s correlation wavefield: late coda correlation. Geophysical Research Letters, 45(7):3035–3042, 2018. [35] Huajian Yao, Robert D van Der Hilst, and Maarten V De Hoop. Surfacewave array tomography in se tibet from ambient seismic noise and twostation analysis—i. phase velocity maps. Geophysical Journal International, 166(2):732–744, 2006. [36] Yinhe Luo, Yingjie Yang, Yixian Xu, Hongrui Xu, Kaifeng Zhao, and Kai Wang. On the limitations of interstation distances in ambient noise tomography. Geophysical Journal International, 201(2):652–661, 2015. [37] Zack Spica, Mathieu Perton, and Denis Legrand. Anatomy of the colima volcano magmatic system, mexico. Earth and Planetary Science Letters, 459:1–13, 2017. [38] MP Barmin, MH Ritzwoller, and AL Levshin. A fast and reliable method for surface wave tomography. In Monitoring the Comprehensive NuclearTest-Ban Treaty: Surface Waves, pages 1351–1375. Springer, 2001. [39] Eric Debayle and Malcolm Sambridge. Inversion of massive surface wave data sets: model construction and resolution assessment. Journal of Geophysical Research: Solid Earth, 109(B2), 2004. [40] Zack Spica, Mathieu Perton, Marco Calò, Denis Legrand, Francisco Córdoba-Montiel, and Arturo Iglesias. 3-d shear wave velocity model of mexico and south us: bridging seismic networks with ambient noise crosscorrelations (c1) and correlation of coda of correlations (c3). Geophysical Journal International, 206(3):1795–1813, 2016. 34 [41] Yiran Ma and Robert W Clayton. The crust and uppermost mantle structure of southern peru from ambient noise and earthquake surface wave analysis. Earth and Planetary Science Letters, 395:61–70, 2014. [42] Charles J Ammon, Minoo Kosarian, Robert B Herrmann, Michael E Pasyanos, William R Walter, and Hrvoje Tkalcic. Simultaneous inversion of receiver functions, multi-mode dispersion, and travel-time tomography for lithospheric structure beneath the middle east and north africa. In Proceedings of the 26th Seismic Research Review—Trends in Nuclear Explosion Monitoring, pages 17–28, 2004. [43] GHF Gardner, LW Gardner, and AR Gregory. Formation velocity and density—the diagnostic basics for stratigraphic traps. Geophysics, 39(6): 770–780, 1974. [44] Adam M Dziewonski and Don L Anderson. Preliminary reference earth model. Physics of the earth and planetary interiors, 25(4):297–356, 1981. [45] Göran Ekström and Adam M Dziewonski. The unique anisotropy of the pacific upper mantle. Nature, 394(6689):168, 1998. [46] Martin L Smith and FA Dahlen. The azimuthal dependence of love and rayleigh wave propagation in a slightly anisotropic medium. Journal of Geophysical Research, 78(17):3321–3333, 1973. [47] George Backus. A geometrical picture of anisotropic elastic tensors. Reviews of geophysics, 8(3):633–671, 1970. [48] Nicholas Harmon, Peter Gerstoft, Catherine A Rychert, Geoffrey A Abers, Mariela Salas de La Cruz, and Karen M Fischer. Phase velocities from seismic noise using beamforming and cross correlation in costa rica and nicaragua. Geophysical Research Letters, 35(19), 2008. [49] Beatriz Gaite, Arturo Iglesias, Antonio Villaseñor, Miguel Herraiz, and Javier F Pacheco. Crustal structure of mexico and surrounding regions from seismic ambient noise tomography. Geophysical Journal International, 188(3):1413–1424, 2012. [50] Francisco Córdoba Montiel, Arturo Iglesias Mendoza, Shri Krishna Singh, Zack Spica, and Denis Legrand. Tomografía de velocidad de grupo de ondas de rayleigh para el oriente de méxico y el istmo de tehuantepec. Boletín de la Sociedad Geológica Mexicana, 66(3):441–457, 2014. [51] Meijian An, Douglas A Wiens, Yue Zhao, Mei Feng, Andrew A Nyblade, Masaki Kanao, Yuansheng Li, Alessia Maggi, and Jean-Jacques Lévêque. S-velocity model and inferred moho topography beneath the antarctic plate from rayleigh waves. Journal of Geophysical Research: Solid Earth, 120(1):359–383, 2015. 35 [52] Kristoffer T Walker, Andrew A Nyblade, Simon L Klemperer, Götz HR Bokelmann, and Thomas J Owens. On the relationship between extension and anisotropy: Constraints from shear wave splitting across the east african plateau. Journal of Geophysical Research: Solid Earth, 109(B8), 2004. [53] Matthew J Fouch and Stéphane Rondenay. Seismic anisotropy beneath stable continental interiors. Physics of the Earth and Planetary Interiors, 158(2-4):292–320, 2006. [54] Adolphe Nicolas and Nikolas I Christensen. Formation of anisotropy in upper mantle peridotites-a review. Composition, structure and dynamics of the lithosphere-asthenosphere system, pages 111–123, 1987. [55] Vladislav Babuska and Michel Cara. Seismic anisotropy in the Earth, volume 10. Springer Science & Business Media, 1991. [56] Kairly Jaxybulatov, NM Shapiro, I Koulakov, A Mordret, M Landès, and C Sens-Schönfelder. A large magmatic sill complex beneath the toba caldera. science, 346(6209):617–619, 2014. [57] Matthew M Haney, T Dylan Mikesell, Kasper van Wijk, and Hisashi Nakahara. Extension of the spatial autocorrelation (spac) method to mixed-component correlations of surface waves. Geophysical Journal International, 191(1):189–206, 2012. [58] H Jödicke, A Jording, L Ferrari, J Arzate, Klaus Mezger, and Lars Rüpke. Fluid release from the subducted cocos plate and partial melting of the crust deduced from magnetotelluric studies in southern mexico: Implications for the generation of volcanism and subduction dynamics. Journal of Geophysical Research: Solid Earth, 111(B8), 2006. [59] Allen L Husker, Vladimir Kostoglodov, Victor M Cruz-Atienza, Denis Legrand, Nikolai M Shapiro, Juan S Payero, Michel Campillo, and Eduardo Huesca-Pérez. Temporal variations of non-volcanic tremor (nvt) locations in the mexican subduction zone: Finding the nvt sweet spot. Geochemistry, Geophysics, Geosystems, 13(3), 2012. [60] Víctor M Cruz-Atienza, Allen Husker, Denis Legrand, Emmanuel Caballero, and Vladimir Kostoglodov. Nonvolcanic tremor locations and mechanisms in guerrero, mexico, from energy-based and particle motion polarization analysis. Journal of Geophysical Research: Solid Earth, 120 (1):275–289, 2015. [61] Jorge Castellanos, Xyoli Pérez-Campos, Raúl Valenzuela, Allen Husker, and Luca Ferrari. Crust and upper-mantle seismic anisotropy variations from the coast to inland in central and southern mexico. Geophysical Journal International, 210(1):360–374, 2017. 36 [62] Don L Anderson. Recent evidence concerning the structure and composition of the earth’s mantle. Physics and Chemistry of the Earth, 6:1–131, 1965. [63] John D West, Matthew J Fouch, Jeffrey B Roth, and Linda T ElkinsTanton. Vertical mantle flow associated with a lithospheric drip beneath the great basin. Nature Geoscience, 2(6):439, 2009. [64] Luca Ferrari, Chiara M Petrone, and Lorella Francalanci. Generation of oceanic-island basalt–type volcanism in the western trans-mexican volcanic belt by slab rollback, asthenosphere infiltration, and variable flux melting. Geology, 29(6):507–510, 2001. [65] Donna K Blackman, J Kendall, et al. Seismic anisotropy in the upper mantle 2. predictions for current plate boundary flow models. Geochemistry, Geophysics, Geosystems, 3(9), 2002. [66] David Mainprice and Benoit Ildefonse. Seismic anisotropy of subduction zone minerals–contribution of hydrous phases. In Subduction zone geodynamics, pages 63–84. Springer, 2009. [67] Haemyeong Jung and Shun-ichiro Karato. Water-induced fabric transitions in olivine. Science, 293(5534):1460–1463, 2001. [68] Erik A Kneller, Peter E Van Keken, Shun-ichiro Karato, and Jeffrey Park. B-type olivine fabric in the mantle wedge: Insights from high-resolution non-newtonian subduction zone models. Earth and Planetary Science Letters, 237(3-4):781–797, 2005. [69] Vlad Constantin Manea, Marina Manea, Vladimir Kostoglodov, and Granville Sewell. Thermo-mechanical model of the mantle wedge in central mexican subduction zone and a blob tracing approach for the magma transport. Physics of the Earth and Planetary Interiors, 149(1-2):165– 186, 2005. [70] Leslie A Bernal-López, Berenice R Garibaldi, Gerardo León Soto, Raúl W Valenzuela, and Christian R Escudero. Seismic anisotropy and mantle flow driven by the cocos slab under southern mexico. Pure and Applied Geophysics, 173(10-11):3373–3393, 2016. [71] Gerardo León Soto, James F Ni, Stephen P Grand, Eric Sandvol, Raúl W Valenzuela, Marco Guzmán Speziale, Juan M Gómez González, and Tonatiuh Domínguez Reyes. Mantle flow in the rivera—cocos subduction zone. Geophysical Journal International, 179(2):1004–1012, 2009. [72] Ting Yang, Stephen P Grand, David Wilson, Marco Guzman-Speziale, Juan Martin Gomez-Gonzalez, Tonatiuh Dominguez-Reyes, and James 37 Ni. Seismic structure beneath the rivera subduction zone from finitefrequency seismic tomography. Journal of Geophysical Research: Solid Earth, 114(B1), 2009. [73] Jing-Yi Lin, Shu-Kun Hsu, and Jean-Claude Sibuet. Melting features along the western ryukyu slab edge (northeast taiwan): Tomographic evidence. Journal of Geophysical Research: Solid Earth, 109(B12), 2004. [74] Yildirim Dilek and Şafak Altunkaynak. Geochemical and temporal evolution of cenozoic magmatism in western turkey: mantle response to collision, slab break-off, and lithospheric tearing in an orogenic belt. Geological Society, London, Special Publications, 311(1):213–233, 2009. [75] Lee Siebert and Gerardo Carrasco-Núñez. Late-pleistocene to precolumbian behind-the-arc mafic volcanism in the eastern mexican volcanic belt; implications for future hazards. Journal of Volcanology and Geothermal Research, 115(1):179–205, 2002. [76] Arturo Gómez-Tuena, Alexandra B LaGatta, Charles H Langmuir, Steven L Goldstein, Fernando Ortega-Gutiérrez, and Gerardo CarrascoNúñez. Temporal control of subduction magmatism in the eastern transmexican volcanic belt: Mantle sources, slab contributions, and crustal contamination. Geochemistry, Geophysics, Geosystems, 4(8), 2003. [77] J Huw Davies and Friedhelm von Blanckenburg. Slab breakoff: a model of lithosphere detachment and its test in the magmatism and deformation of collisional orogens. Earth and Planetary Science Letters, 129(1-4):85–102, 1995. [78] Christele Guivel, Diego Morata, Ewan Pelleter, Felipe Espinoza, René C Maury, Yves Lagabrielle, Mireille Polvé, Hervé Bellon, Joseph Cotten, Mathieu Benoit, et al. Miocene to late quaternary patagonian basalts (46–47 s): geochronometric and geochemical evidence for slab tearing due to active spreading ridge subduction. Journal of Volcanology and Geothermal Research, 149(3):346–370, 2006. [79] Julia M Ribeiro, René C Maury, and Michel Grégoire. Are adakites slab melts or high-pressure fractionated mantle melts? Journal of Petrology, 57(5):839–862, 2016. [80] WP Schellart, Justin Freeman, DR Stegman, L Moresi, and David May. Evolution and diversity of subduction zones controlled by slab width. Nature, 446(7133):308, 2007. [81] Zohar Gvirtzman and Amos Nur. The formation of mount etna as the consequence of slab rollback. Nature, 401(6755):782, 1999. 38 [82] Stephen A Nelson and Erika González-Caver. Geology and k-ar dating of the tuxtla volcanic field, veracruz, mexico. Bulletin of Volcanology, 55 (1-2):85–96, 1992. [83] Surendra P Verma. Extension-related origin of magmas from a garnetbearing source in the los tuxtlas volcanic field, mexico. International Journal of Earth Sciences, 95(5):871, 2006. [84] Vlad C Manea, Marina Manea, and Luca Ferrari. A geodynamical perspective on the subduction of cocos and rivera plates beneath mexico and central america. Tectonophysics, 609:56–81, 2013. [85] Stephen Kirby, Robert E Engdahl, and Roger Denlinger. Intermediatedepth intraslab earthquakes and arc volcanism as physical expressions of crustal and uppermost mantle metamorphism in subducting slabs. Subduction top to bottom, pages 195–214, 1996. [86] Jian Wang, Dapeng Zhao, and Zhenxing Yao. Seismic anisotropy evidence for dehydration embrittlement triggering intermediate-depth earthquakes. Scientific reports, 7(1):2613, 2017. [87] Bradley R Hacker, Simon M Peacock, Geoffrey A Abers, and Stephen D Holloway. Subduction factory 2. are intermediate-depth earthquakes in subducting slabs linked to metamorphic dehydration reactions? Journal of Geophysical Research: Solid Earth, 108(B1), 2003. [88] Hallie E Meighan, Uri Brink, and Jay Pulliam. Slab tears and intermediate-depth seismicity. Geophysical Research Letters, 40(16):4244– 4248, 2013. [89] Bryan Isacks and Peter Molnar. Mantle earthquake mechanisms and the sinking of the lithosphere. Nature, 223(5211):1121, 1969. [90] Andrew D Bottrill, Jeroen van Hunen, and Mark B Allen. Insight into collision zone dynamics from topography: numerical modelling results and observations. Solid Earth, 3(2):387, 2012. [91] Masayuki Obayashi, Junko Yoshimitsu, and Yoshio Fukao. Tearing of stagnant slab. Science, 324(5931):1173–1175, 2009. [92] Paul E Sacks and Donald T Secor Jr. Delamination in collisional orogens. Geology, 18(10):999–1002, 1990. [93] R Govers and MJR Wortel. Lithosphere tearing at step faults: Response to edges of subduction zones. Earth and Planetary Science Letters, 236 (1-2):505–523, 2005. [94] Saskia Goes and Suzan van der Lee. Thermal structure of the north american uppermost mantle inferred from seismic tomography. Journal of Geophysical Research: Solid Earth, 107(B3), 2002. 39 Chapter 3 THE DYNAMIC INTERACTION BETWEEN THE CRUST AND THE MANTLE: A CASE STUDY IN THE PACIFIC NW. This chapter was adapted from: J. C. Castellanos, J. Perry-Houts, Robert W. Clayton, Y. Kim, A. C. Stanciu, B. Niday, and E. Humphreys (2020). “Seismic anisotropy reveals crustal flow driven by mantle vertical loading in the Pacific NW.” In: Science Advances 6.28. doi: 10.1126/sciadv.abb0476. 3.1 Abstract Buoyancy anomalies within the Earth’s mantle create different styles of convective currents that are thought to control the evolution of the lithosphere. While tectonic plate motions provide evidence for this relation, the mechanism by which mantle processes influence near-surface tectonics remains elusive. Here, we present an azimuthal anisotropy model for the Pacific NW crust that strongly correlates with high-velocity structures in the underlying mantle, but shows no association with the regional mantle flow field. We suggest that the crustal anisotropy is decoupled from horizontal basal tractions and instead is created by upper mantle vertical loading, which in turn generates pressure gradients that drive channelized flow in the ductile mid-lower crust. We then demonstrate the interplay between small-scale mantle heterogeneities and lithosphere dynamics by predicting the viscous crustal flow that is driven by local buoyancy sources within the upper mantle. Our findings reveal how mantle vertical load distribution, acting independently from horizontal mantle flow, can actively control crustal deformation on a scale of several hundred kilometers. 3.2 Introduction Geodynamic models commonly describe the relation between buoyancy-driven mantle convection and plate tectonics with two components of traction applied to the base of the lithosphere—vertical tractions giving rise to dynamic topography [1, 2], and horizontal basal tractions driving plate motion and tectonic deformation [3]. Although lithospheric stress field measurements [4, 5] and 40 mantle flow patterns [6, 7] provide critical constraints on the dynamics of these interactions, attempts to isolate their relative influence on near-surface tectonics often yield ambiguous results. In most cases, the difficulty derives from our imperfect knowledge of the mantle density structure and the high variability in the material strength of the lithospheric rocks, which greatly influences the degree of mechanical coupling between the tectonic plates and the convective flow [8]. While substantial advancements in seismic imaging have permitted the construction of high-resolution models of the mantle’s mass distribution, an ability to accurately quantifying the degree of coupling between the mantle and the lithosphere remains underdeveloped. This limitation, in combination with a paucity of observational constraints, has prevented any reliable assessment of how mantle-based forces interact with plate-scale processes to give rise to the tectonic stresses that drive surface deformation. Here we show that, under certain rheological conditions, crustal anisotropy is transparent to the structural complications of the crust and can reveal a crustal flow driven by the vertical coupling of the mantle and the lithosphere. The process for detecting mantle-induced vertical deformation generally involves identifying regions that have experienced rapid surface uplift or subsidence [9, 10] or areas with sharp elevation contrasts that are difficult to explain with simple isostatic models [11]. The Wallowa Mountain block in northeastern Oregon for example represents a remarkable regime where both mantle- and crustal-based stresses appear to have played an essential role in lithosphere dynamics. These mountains are comprised of a sizable granitic batholith that rapidly rose 2 km above the surrounding area shortly after the deposition of the Columbia River flood basalts (CRB) 16Ma [12, 13], creating an impressive topographic bullseye centered on the Wallowa Mountains (Figure 3.1). The compact and isolated uplift of the granitic Wallowa batholith suggests the foundering [13] of a dense garnet-rich [14] pluton root during or shortly after the CRB eruptions. However, the regional post-CRB uplift of the entire topographic bullseye region [14] indicates the existence of a larger-scale mechanism that dynamically drives crustal deformation around the site of the inferred foundering event. In the mantle beneath the Wallowa Mountains, high-resolution tomographic images persistently reveal the presence of a major high-velocity anomaly (the Wallowa anomaly) that is circular in map view and extends to a depth of 350 km [15, 16]. This structure appears to be part of a system of ancient slab fragments that are dangling beneath the Pacific NW 41 and, together with buoyant plumes of rising asthenosphere, is hypothesized to drive the small-scale mantle convection that actively modifies the western U.S. lithosphere [16]. Although the precise role of the Wallowa anomaly and other nearby mantle heterogeneities in shaping the topography of NE Oregon remains unclear, recent seismic imaging studies have revealed that the crust just north of the Wallowas is about 20 km thicker than the surrounding area [17, 18]. The correlation between these two puzzling features seems to suggest that the negatively buoyant Wallowa anomaly is responsible for the localized pull-down on the Moho. Figure 3.1: Regional and location maps for northeastern (NE) Oregon. (A) Regional map showing the broadband seismic stations used in this study (black inverted triangles). The dashed blue line depicts the Snake River Plain (SRP). WA, Washington; OR, Oregon; ID, Idaho. (B) Global map centered in NE Oregon. The thick black line encloses the region shown in (A). (C) Elevation map of the topographic bullseye region (red area in A). The dashed larger ellipse is the outer limit of the bullseye, whereas the inner ellipse locates the Wallowa batholith. In this study, we developed an azimuthal anisotropy model for the crust of NE Oregon and its surrounding regions using short-period (3-17 s) Rayleigh waves extracted from ambient seismic noise cross-correlations (Figure 3.2). With the 42 concentration of broadband stations in the area and wide azimuthal interstation path coverage, we can reliably resolve the lateral variations of seismic anisotropy for the uppermost 35 km of the crust. To measure the anisotropy, we implemented a beamforming scheme that allows us to characterize the seismic wavefield’s velocity dependence with propagation direction beneath each station [19, 20]. The reliability of our model was then verified by comparing our surface wave anisotropy measurements to those that were obtained by characterizing the azimuthal dependence of receiver functions at stations surrounding the Wallowa Mountains. Figure 3.2: Example of beamformer outputs and final azimuthal anisotropy model. (A) Two-dimensional histograms over the azimuth-velocity space for the 3- to 17-s period band with the best-fitting anisotropy model (red dashed lines). The green bars on top of each panel indicate the number of noise cross-correlations available for each azimuth. (B) Azimuthal anisotropy model for the crust of the Pacific NW. Bar orientation gives the fast direction of azimuthal anisotropy, and bar length is proportional to anisotropy amplitude. The background color represents the intersection density of the anisotropy vectors assuming that they are of infinite length (i.e., projected to the bounds of the study region). The green and blue dots indicate the location of the two stations beamformed in (A). 43 3.3 Results The azimuthal anisotropy model for the crust underlying this region does not correlate with the surface geology, structural trends, or the mapped crustal stress field [21]. However, it instead holds a remarkable connection to the upper mantle velocity distribution (Figure 3.3A). The fast directions of anisotropy, which are generally thought to reflect the coherent deformation of small-scale structures and preferred alignment of anisotropic minerals [22], show a simple and well-defined radial pattern that strongly correlates with the Wallowa anomaly. Moreover, the northern and easternmost anisotropy vectors display a subtle fan-like pattern that correlates strikingly well with the geographic extent of the Siletzia slab curtain beneath Idaho [15, 16]. The amplitude of the azimuthal anisotropy also decreases to near-zero values for seismic stations above the Wallowa anomaly and slab curtain kink, where the geometry of anisotropy may be transitioning into one that is null in the horizontal plane. The connection between crustal anisotropy and upper mantle velocity structure suggests that mantle gravitational loads actively induce vertical stresses on the overlying material and, through this relation, control crustal deformation in NE Oregon and its adjacent regions. Based on the spatial coherence of our measurements, we hypothesize that the azimuthal anisotropy in NE Oregon results from the lattice preferred orientation (LPO) of anisotropic minerals with the subhorizontal flow of the mid-lower crust [23, 24]. Recent numerical studies show that stresses transmitted upwards from the underlying mantle can induce significant amounts of intraplate deformation through Poiseuille and Couette flow due to lateral pressure variations and basal shear [25]. This style of deformation requires the lithospheric rocks to have low viscous strength in order to form a channeled ductile flow system in the mid-lower crust that decouples the upper crustal and upper mantle stress fields. Because of the relatively recent magmatic activity in NE Oregon, the crust beneath this region can achieve the adequate thermal conditions (700-1000 °C) to create such a ductile and mobile environment, especially along the Snake River Plain and beneath the Wallowa Mountains [26, 27]. The existence of a mid-lower crustal weak channel would then allow the crust to flow in response to lateral pressure gradients and accommodate the vertical stresses exerted by the underlying mantle. Because of the high local Moho temperature, the mapped crustal anisotropy is most likely to be dominated by the LPO of type II and III fabrics in amphibole, for which the 44 Figure 3.3: Comparison of seismic and geodynamic results with crustal anisotropy. (A) Azimuthal anisotropy model for the crust of the Pacific NW overlying a depth slice through the Vp tomography model at 250 km [16]. The red dashed lines depict the Wallowa anomaly and the Siletzia slab curtain. (B) Modeled Moho stress and mid-lower crustal flow velocity for the Pacific NW. The colored contours represent the vertical stress at the Moho based on a global geodynamic model driven by density anomalies derived from the P-wave velocity structure. The black arrows denote the predicted mid-lower crustal flow velocity that results from the application of the modeled Moho stress to a viscously heterogeneous crust. The red bars represent the anisotropy measurements derived from this study. fast direction of anisotropy is sub-parallel to the flow direction [28]. The alignment of mica crystals may also contribute to the overall observed anisotropy; however, global compilations of the structure of the continental crust suggest that its deep portion contains rather little mica, and that amphibole takes a larger fraction of its composition [29]. At sublithospheric depths, shear-wave splitting observations reveal that the asthenospheric flow in the Pacific NW is primarily controlled by a combination of North American plate motion and the sinking of the Juan de Fuca and Farallon slab systems [30, 31]. These measurements also reveal that there is little, if any, mantle deformation caused by the downward movement of the 45 Wallowa anomaly. As a matter of fact, the most recent SKS splitting observations in NE Oregon indicate that mantle materials flow smoothly around the lateral boundaries of the Wallowa anomaly rather than converging on the site of lithospheric load [32] (Figure 3.4). The lack of a strong disturbance in the mantle flow field beneath this area suggests that the asthenospheric strain that is created by the downwelling velocity of the Wallowa anomaly is not strong enough to perturb the current LPO that has been established by the long-term movement of the tectonic plates. This observation leads us to the notion that whatever vertical forces are being exerted by the upper mantle and driving the crustal flow are almost entirely derived from the negative buoyancy of its dense structures rather than the weak vertical asthenosphere flow that is excited by their vertical movement. Furthermore, the absence of correlation between the crustal and the upper mantle deformation fields suggests that the mantle lithospheric strength isolates horizontal asthenospheric flow from that in the crust such that there is insignificant basal shearing by the underlying mantle (i.e. mantle flow is not driving Couette flow in the viscous crust). This type of decoupling is consistent with previous tomographic findings that show a weak correlation between crustal and upper mantle anisotropy in most regions of the western U.S. [33]. Here, it is important to note that a key to the low strain rate of the dense mantle structures is their greater viscosity, and that their low sinking rate is a result of them being attached to the North American lithosphere, as seismically imaged [16, 17]. In the ductile regime, viscous strain rate preferentially orients minerals relative to one another, generating seismic anisotropy that is aligned with the flow direction [36]. Within this framework, we model the crustal deformation induced by mantle loading through scaling the seismically-imaged mantle velocity anomalies [16] to density structure [37] and predicting vertical tractions at Moho depths. Here, we exclude the subducting Juan de Fuca slab and North American craton since these structures have been in long-term steady state relative to the more recent mantle structures of the interior Pacific NW and are, therefore, unlikely to play a crucial role in shaping the present-day crustal strain field. We also impose a small-scale load, with moderate stress magnitude, to the predicted Moho traction in the Wallowas area to account for the foundering of the mountain’s pluton root. This last step is taken because this event is a rather short-lived transient phenomenon that is not captured by the seismic tomography and may also contribute to the observed anisotropy; the 46 Figure 3.4: Station averaged shear wave splitting measurements for the Pacific NW. (A) SKS splitting measurements for the entire western U.S. [31]. The red arrow depicts the relative motion between the North American plate (NA) and the hotspot reference frame (HS) [34]. (B) SKS splitting measurements for our study region (red area in A). The thick blue vectors depict the measurements of Niday and Humphreys [32], and the black vectors are from the database of Becker et al. [31]. The orientation of the vectors gives the angle of the fast polarization, and the length of the bars is proportional to the magnitude of the shear wave splitting. The white trajectories through the anisotropy field lines in (B) are used to represent the streamlines of the mantle flow assuming an east-oriented flow [35]. The red region in the background marks the location of the Wallowa anomaly. Note how mantle materials appear to flow smoothly around the lateral boundaries of the Wallowa anomaly. destabilization and subsequent removal of the root would cause the weak midlower crust around the Wallowas to flow towards the vacated root region [25]. Moho traction calculations made without incorporating the localized Wallowa load reveal that the first-order vertical stress distribution is not significantly altered. The final mantle-derived vertical tractions are then used to drive viscous Stokes flow in a thin crustal channel using surface heat flow as a modulator of crustal viscosity [38]. Figure 3.3B shows a comparison between the observed crustal anisotropy, the estimated vertical stress at the Moho, and the 47 predicted mid-lower crustal flow. Both the relative amplitude and orientation of the crustal flow velocity agree remarkably well with the measured crustal anisotropy within the main study area, displaying a dominant radial flow pattern centralized at the Wallowa Mountains site. Such flow would lead to crustal thickening in the low-pressure regions, which is consistent with the nearly circular ∼20-km Moho depression that is observed above the Wallowa mantle anomaly [17, 18]. Note, however, that there exists a disagreement between our modeled and measured anisotropy in the southwestern part of our study region. This discrepancy may be the result of either non-mantle processes that are not included in our numerical model (i.e. tectonic strain in the active Basin and Range), or due to the inherent shortcoming of the beamforming technique to resolve lateral sharp changes in the anisotropic structure. Nonetheless, the general agreement between the crustal flow predicted by our simple model and the seismic observables strongly suggests that mantle-induced stresses can, in some cases, have more significant control on intraplate deformation than those transmitted laterally from nearby active plate boundaries [2]. 3.4 Discussion and Conclusion On the basis that the lithosphere is rheologically stratified, we propose an upside-down water-bed model in which vertical mantle loads cause the ductile rocks inside the weak mechanical layer to migrate horizontally towards low-pressure regions through Poiseuille flow, involving little mantle deformation (Figure 3.5). The mechanics of the channelized flow that is induced by this model are similar to the ones that are typically invoked to explain nearsurface deformation in extreme tectonic environments such as the Tibetan plateau [37] or the Altiplano in the Bolivian Andes [39]. The difference, however, lies on the fact that crustal flow in such regimes is generally thought to be driven tectonically or gravitationally as a response of the buoyancy forces that arise from differential crustal densities [40] or the pressure differences caused by varying crustal thickness [41]. Therefore, the evidence that mantle gravitational loads are capable of displacing weak crustal materials in a comparable manner not only refines our understanding of the interaction between crustal tectonics and mantle dynamics, but also brings to light another source of deformation that might be necessary to explain the state of stress in the crust of other tectonically enigmatic regions. Figure 3.6 for instance shows the crustal anisotropy measurements of Lin et al. [33] around south and 48 central California, where other dense mantle anomalies have been imaged by different tomographic studies. Similar to the case of the Wallowas, the strong correlation between the crustal anisotropy and upper mantle velocity structure in this region suggests that the upside-down water-bed model is playing a crucial role in driving the evolution of the lithosphere. However, because of the unique tectonic of California, deciphering the precise role and contribution of mantle-based stresses in its surficial processes would require a more complete modeling that incorporates the rather elevated horizontal strains that are exerted from the active plate margin. Figure 3.5: Schematic representation of the upside-down water-bed model. The load of the mantle lithosphere is a force creating vertical stresses (σzz ) on the Moho. The lithospheric load pulls down on the crust, which creates a lateral pressure gradient that drives Poiseuille flow in the ductile mid-lower crust. The asthenosphere flows independently (as evidenced by its independent anisotropy field; Figure 3.4), creating a local Couette flow that is decoupled from the mid-lower crust by the mantle lithosphere. In general, the crustal anisotropy that results from the upside-down water bed model is most sensitive to recent deformation, and hence relevant for addressing young tectonic evolution. For the case of NE Oregon, the temporal sequence of the Siletzia slab curtain formation 53 Ma and the Wallowa anomaly delamination 16 Ma [15] may well explain the apparent dominance of the Wallowa anomaly in aligning the fast directions of anisotropy. This argument is supported by the fact that, although the additional load imposed at the site of the Wallowa batholith was initially designed to incorporate the effects of the foundering of its root, the localized Wallowa enhancement would still be required to represent the latest stage of upper mantle vertical forcing and achieve the remarkable centralized flow pattern that is illuminated 49 Figure 3.6: Crustal anisotropy and upper mantle velocity structure of California. The black vectors depict the surface wave anisotropy measurements at 12 s period [33]. Bar orientation gives the fast direction of azimuthal anisotropy, and bar length is proportional to anisotropy amplitude. The background color corresponds to a depth slice through the Vp tomography model at 195 km [16]. The red dashed lines denote the two seismically fast and likely dense mantle anomalies. by the crustal anisotropy. The flow that is predicted by mantle loading alone thus suggests that the amount of strain exerted by the delamination of the Wallowa mantle anomaly is enough to effectively align mid-crustal minerals and overprint on any pre-strained fabric. An alternative and simple interpretation is that the viscous strength of the crustal rocks beneath and around the Wallowa Mountains is weaker because of the recent CRB eruptions [27], and have deformed easily, flowing toward the site of the mantle loading and observed crustal thickening. Regardless of the relative level of contribution of each of these mechanisms, our findings provide strong observational evidence of regional-scale mantle-crust vertical coupling and highlight the fundamental importance of upper-mantle buoyancy in understanding near-surface tectonics. 50 References [1] Bradford H Hager, Robert W Clayton, Mark A Richards, Robert P Comer, and Adam M Dziewonski. Lower mantle heterogeneity, dynamic topography and the geoid. Nature, 313(6003):541–545, 1985. [2] Emily A Neil and Gregory A Houseman. Rayleigh–taylor instability of the upper mantle and its role in intraplate orogeny. Geophysical Journal International, 138(1):89–107, 1999. [3] NJ Kusznir and RG Park. Intraplate lithosphere deformation and the strength of the lithosphere. Geophysical Journal International, 79(2): 513–538, 1984. [4] Bernhard Steinberger, Harro Schmeling, and Gabriele Marquart. Largescale lithospheric stress field and topography induced by global mantle circulation. Earth and Planetary Science Letters, 186(1):75–91, 2001. [5] Peter Bird, Zhen Liu, and William Kurt Rucker. Stresses that drive the plates from below: Definitions, computational path, model optimization, and error analysis. Journal of Geophysical Research: Solid Earth, 113 (B11), 2008. [6] Carolina Lithgow-Bertelloni and Paul G Silver. Dynamic topography, plate driving forces and the african superswell. Nature, 395(6699):269– 272, 1998. [7] Thorsten W Becker and Richard J O’Connell. Predicting plate velocities with mantle circulation models. Geochemistry, Geophysics, Geosystems, 2(12), 2001. [8] Attreyee Ghosh and William E Holt. Plate motions and stresses from global dynamic models. Science, 335(6070):838–843, 2012. [9] MJ Hoggard, Nicholas White, and David Al-Attar. Global dynamic topography observations reveal limited influence of large-scale mantle flow. Nature Geoscience, 9(6):456–463, 2016. [10] Michael Gurnis. Rapid continental subsidence following the initiation and evolution of subduction. Science, 255(5051):1556–1558, 1992. [11] Jean Braun. The many surface expressions of mantle dynamics. Nature Geoscience, 3(12):825–833, 2010. [12] JA Wolff, FC Ramos, SP Reidel, VE Camp, ME Ross, BS Martin, TL Tolan, and RE Wells. Source materials for the main phase of the columbia river basalt group: Geochemical evidence and implications for magma storage and transport. The Columbia River Flood Basalt Province: Geological Society of America Special Paper, 497:273–291, 2013. 51 [13] TC Hales, DL Abt, ED Humphreys, and Joshua J Roering. A lithospheric instability origin for columbia river flood basalts and wallowa mountains uplift in northeast oregon. Nature, 438(7069):842–845, 2005. [14] Leland J O’Driscoll and K Johnson. Geochemical evidence for a dense mafic root of units in the wallowa batholith, wallowa mountains, northeast oregon. In GSA Annual Meeting, Houston, Texas, volume 9, 2008. [15] Amberlee Darold and Eugene Humphreys. Upper mantle seismic structure beneath the pacific northwest: A plume-triggered delamination origin for the columbia river flood basalt eruptions. Earth and Planetary Science Letters, 365:232–242, 2013. [16] Brandon Schmandt and Eugene Humphreys. Complex subduction and small-scale convection revealed by body-wave tomography of the western united states upper mantle. Earth and Planetary Science Letters, 297 (3-4):435–445, 2010. [17] Haiying Gao, Eugene D Humphreys, Huajian Yao, and Robert D van der Hilst. Crust and lithosphere structure of the northwestern us with ambient noise tomography: Terrane accretion and cascade arc development. Earth and Planetary Science Letters, 304(1-2):202–211, 2011. [18] Haiying Gao. Crustal seismic structure beneath the source area of the columbia river flood basalt: Bifurcation of the moho driven by lithosphere delamination. Geophysical Research Letters, 42(22):9764–9771, 2015. [19] Nima Riahi, Götz Bokelmann, Paola Sala, and Erik H Saenger. Timelapse analysis of ambient surface wave anisotropy: A three-component array study above an underground gas storage. Journal of Geophysical Research: Solid Earth, 118(10):5339–5351, 2013. [20] Katrin Löer, Nima Riahi, and Erik H Saenger. Three-component ambient noise beamforming in the parkfield area. Geophysical Journal International, 213(3):1478–1491, 2018. [21] Oliver Heidbach, Mojtaba Rajabi, Xiaofeng Cui, Karl Fuchs, Birgit Müller, John Reinecker, Karsten Reiter, Mark Tingay, Friedemann Wenzel, Furen Xie, et al. The world stress map database release 2016: Crustal stress pattern across scales. Tectonophysics, 744:484–498, 2018. [22] Vladislav Babuska and Michel Cara. Seismic anisotropy in the Earth, volume 10. Springer Science & Business Media, 1991. [23] Dan McKenzie and James Jackson. Conditions for flow in the continental crust. Tectonics, 21(6):5–1, 2002. 52 [24] Christopher Beaumont, Rebecca A Jamieson, Mai H Nguyen, and Sergei Medvedev. Crustal channel flows: 1. numerical models with applications to the tectonics of the himalayan-tibetan orogen. Journal of Geophysical Research: Solid Earth, 109(B6), 2004. [25] Huilin Wang and Claire A Currie. Crustal deformation induced by mantle dynamics: Insights from models of gravitational lithosphere removal. Geophysical Journal International, 210(2):1070–1091, 2017. [26] Nadine McQuarrie and David W Rodgers. Subsidence of a volcanic basin by flexure and lower crustal flow: The eastern snake river plain, idaho. Tectonics, 17(2):203–220, 1998. [27] Derek L Schutt, Anthony R Lowry, and Janine S Buehler. Moho temperature and mobility of lower crust in the western united states. Geology, 46(3):219–222, 2018. [28] Byeongkwan Ko and Haemyeong Jung. Crystal preferred orientation of an amphibole experimentally deformed by simple shear. Nature communications, 6(1):1–10, 2015. [29] DJ Tatham, GE Lloyd, RWH Butler, and M Casey. Amphibole and lower crustal seismic properties. Earth and Planetary Science Letters, 267(1-2): 118–128, 2008. [30] Federica Marone and Barbara Romanowicz. The depth distribution of azimuthal anisotropy in the continental upper mantle. Nature, 447(7141): 198–201, 2007. [31] Th W Becker, S Lebedev, and MD Long. On the relationship between azimuthal anisotropy from shear wave splitting and surface wave tomography. Journal of Geophysical Research: Solid Earth, 117(B1), 2012. [32] William Niday and Eugene Humphreys. Complex upper mantle anisotropy in the pacific northwest: Evidence from sks splitting. Earth and Planetary Science Letters, 540:116264, 2020. [33] Fan-Chi Lin, Michael H Ritzwoller, Yingjie Yang, Morgan P Moschetti, and Matthew J Fouch. Complex and variable crustal and uppermost mantle seismic anisotropy in the western united states. Nature Geoscience, 4(1):55–61, 2011. [34] Alice E Gripp and Richard G Gordon. Young tracks of hotspots and current plate velocities. Geophysical Journal International, 150(2):321– 361, 2002. [35] PG Silver and WE Holt. The mantle flow field beneath western north america. Science, 295(5557):1054–1057, 2002. 53 [36] Jeffrey Park and Vadim Levin. Seismic anisotropy: tracing plate dynamics in the mantle. Science, 296(5567):485–489, 2002. [37] Leigh H Royden, B Clark Burchfiel, Robert W King, Erchie Wang, Zhiliang Chen, Feng Shen, and Yuping Liu. Surface deformation and lower crustal flow in eastern tibet. science, 276(5313):788–790, 1997. [38] David Blackwell, Maria Richards, Zachary Frone, Joe Batir, Andrés Ruzo, Ryan Dingwall, and Mitchell Williams. Temperature-at-depth maps for the conterminous us and geothermal resource estimates. Technical report, Southern Methodist University Geothermal Laboratory, Dallas, TX (United States), 2011. [39] L Husson and Thierry Sempere. Thickening the altiplano crust by gravitydriven crustal channel flow. Geophysical Research Letters, 30(5), 2003. [40] Susan Ellis. Crustal flow modes in large hot orogens. Channel flow, ductile extrusion and exhumation in continental collision zones, (268):91, 2006. [41] Marin Kristen Clark and Leigh Handy Royden. Topographic ooze: Building the eastern margin of tibet by lower crustal flow. Geology, 28(8): 703–706, 2000. 54 Chapter 4 THE DYNAMIC INTERACTION BETWEEN THE CRUST AND THE MANTLE: IN SEARCH FOR THE UPSIDE-DOWN WATER BED MODEL. This chapter was adapted from: J. C. Castellanos, E. Humphreys, and R. W. Clayton (In review). “Evidence for mantle-based deformation across the western US.” In: Geophysical Research Letters. 4.1 Abstract We investigate the role that upper mantle buoyancy anomalies play in determining the behavior of the crust. Recently, Castellanos et al. [1] observed that the anisotropy of the Pacific NW crust correlates with the upper-mantle velocity structure and, from this connection, suggested that vertical loads in the upper mantle can drive crustal flow on a regional scale. To provide further insight into this relation, we resolve the crustal anisotropy in regions where near-surface mantle-based deformation might be have occurred or is presently occurring. Specifically, we focus on the crust around the Rocky Mountains and around California since high-resolution tomographic images reveal the presence of mantle structures similar to the ones that are thought to be driving the crust in the Pacific NW. Our results reveal crustal flow driven by mantle vertical loading in both regions and suggest that this mechanism may be key in maintaining crustal isostasy during an orogeny. 4.2 Introduction Unraveling the influence of deep geodynamic processes on the Earth’s surface stress field is critical for understanding the driving forces of tectonic deformation. To date, it is well-established that there are two main sources of stress in the lithosphere: (1) internal buoyancy forces arising from lateral compositional variations within the crust and lithospheric mantle [e.g., density and thickness; 2, 3], and (2) vertical and horizontal basal tractions arising from buoyancydriven mantle convection below the lithosphere [4, 5]. While substantial work 55 has been done to define the kinematics of these two sources, their relative contribution on both the long-term stability of continents and their state of stress is largely unknown. Here, we investigate how mantle-based stresses affect the dynamics of the lithosphere through the analysis of crustal anisotropy. In general, the difficulty of elucidating the origin of lithospheric stresses stems from our imperfect knowledge of the physical properties of the crust and the lack of constraints on the degree of coupling between the tectonic plates and the convective flow of the mantle. Over the last few decades, numerous studies have aimed at constraining the compositional structure of the crust. These efforts typically involve the modeling of the Earth’s topographic response to tectonic loading [e.g., 6, 7] or the use of earthquake data [e.g., 8] to derive estimates of crustal viscosity and temperature. Findings show, for instance, that there can exist large compositional lateral variations across a single craton [e.g., 9], and that certain regions around the world have the conditions for the lower crust to act as a weak viscous layer capable of accommodating the lateral pressure gradients within the lithosphere [e.g., 10, 11]. Methods aimed at constraining the degree of coupling between the crust and the mantle have also been developed. These techniques include modeling the lithospheric stress field and the prediction of the tectonic plate’s motion using mantle circulation models [e.g., 5, 12]. Although these investigations have made tremendous advancements in explaining a large part of the Earth’s surface observables, there are still many regions where agreement between observed and predicted stresses is poor. These discrepancies may reflect the existence of sources of deformation that are not accounted for in the modeling, or the fact that the crust might have sections of complex layering such that stresses measured at the surface do not necessarily represent the state of stress of the entire plate [e.g., 13]. Consequently, the possibility of these scenarios demand for more observational data capable of constraining the orientation and depth of deformation in the Earth. Seismic anisotropy–the dependence of seismic wave speeds on propagation direction and polarization–is a useful tool to detect the existence of coherently deformed structures and preferentially oriented anisotropic minerals [14, 15]. As such, the characterization of seismic anisotropy represents an ideal observational method to map how stress is accommodated within the Earth. In a previous study, Lin et al. [16] used stations from the EarthScope Transportable 56 Array (TA) to investigate the crustal and mantle anisotropy beneath the western US. In their analysis, they confirm that the asthenospheric flow beneath this region is primarily controlled by a combination of North American plate motion and the sinking of the Juan de Fuca and Farallon slabs. But, more interestingly, they observed that mid- to lower-crustal seismic anisotropy is regionally coherent yet largely uncorrelated with the mantle anisotropy, suggesting that these two layers deformed independently. While the cause of the crustal anisotropy is agreed to be the result of the latticed preferred orientation of anisotropic minerals, the mechanisms that underlie the creation of some of the observed patterns remain puzzling. In a more recent study, Castellanos et al. [1] suggest vertical loads in the upper mantle can drive crustal deformation on a regional scale that is independent of the asthenosphere’s flow field. In detail, they proposed that the connection between mantle density structure and surface deformation occurs because the mantle buoyancy-created stresses displace the Moho and drive viscous flow in the ductile lower crust. One example of this mechanism was found in the Pacific northwest (PNW; dashed region in Figure 4.1), where a simple anisotropy field is centered on a well-developed local Moho depression of 15 km beneath the eastern end of the Washington-Oregon state line [17, 18]. The crustal down-warp is itself centered above the seismically fast, compact “Wallowa” mantle anomaly [19], and is calculated to hold ∼150,000 km3 of crust, presumably supplied by crustal inflow. The surrounding crustal anisotropy field is radial in form, extending to 250 km from the anomaly center. To obtain a broader perspective on mantle-driven viscous crustal flow, we extend the analyses of Lin et al. [16] and Castellanos et al. [1] and use the dense station coverage provided by the TA, and a few other temporal networks, to reliably resolve the lateral variations of crustal anisotropy in regions where near-surface mantle-based deformation might have occurred or is presently occurring. In particular, we focus on areas around the Rocky Mountains and California since high-resolution tomographic images [20] reveal the presence of seismically fast mantle structures similar to the ones that are thought to have driven crustal flow in the PNW (Figure 4.1). 4.3 Data and Methods We cross-correlate the ambient seismic noise field recorded by the TA with the one recorded by every publicly available broadband seismometer that operated 57 Figure 4.1: Upper mantle velocity structure of the western US and geography of the target regions. The main map shows a depth slice through the Vp tomography model at 195 km [20]. The dashed black rectangle delimits the region of study of Castellanos et al. [1], whereas the continuous black rectangles delimit the regions that are focused in this investigation. The zoomed maps show the distribution of broadband seismic stations used to resolve the crustal anisotropy beneath each region (red inverted triangles). within our two regions of interest. This process allowed us to extract clear fundamental mode surface waves traveling between different pairs of stations that were then used to construct the azimuthal anisotropy models. Because of the large differences in the crustal properties between the western and central US, we performed our surface wave study in two different frequency bands. The limits of these filters were determined by using a 3-D velocity model of the US [Figure 4.2AB; 21], and computing the Rayleigh wave sensitivity kernels for a wide range of frequencies at each (x-y) coordinate of the domain. We then used a 3-D crustal thickness model of the US [Figure 4.2C; 22] to define the depth of the middle and lower crust at each location and, with it, find the shortest and longest period in which the majority of the sensitivity kernel’s amplitude lies between the two layers. Figure 4.2EF shows the shortest period that is sensitive to the middle crust and the longest period that is sensitive to the lower crust, respectively. With this analysis, we determine that, in order to map the average anisotropic properties of the lower crust, the surface wave analyses must be made between the 18-31 s period band for the Rocky Mountain Complex, and between the 12-23 s period band for the California area. 58 Figure 4.2: Crustal properties of the western US crust. (A-B) Show the midcrust and lower crust shear wave velocity structure as derived from ambient seismic noise and earthquake tomography, receiver functions, and Rayleigh wave ellipticity (H/V) measurements [21]. (C) Shows the crustal thickness as derived from teleseismic P-to-S receiver functions, and Rayleigh wave phase velocities and ellipticity from noise interferometry and earthquakes [22]. (D) Shows the surface topography [23]. (E-D) Shows the shortest Rayleigh wave period that is sensitive to the middle crust and the longest Rayleigh period that is sensitive to the lower crust, respectively. The average upper and lower period for the two target regions is shown in the bottom of each of the geographic borders. To build the anisotropy models, we adopt a seismic beamforming scheme [24]. Within this framework, we use every station that is in our target areas and create several subarrays so that the coherent energy moving through each group of stations can be translated into a local phase velocity and direction of propagation. Here, the number of instruments that composed each subarray 59 is variable, but their radii, and hence their resolution, are close to constant as we set a 1-wavelength minimum and a 2-wavelength maximum threshold, thereby avoiding spatial aliasing [25, 26]. To ensure the exclusive use of highquality waveform data, we only use cross-correlation functions that have a SNR larger than 5 and an offset of at least 3-wavelengths away from the geographic center of the subarrays. Once the band-passed energy of all the virtual sources are beamformed, we collect all phase velocity and backazimuth measurements that were made at each subarray and characterize the wavefield’s azimuthal dependence using the generalized model of Smith and Dahlen [27] for surface waves in a weakly anisotropic media. This parametrization allows us to extract a fast direction term and an amplitude term that are, in principle, related to the anisotropic properties of the lower crust. 4.4 Results Figure 4.3AC shows the lower crust azimuthal anisotropy model for both target areas plotted on top of the mantle P-wave structure at 195-km depth. The bars are oriented in the fast azimuth direction and their length is proportional to the magnitude of anisotropy. In a similar representation, Figure 4.3BD, shows the SKS-derived upper mantle anisotropy measurements compiled by Becker et al. [28]. Below we analyze the relation between these two anisotropy fields and discuss the mechanisms which may underlie the creation of the lower crust anisotropy of both regions. Mantle anisotropy beneath the Rocky Mountain region appears to be fully dissociated from the anisotropy in the overlying crust. Nonetheless, and similar to that seen in the PNW, the crustal anisotropy and upper mantle seismic velocity structures seem to relate to one another (Figure 3A). Fast azimuth orientations tend to be perpendicular to a SW-trending seismically fast upper mantle structure beneath Wyoming and parts of Colorado and Utah, with near-zero amplitudes at almost its geographic center. Both of these features suggests that the anisotropy of this region was created by crustal flow roughly perpendicular to this trend. Consideration of the tectonic history of the crust around Rocky Mountains offers insight into the creation of its anisotropy field. After residing near sea level for 100s of m.y., the swath of crust that now resides above a SW-trending zone of deep fast mantle beneath Wyoming became anomalously depressed in 60 Figure 4.3: Lower crust and mantle anisotropy for the two target regions. (A) Lower crust azimuthal anisotropy around the Rocky Mountains. (B) Stationaveraged shear-wave splitting measurements around the Rocky Mountains [28]. (C) Lower crust azimuthal anisotropy around California. (D) Station-averaged shear-wave splitting measurements around California [28]. The bar orientation gives the fast direction of anisotropy, and the bar length is proportional to anisotropy amplitude. The background color represents the mantle P-wave structure at 195 km depth [20]. In panels (A) and (B) the green arrows depict our preferred interpretation of crustal flow. The blue arrows in (C) depict the relative motions between the Pacific (PA) and North American plates (NA), the Juan de Fuca plate (JdF) and NA. The purple lines in (C) mark the plate boundaries. ISA stands for the Isabella mantle anomaly. 61 the mid-Cretaceous. This event is evidenced by a series of deep marine basins that represent a local mantle loading along this trend [29, 30]. In this setting, and assuming that the rheological conditions of the crust had sufficiently low viscosity, crustal materials could have flowed towards the Moho depression and created the anisotropy that is observed at present (Figure 4.4A). In an alternate scenario, subsequent uplift of the region represents an increase in mantle buoyancy that progressively elevated the continental interior–including the depressed trend–into the broad western US uplift (Figure 4.2D). As a result, the trend of Cretaceous basins became no longer a zone of anomalous topography or crustal thickness [21], and the early crustal thickening could have flowed away from the area of former Moho depression. This outflow would have also developed the anisotropy that is observed at present. Potentially complicating the simple interpretation of this crustal flow field, is the relatively recent emplacement of slow mantle beneath NE Wyoming associated with Yellowstone [31] and beneath central Colorado [32, 33]. In each case, the mantle is thought to be buoyant not only because it is seismically slow but also because each area has experienced substantial young uplift. This scenario might explain why crustal anisotropy is roughly radial around both of these mantle slow volumes (Figure 4.4B). California mantle anisotropy orientation south of the Juan de Fuca slab (JdF) rotates clockwise from 70-90◦ in eastern California to 100-130◦ in western California, near the San Andreas Fault (SAF). Compared to the mantle, the crustal anisotropy away from the Isabella mantle anomaly (ISA) and south of San Francisco has a similar azimuth in eastern California, but near the SAF the azimuth is further counter-clockwise at 110-150◦ [34]. Away from perturbations, anisotropy is expected to rotate progressively toward the direction of simple shear, i.e., it should be slightly counter-clockwise of the SAF. Assuming this behavior, there are two points to be made: (1) the crustal anisotropy is ahead of the mantle [35], suggesting either that the crust is driving the mantle or the crustal deformation zone is narrower (and more highly strained) than that of the mantle (Figure 4.4C), and (2) the anisotropy near the ISA mantle anomaly appears to be a simple "sinker" perturbation to the SAF-related anisotropy (Figure 4.4A). In northern California, the crustal anisotropy is coherent and organized around the southern JdF slab. The interpretation of these results is still not to- 62 Figure 4.4: Proposed mechanisms that give rise to the observed crustal anisotropy. (A) Negatively buoyant mantle anomalies. The load of the mantle lithosphere creates vertical stresses on the Moho and pull the crust down, creating a lateral pressure gradient that drives Poiseuille flow in the ductile mid-lower crust towards the mantle anomaly. The green arrow shows the sign of the lithospheric load. (B) Positively buoyant mantle anomalies. The hot and buoyant mantle rises and pushes the crust upward, creating a lateral pressure gradient that drives Poiseuille flow in the ductile mid-lower crust outwards from the mantle anomaly. The green arrow shows the sign of the lithospheric load. (C) Plate interactions. For the case of a transform fault, the deformation in the crust and mantle lithosphere is dominated by strain in a narrow shear zone, whereas in the asthenosphere deformation is in a zone whose width increases with depth. The green arrows show velocity relative to the fault, so that velocity is zero beneath the fault. Motion is in opposite directions on either side of the fault, and far from the fault the velocity is at plate rate. (D) No crustal anisotropy. There are no stresses from plate boundaries and no mantle-buoyancy anomalies. The mantle lithospheric strength isolates the crust from the sub-horizontal asthenospheric flow (i.e. mantle flow does not drive Couette flow in the mid-lower crust). tally clear, but the strong and complex crustal dynamics of the Mendocino triple junction area has been recognized and discussed by Furlong and Govers [36] and Liu et al. [37]. Between the southern JdF slab and the ISA mantle anomaly, the margin-normal crustal anisotropy is nearly E-W in orientation. A straightforward interpretation is crustal flow away from the high-standing Great Basin, much like that observed beneath the Tibetan region [38]. This flow occurs largely between the two seismically fast mantle features, perhaps owing to cooler lower crustal temperature and higher crustal viscosity above 63 these cooler mantle structures [which comes up to the Moho; 39]. 4.5 Discussion and Conclusions Including the observations of Castellanos et al. [1], we have three areas that provide different tectonic settings for study. Yet their crustal anisotropy fields can be explained with a few simple kinematic processes. From these observations we can make some general conclusions about continental crust during an orogeny. In the PNW and Rocky Mountains, crustal azimuthal anisotropy is roughly perpendicular to central seismically fast upper mantle structures (Castellanos et al., 2020; Figure 4.3A). The anisotropy is attributed to crustal flow that is approximately perpendicular to the location of seismically fast mantle. In each case there is no obvious topographic expression (the Wallowa Mountains are about 60 km south of the area of thick crust). For the Wallowa mantle anomaly, excitation is young and ongoing; in the Rocky Mountains, the anisotropy is Laramide in age and it appears to have preserved this last major strain event. For the crustal flow, we cannot distinguish between the inflow and outflow scenarios. Nonetheless, crustal flow appears to be poloidal in response to local buoyancy structure, i.e., the inferred flow is curl-free, centered on the buoyancy anomaly. Apparent outflow from the elevated Great Basin is similar, but it is associated with a strong topographic gradient. In contrast to the above examples, southern California crustal anisotropy is sub-parallel to the well-developed San Andreas shear zone and the crustal flow appears to be approximately toroidal, i.e. divergence-free, which is expected for a passive transform margin. The obvious anisotropy orientation perturbation near the ISA mantle anomaly is easily attributed to the superimposed poloidal flow, with crust converging on the site of mantle loading created by the ISA lithospheric load (Figure 4.3C). With respect to continental orogenies, we infer a few important consequences: (i) With mantle loading, isostatic balance is fundamentally maintained between the mantle load and compensating change in crustal thickness accommodated by Moho warping. However, if loading occurs more rapidly than crustal flow can supply compensating crust, a surface deflection will occur, such as the 80 Ma creation of basins in mechanism. With continued crustal flow and Moho adjustment, isostasy is ultimately maintained without a sur- 64 face expression and could be topographically invisible. This mechanism differs from tectonically-driven changes in crustal thickness, in which crustal thickness variations are associated with (and compensated by) topographic variations. (ii) Because the density difference across the Moho is small, modest mantle loads can create significant Moho warping and drive large amounts of crustal flow. This mechanism may be common in orogenies, and it is an effective way for straining lower crust to distances relatively far from the site of loading. In contrast, simple shear deformation (e.g., below the SAF or an area of detachment faulting) can accommodate large strain within the small volume of a shear zone, and anisotropy may be developed within only a relatively small volume of rock. (iii) The observed common occurrence of crustal anisotropy related to mantle loading suggests a widely distributed low-viscosity lower crust. Such crustal conditions may require an orogeny to supply crustal hydration or heating. Low crustal viscosity will tend to mechanically decouple the crust from horizontal mantle flow, in which case horizontal mantle flow would not play an important role in crustal tectonics. References [1] Jorge C Castellanos, Jonathan Perry-Houts, Robert W Clayton, YoungHee Kim, A Christian Stanciu, Bart Niday, and Eugene Humphreys. Seismic anisotropy reveals crustal flow driven by mantle vertical loading in the pacific nw. Science advances, 6(28):eabb0476, 2020. [2] Arthur H Lachenbruch, John Harvey Sass, and SP Galanis Jr. Heat flow in southernmost california and the origin of the salton trough. Journal of Geophysical Research: Solid Earth, 90(B8):6709–6736, 1985. [3] Arthur H Lachenbruch and Paul Morgan. Continental extension, magmatism and elevation; formal relations and rules of thumb. Tectonophysics, 174(1-2):39–62, 1990. [4] Bradford H Hager, Robert W Clayton, Mark A Richards, Robert P Comer, and Adam M Dziewonski. Lower mantle heterogeneity, dynamic topography and the geoid. Nature, 313(6003):541–545, 1985. [5] Bernhard Steinberger, Harro Schmeling, and Gabriele Marquart. Largescale lithospheric stress field and topography induced by global mantle circulation. Earth and Planetary Science Letters, 186(1):75–91, 2001. [6] Shimon Wdowinski and Gary J Axen. Isostatic rebound due to tectonic denudation: A viscous flow model of a layered lithosphere. Tectonics, 11 (2):303–315, 1992. 65 [7] Peter S Kaufman and Leigh H Royden. Lower crustal flow in an extensional setting: Constraints from the halloran hills region, eastern mojave desert, california. Journal of Geophysical Research: Solid Earth, 99(B8): 15723–15739, 1994. [8] Derek L Schutt, Anthony R Lowry, and Janine S Buehler. Moho temperature and mobility of lower crust in the western united states. Geology, 46(3):219–222, 2018. [9] Magdala Tesauro, Mikhail K Kaban, Walter D Mooney, and Sierd APL Cloetingh. Density, temperature, and composition of the n orth a merican lithosphere—new insights from a joint analysis of seismic, gravity, and mineral physics data: 2. thermal and compositional model of the upper mantle. Geochemistry, Geophysics, Geosystems, 15(12):4808–4830, 2014. [10] Lisa Block and Leigh H Royden. Core complex geometries and regional scale flow in the lower crust. Tectonics, 9(4):557–567, 1990. [11] Peter Bird. Lateral extrusion of lower crust from under high topography in the isostatic limit. Journal of Geophysical Research: Solid Earth, 96 (B6):10275–10286, 1991. [12] Peter Bird, Zhen Liu, and William Kurt Rucker. Stresses that drive the plates from below: Definitions, computational path, model optimization, and error analysis. Journal of Geophysical Research: Solid Earth, 113 (B11), 2008. [13] Carolina Lithgow-Bertelloni and Jerome H Guynn. Origin of the lithospheric stress field. Journal of Geophysical Research: Solid Earth, 109 (B1), 2004. [14] Vladislav Babuska and Michel Cara. Seismic anisotropy in the Earth, volume 10. Springer Science & Business Media, 1991. [15] MK Savage. Seismic anisotropy and mantle deformation: what have we learned from shear wave splitting? Reviews of Geophysics, 37(1):65–106, 1999. [16] Fan-Chi Lin, Michael H Ritzwoller, Yingjie Yang, Morgan P Moschetti, and Matthew J Fouch. Complex and variable crustal and uppermost mantle seismic anisotropy in the western united states. Nature Geoscience, 4(1):55–61, 2011. [17] Haiying Gao, Eugene D Humphreys, Huajian Yao, and Robert D van der Hilst. Crust and lithosphere structure of the northwestern us with ambient noise tomography: Terrane accretion and cascade arc development. Earth and Planetary Science Letters, 304(1-2):202–211, 2011. 66 [18] Haiying Gao. Crustal seismic structure beneath the source area of the columbia river flood basalt: Bifurcation of the moho driven by lithosphere delamination. Geophysical Research Letters, 42(22):9764–9771, 2015. [19] A Christian Stanciu and Eugene D Humphreys. Upper mantle tomography beneath the pacific northwest interior. Earth and Planetary Science Letters, 539:116214, 2020. [20] Brandon Schmandt and Eugene Humphreys. Complex subduction and small-scale convection revealed by body-wave tomography of the western united states upper mantle. Earth and Planetary Science Letters, 297 (3-4):435–445, 2010. [21] Weisen Shen and Michael H Ritzwoller. Crustal and uppermost mantle structure beneath the united states. Journal of Geophysical Research: Solid Earth, 121(6):4306–4342, 2016. [22] Brandon Schmandt, Fan-Chi Lin, and Karl E Karlstrom. Distinct crustal isostasy trends east and west of the rocky mountain front. Geophysical Research Letters, 42(23):10–290, 2015. [23] William BF Ryan, Suzanne M Carbotte, Justin O Coplan, Suzanne O’Hara, Andrew Melkonian, Robert Arko, Rose Anne Weissel, Vicki Ferrini, Andrew Goodwillie, Frank Nitsche, et al. Global multi-resolution topography synthesis. Geochemistry, Geophysics, Geosystems, 10(3), 2009. [24] Jorge C Castellanos and Robert W Clayton. The fine-scale structure of long beach, california, and its impact on ground motion acceleration. Journal of Geophysical Research: Solid Earth, In-Review. [25] Florent Brenguier, Philippe Kowalski, N Ackerley, N Nakata, P Boué, M Campillo, E Larose, S Rambaud, C Pequegnat, T Lecocq, et al. Toward 4d noise-based seismic probing of volcanoes: Perspectives from a largen experiment on piton de la fournaise volcano. Seismological Research Letters, 87(1):15–25, 2015. [26] Nori Nakata, Pierre Boué, Florent Brenguier, Philippe Roux, Valérie Ferrazzini, and Michel Campillo. Body and surface wave reconstruction from seismic noise correlations between arrays at piton de la fournaise volcano. Geophysical Research Letters, 43(3):1047–1054, 2016. [27] Martin L Smith and FA Dahlen. The azimuthal dependence of love and rayleigh wave propagation in a slightly anisotropic medium. Journal of Geophysical Research, 78(17):3321–3333, 1973. [28] Th W Becker, S Lebedev, and MD Long. On the relationship between azimuthal anisotropy from shear wave splitting and surface wave tomography. Journal of Geophysical Research: Solid Earth, 117(B1), 2012. 67 [29] Shaofeng Liu, Dag Nummedal, and Lijun Liu. Migration of dynamic subsidence across the late cretaceous united states western interior basin in response to farallon plate subduction. Geology, 39(6):555–558, 2011. [30] Zhiyang Li and Jennifer Aschoff. Constraining the effects of dynamic topography on the development of late cretaceous cordilleran foreland basin, western united states. GSA Bulletin, 2021. [31] Kenneth L Pierce and Lisa A Morgan. Is the track of the yellowstone hotspot driven by a deep mantle plume?—review of volcanism, faulting, and uplift in light of new data. Journal of Volcanology and Geothermal Research, 188(1-3):1–25, 2009. [32] Andres Aslan, Karl E Karlstrom, Laura J Crossey, Shari Kelley, Rex Cole, Greg Lazear, Andy Darling, LA Morgan, and SL Quane. Late cenozoic evolution of the colorado rockies: Evidence for neogene uplift and drainage integration. Through the generations: Geologic and anthropogenic field excursions in the Rocky Mountains from modern to ancient: Geological Society of America Field Guide, 18:21–54, 2010. [33] Greg Lazear, Karl Karlstrom, Andres Aslan, and Shari Kelley. Denudation and flexural isostatic response of the colorado plateau and southern rocky mountains region since 10 ma. Geosphere, 9(4):792–814, 2013. [34] Minoo Kosarian, Paul M Davis, Toshiro Tanimoto, and Robert W Clayton. The relationship between upper mantle anisotropic structures beneath california, transpression, and absolute plate motions. Journal of Geophysical Research: Solid Earth, 116(B8), 2011. [35] Sylvain Barbot. Mantle flow distribution beneath the california margin. Nature communications, 11(1):1–14, 2020. [36] Kevin P Furlong and Rob Govers. Ephemeral crustal thickening at a triple junction: The mendocino crustal conveyor. Geology, 27(2):127–130, 1999. [37] Kaijian Liu, Alan Levander, Yongbo Zhai, Robert W Porritt, and Richard M Allen. Asthenospheric flow and lithospheric evolution near the mendocino triple junction. Earth and Planetary Science Letters, 323: 60–71, 2012. [38] Jiayi Xie, Michael H Ritzwoller, Weisen Shen, and Weitao Wang. Crustal anisotropy across eastern tibet and surroundings modeled as a depthdependent tilted hexagonally symmetric medium. Geophysical Journal International, 209(1):466–491, 2017. [39] A Christian Stanciu and Eugene Humphreys. An improved tomographic model for southern california upper mantle seismic structure. In AGU Fall Meeting 2019. AGU, 2019. 68 Chapter 5 THE FINE-SCALE STRUCTURE AT LONG BEACH, CALIFORNIA: SURFACE WAVES. This chapter was adapted from: J. C. Castellanos, and R. W. Clayton (In review). “The fine-scale structure of Long Beach, California, and its impact on ground motion acceleration.” In: Journal of Geophysical Research. 5.1 Abstract The metropolitan Los Angeles region represents a zone of high-seismic risk due to its proximity to several fault systems, including the San Andreas fault. Adding to this problem is the fact that Los Angeles and its surrounding cities are built on top of soft sediments that tend to trap and amplify seismic waves generated by earthquakes. In this study, we use three dense petroleum industry surveys deployed in a 16x16-km area at Long Beach, California, to produce a high-resolution model of the top kilometer of the crust and investigate the influence of its structural variations on the amplification of seismic waves. Our velocity estimates reveal substantial lateral contrasts and correlate remarkably well with the geological background of the area, illuminating features such as the Newport-Inglewood Fault, the Silverado Aquifer, and the San Gabriel River. We then use computational modeling to show that the presence of these small-scale structures have a clear impact on the intensity of the expected shaking, and can cause ground-motion acceleration to change by several factors over a sub-kilometer horizontal scale. These results shed light onto the scale of variations that can be expected in this type of tectonic setting and highlight the importance of resolution in modern-day seismic hazard estimates. 5.2 Introduction Los Angeles County and its surrounding areas are known for their susceptibility to earthquake shaking and for having several crossing faults that are capable of producing major damage to their infrastructure. The 1933 Mw 6.4 Long Beach earthquake, for instance, caused widespread damage throughout the southern part of Los Angeles Basin, resulting in more than 100 civilian fatalities [1] 69 and over 40 million dollars in economic loss [2]. The 1994 Mw 6.7 Northridge earthquake caused several gas leaks and collapsed numerous roadways and buildings across the region, making it one of the most costly disasters in US history [3]. More recently, the 2019 Ridgecrest earthquake sequence struck the Californian Mojave Desert, causing significant structural damage near its epicenters and noticeable shaking in the greater Los Angeles area. Modern seismic hazard assessments are generally performed using deterministic approaches, in which a realistic representation of both the source process and the subsurface properties are used to obtain an estimate of the ground motion due to a large earthquake [e.g., 4–6]. As our knowledge of the crustal velocity distribution across southern California has evolved, these calculations have made it possible to identify zones in which the heaviest damage can be expected and have ultimately allowed the generation of building codes to make structures more resistant to earthquakes [7–9]. However, with the ever increasing amount of seismic records, it has become clear that there are still significant differences between modeled and observed ground-motions across multiple scales that tend to produce large errors in seismic hazard estimates. Such discrepancies appear to stem from the fact that the shallow structure in Los Angeles Basin is not well-characterized by the 3D models used in the simulations and that, as a result, its influence on the propagation of seismic waves is not entirely captured [10, 11]. These shortcomings have important societal impacts as the underprediction of expected ground motion possess an immediate risk to the inhabitants of the area, whereas the overprediction leads to more costly structural design and construction [12]. Substantial work has been done in order to improve the accuracy of the velocity models of southern California. Examples of such efforts include the development of sophisticated imaging methods [e.g., 13–15], the incorporation of more seismic data [e.g., 16–18], and the deployment of temporary nodal arrays that allows us to explore the complex architecture of Los Angeles Basin with refined detail [e.g., 19]. In particular, the introduction of oil-industry surveys into the field of crustal geophysics has dramatically increased the resolution of regional velocity models and, with it, improved the prediction of several seismic observables [20, 21]. Because of their unprecedented density, these seismic arrays are perhaps the most reliable way to determine the finescale structure of the crust and investigate the many complicated phenomena 70 that can occur during seismic wave propagation. However, because of their financial cost, only a few regions around the world have hosted this type of experiments, resulting in a lack of knowledge on both the scale of structural variations and the spatial variability in shaking intensity that we can expect in this type of tectonic settings. In this study, we use three of the highest resolution seismic arrays ever deployed in southern California to probe the elastic properties of the Long Beach crust and investigate the seismic response of the ground beneath the surveys. Within this scope, we present a newly developed framework that allows us to perform reliable surface wave phase velocity measurements and construct a high-resolution 3-D model of the crust. We then report the peak ground motion acceleration that is expected at different sections of the experiment as obtained by propagating numerous wavefields through our velocity estimates. Finally, we discuss the potential of these type of arrays as well as their impact on future hazard assessment studies. 5.3 Data and Methods Ambient Noise Data Three dense arrays, with a total of ∼13,000 seismic stations, were deployed as part of different petroleum industry surveys in Long Beach, CA, across the northwest-southeast trending Newport-Inglewood Fault (NIF) system (Figure 5.1). Each of these networks (commonly known as the Long Beach, the Extended Long Beach and the Seal Beach arrays) consisted of a large group of 100-meter-equispaced high-frequency velocity sensors that were designed to illuminate the shallow oil deposits associated with the faulting of the area [22]. Fortunately, during their time of operation, these instruments recorded continuously and, therefore, not only captured the active source component of the surveys, but also several passive sources including the ambient noise field. This last characteristic of the experiments allows us to move beyond the focus of traditional petroleum seismology and use this type of instrumentation to investigate local microseismicity [e.g., 23], the mechanics of active fault zones [e.g., 24], and the structure of the deeper crust [e.g., 20, 25–27]. In recent years, ambient-noise tomography has become a well-established imaging tool for investigating the elastic properties of the subsurface. This technique starts by cross-correlating the Earth’s background vibrations recorded 71 Figure 5.1: Regional map of southern California and the Continental Borderland. The black rectangle marks the area where the three petroleum industry surveys were deployed and the inset map shows the distribution of the instruments (with the Long Beach array in red, the Extended Long Beach in green, and the Seal Beach array in blue). The stations of the Southern California Seismic Network (SCSN) are shown as inverted triangles and the regional mapped faults are delineated with thick gray lines. Seismic station LAF is shown as an inverted red triangle. at a pair of receivers to approximate the Green’s function between the two instruments [28–30]. Once a sufficiently long segment of ambient noise is cross-correlated, and the correlation function converges, classical tomographic techniques are applied to the emerged signals to extract critical information of the average velocity structure of the Earth along the path that waves have traveled [31–34]. Because the resolution of ambient-noise tomography depends only on the geometry and distribution of seismic instruments, the analysis of noise-derived ballistic waves has enabled the construction of velocity models with unprecedented detail, and even allowed the mapping of the near-surface structure of seismically quiescent areas [e.g., 35]. Here, we expand on the work of Lin et al. [20] and use ambient noise data 72 from all three temporary oil industry surveys to extract the Green’s function between nodal instruments. Because these arrays operated at different times and, therefore, recorded asynchronously, each network was only crosscorrelated with itself. The processing scheme used to process the noise recordings closely follows that of Bensen et al. [30] and resulted in more than 30 million correlograms that show clear surface waves traveling between the different pairs of stations (Figure 5.2). Since the dense arrays were composed of vertical geophones only, the emerged signals likely correspond to Rayleigh waves. These arrivals are the base of our investigation and are the tool with which we image the structure beneath the arrays. Figure 5.2: Wavefield snapshots of the ambient noise correlation functions from a single virtual source of the Long Beach array (left), Extended Long Beach array (middle), and the Seal Beach array (right). Clear ballistic waves propagating away from all three of the sources can be observed. The number of instruments and time of operation of each survey is presented on top of each panel. Inversion for a 1-D Velocity Model As with every structural investigation, it is desirable have some a-priori knowledge of the average velocity structure of the target area before characterizing the small deviations from it. To this end, we use all of the available crosscorrelations to generate a stacked gather for each of the surveys. We then convert these record sections to the T au-p domain via a slant-stack scheme to obtain the frequency-phase-velocity representation of the signals [36]. The collective analysis of multiple-offset traces allows us to retrieve clear fundamental mode (FM) and first overtone (FO) dispersion curves for the 0.5-4 Hz frequency band, which we use to construct a shear-wave velocity profile that 73 represents the mean structure of the top kilometer of the crust beneath the arrays. For the actual inversion, we average all of the extracted phase velocities to build a single FM and FO dispersion curve, and use Haney and Tsai [37] perturbational algorithm to jointly invert the two modes for the best-fitting shear-wave velocity model. The average correlation gather, dispersion curves, modes sensitivities, and inverted velocity profile are shown in Figure 5.3a-e. Figure 5.3: Stacked gather of the ambient noise cross-correlations and inversion for an average 1-D shear-wave velocity model. (a) Shows the first 30 seconds of the stacked cross-correlations from all three surveys. The bin-size of the spatial stacking is of 100 meters. The waveforms are band-passed filtered between 0.5 and 4 Hz. Clear fundamental (FM) and first overtone (FO) surface waves traveling at velocities between 0.3 and 1.5 km/s can be observed for the entire distance range. (b) Shows the average FM and FO dispersion curves that are obtained through the slant stack analysis of each survey’s dataset. The continuous red lines represent the predicted dispersion curve that results from our inversion process. (c-d) Shows the sensitivity kernels for the FM and FO, respectively. (e) Shows a comparison between shear wave velocity profile that is obtained from jointly inverting dispersion curves in (b) (black line) and a profile that is extracted from the CVM-S4 model [38] at the geographic center of the nodal stations (red line). The black dashed line on top of panels (c-e) represents the free surface (FS). Surface Wave Relative Traveltime Measurements Although the slant-stack analysis of entire record sections allows us to retrieve robust velocity estimates, as evidenced by the agreement between our inversion 74 results with the CVM-S4 reference model (Figure 5.3e), our ultimate goal is to analyze the small velocity variations across the seismic arrays. Such a task, however, relies on our capability to extract reliable velocity measurements from individual cross-correlation functions, which can be problematic due to scattering, attenuation, and the inherent high-noise level of the signals. For this reason, we implement an automatic neighborhood-based cross-correlation method for phase arrival picking that allows us to deal with complicated waveforms [39, 40]. This technique begins by band-pass filtering all of the existing traces of a given source to a particular narrow-band frequency and measuring their differential delay times, ∆t, with respect to the waveforms recorded at neighboring stations (Figure 5.4). To perform this calculation, we window the correlograms around the expected arrival time of the ballistic wave using a cosine-tapered window that is flat for 5 times the center frequency of the band-pass filter. Moreover, to address the problem of cycle-skipping, we limit the radius of the interstation cross-correlation to one wavelength long, as determined by the reference dispersion curve that was extracted from the stacked correlation gather in Figure 5.3. Once all of the cross-correlation derived differential delay times are collected, we arrange them into the form:  1  1   1  .  ..   1   0  . . .  1 −1 0 · · · 0 −1 · · · 0 0 ··· .. .. . . . . . 0 0 ··· 1 −1 · · · .. .. . . . . . 1 1 ···     tR1  ∆t 0    R1 R2      tR2     ∆tR1 R3    0        tR3     ∆tR1 R4   0      tR4   ..  ..      .  .     ,   tR5  =       −1   ∆tR1 RN    tR6      ∆t   0   R2 R3   t     ..   R7   ..  .   .   .     ..     0 1 tRN  (5.1) so that the relative arrival time, t, of the wavefront recorded by every RN receivers can be inverted using a simple least squares method [41]. To ensure that the traveltime field that is obtained from this operation is locally smooth, we introduce the regularization term: ! |r − r’| , S(r, r’) = exp − 2σ 2 (5.2) 75 Z S S(r, r’)dr’ = 1, (5.3) where r is the position vector of the stations and σ is the spatial smoothing width or correlation length [42]. Here, we set the value of this term to a half-wavelength of the filter’s center frequency, and determine the overall strength of the smoothing kernel, S, from an L-curve analysis (Figure 5.5a). For the actual inversion, we chose to consider only high-quality waveforms and remove all relative delay times with a correlation coefficient smaller than 0.90. Figure 5.5b-c shows the inverted relative traveltime measurements for the 1-s period fundamental model Rayleigh wave that is derived from a single virtual source located at the northwestern end of the Long Beach Array. Note that, to honor the far-field approximation, we only solve for the relative traveltimes of stations that are more than three wavelengths away from the virtual source [31, 33]. To produce 3-D velocity model of the subsurface, it is necessary to characterize the propagation properties of surface waves that are sampling different depths of the crust. This task, in turn, requires us to extend our analysis to a wide range of frequencies, which can be challenging particularly when moving to a higher spectrum. For this reason, we introduce an intermediate step into our processing scheme, where we use the differential delay times of the long period waveforms to track the correct phase delays to higher frequencies (Figure 5.5d). This implementation allows us to use simpler signals to guide the relative traveltime surfaces to higher frequencies and hence avoid cycle-skipping. In this study, we apply this algorithm to every virtual source available in all three of the temporary surveys to generate relative traveltime surfaces of the fundamental mode Rayleigh waves for the 0.5-4 Hz frequency band. For these calculations, we used the 1-s measurements to track the correct phase delay to higher frequencies, and performed an individual inversion for all the longer periods. We also disregarded the first-overtone surface waves since these wave packets were not always visible in individual waveforms and, for certain offsets and frequencies, the two modes interfered with each other, with the fundamental mode dominating the first overtone. Lastly, we remove all data involving the 300 marine nodes of the Seal Beach array as these recordings were particularly noisy and no clear surface waves from these virtual sources were observed. This operation resulted in almost 13,000 relative traveltime surfaces 76 Figure 5.4: Schematic representation of the neighborhood-based crosscorrelation method for phase arrival picking. Receivers that are less than one wavelength apart (<1λ) from a reference station (R1 ) are grouped together to form a subarray. The differential delay times, ∆t, between the reference waveform and all the other waveforms in the subarray are measured through waveform cross-correlation. for each frequency that was analyzed, which we then used to construct phase velocity maps. Surface Wave Phase Velocity Measurements We use the collection of the relative traveltimes, to invert for the frequencydependent phase velocities. To do so, we implement a time-based beamforming scheme, where, for a given virtual source, we take groups of relative arrival times that are less than a half-wavelength away from a given reference receiver and use a plane-wave approximation to solve for the horizontal slowness vector, u, that can best explain the observations (Figure 5.6a-b). Here, we define said vector as: 77 Figure 5.5: Example of relative phase delay measurements for a virtual source of the Long Beach array. (a) Shows the L-curve analysis that was performed to determine the strength of the regularization term, S, in the inversion of the relative traveltimes. (b) Shows how the resulting relative traveltimes of the 1-s period waveforms vary as a function of distance from the virtual source. The dashed red line on top of the measurements has a corresponding slope of 0.70 km/s. (c) Shows the spatial distribution of the relative traveltimes that are in (b). The contours are drawn in time increments of 1 second. The red star denotes the location of the virtual source and the thick dashed circle around the source marks the three-wavelength minimum distance criterion that was imposed in our analysis. (d) Shows the process by which we extend the differential delay measurements to higher frequencies. The waveforms on the top are two ambient noise cross-correlations functions that are filtered at 1-s period and windowed around the expected time of arrival of the ballistic surface waves. The source-receiver raypath of these two waveforms is shown in (c). The bottom panel shows the correlation function of the two traces as a function of period. To determine the differential delay time of these two waveforms across the entire frequency range, we use the 1-s measurement to track down the phase delay to higher frequencies (red dashed line).     cos(θ)  ux  u = u = , sin(θ) uy (5.4) 78 where u represents the horizontal slowness and θ the direction of propagation of the wave. From this expression, we can then write the relative arrival time of a group of receivers as: t(x, u) = u · (x − xc ), (5.5) where the vector x represents the spatial coordinates of every station in the subarray and xc the center coordinates of the subarray. This formulation allows us to directly solve for the local phase velocity and the wave’s propagation direction at the different sections of the array with very few assumptions. However, the primary reason why we measure phase velocities with this approach is because this operation is efficient enough to allow the use of Gaussian statistics to quantify the robustness of the inverted slownesses and azimuths. To do so, we apply a bootstrapping method to every subarray relative traveltimes using a total of 300 resamples and, for each family of measurements, obtain a new local phase velocity and azimuth (Figure 5.6c-d). We then translate the spread of these quantities into an absolute error value by taking the product of the eigenvalues of their covariance matrix. To show the performance of our scheme, Figure 5.6e-g shows the local phase velocity, backazimuth, and associated error for the 1-s period Rayleigh wave that is derived from a single virtual source located at the southeastern end of the Long Beach Array. It is worth noticing how the stations with the largest error are coincident with the sections of the survey that have a relatively anomalous phase velocity and backazimuth estimation. Once all of the virtual sources have been used to derive the frequency-dependent phase velocity measurements, we weight-stack them to produce an isotropic phase velocity map. For this process, we use the bootstrap-derived errors as a means to weigh the phase velocity values and opt to only retain sites in which more than 100 velocity measurements were made. This step allows us to dramatically reduce the effects of irregular measurements and obtain a robust representation of the phase velocity at each site. Figure 5.7 schematically shows the construction process of the 1-s period isotropic phase velocity map of the Long Beach Array, together with a group of higher-frequency velocity maps that were generated using the longer period relative phase delay times as a reference. 79 Figure 5.6: Example of phase velocity measurements for a virtual source of the Long Beach array at 1-s period. (a) Shows the source and receiver configuration that is used to make a single phase velocity measurement at a particular site. The yellow star marks the reference station that is used to form the receiver subarray and the green circles correspond to all the stations that were grouped together using the half-a-wavelength maximum distance criterion. The aperture of the subarray is approximately one-wavelength-long. (b) Shows the relative traveltime of the receiver subarray in (a) as a function of their distance from the virtual source. Each marker is color-coded by the amount of times that a differential measurement of that station was used to construct the relative traveltime surface. The red circles depict the predicted relative arrival times that result from our time-based beamforming. (c-d) Show the distribution of our bootstrapped backazimuth and velocity measurements for the data shown in (b). (e-g) Shows the spatial distribution of all the velocity, backazimuth, and error measurements that were obtained for the same virtual source. Up to this point, all of the processing steps and examples that have been presented here have been derived using data from the Long Beach survey. Nonetheless, the adaptability of our framework makes it straightforward to extend our analysis to the neighboring arrays. Figure 5.8 shows the phase velocity distribution for the 1-s period Rayleigh wave that is obtained by com- 80 Figure 5.7: Schematic representation of the construction process of isotropic phase velocity maps. (a-c) Show the phase velocities, backazimuths, and error measurements for different virtual sources of the Long Beach array at 1-s period. The errors in (c) are used to weight the averaging of the velocities in (a) to ultimately construct a frequency-dependent isotropic velocity map. The backazimuth measurements in (b) are reserved and used at a later step to characterize the anisotropy across the array. (d) Shows, from left to right, the isotropic phase velocity maps for the 1-, 0.8-, 0.6-, and 0.40-s Rayleigh waves. Each of these maps are generated in the same fashion, with the exception that the relative phase delay measurements of the higher frequencies were derived from the 1-s period waveforms. For periods larger than 1-s, such measurements were done independently. bining all three of the temporary surveys. From this image, we can note the continuity of the phase velocity estimates, as evidenced by the northwestsoutheast trending fast anomaly that is characteristic of the NIF [20]. This feature of our velocity map suggests consistency in our processing scheme. Moreover, we compare our phase velocity measurements across the Seal Beach array with a reflective time slice that was provided by the oil company that processed the active source data [43]. From this comparison, we find that there is a remarkable agreement between our velocity estimates and the main reflection horizons, particularly those associated to the Garden Grove Fault 81 and the prograding clinoforms located near the center of the survey. Figure 5.8: Comparison of our 1-s isotropic phase velocity map (left) with a reflective time slice of the active part of the survey at 0.40 s (right). Despite their difference in nature, the two different quantities in these images appear to illuminate similar tectonic features. The reflection image is courtesy of 3D Seismic Solutions. Inversion for Shear Wave Velocity After all of the velocity maps are constructed, we extract a regionalized dispersion curve at each point across the survey and use the Haney and Tsai [37] perturbational algorithm to map the phase velocity as a function of frequency to shear wave velocity as a function of depth. For each inversion, we use a homogeneous velocity profile that is discretized into 25-m layers to iteratively update each component of the model until a satisfactory fit to the dispersion curve is achieved (Figure 5.9). Throughout this process, we assume a constant Vp /Vs of 1.8 for the entire structure and determine the density from the compressional wave velocity using the empirical relationship of Gardner et al. [44]. Each velocity profile is then integrated into a whole 3-D velocity model that represents the structure of the top kilometer of the Long Beach crust. 5.4 Results and Discussion Figure 5.10a shows three depth-slices of the shear wave velocity model that is obtained through the inversion of the regionalized dispersion curves. It is clear from these images, that there are significant variations in the velocity 82 Figure 5.9: Inversion process of the regionalized dispersion curves. The maps on the left show how the phase velocities vary as a function of period. The black dashed black lines on the right panels correspond to the dispersion curves that are extracted at locations A and B of the maps. The red continuous lines correspond to the dispersion’s best-fit. The number of updates that were required for the inversion to converge in these two examples is shown on the top-right of each panel. structure beneath the surveys. At shallow depths (<100 m), the range of velocity variations across the survey is close to 10%, with the slower velocities being concentrated in the southern part of the Long Beach Array and around the San Gabriel River. In these images the Silverado water bearing unit and a section of the Compton-Los Alamitos Fault system can be seen. Both of these structures are of great importance of the area as the former is the main supplier of the groundwater that is extracted in Long Beach [45], and the latter is part of a blind thrust fault that is known to be capable of generating large magnitude earthquakes (Mw 7.0-7.4) [46]. At deeper depths (∼600 m), the most prominent feature in our velocity maps is the NIF. In agreement with Lin et al. [20], we argue that the fast velocity signature of this structure is related to the presence of deeper (and hence faster) rocks that were exhumed by the transpressional stresses that are acting on this region [22, 47]. This deformation process would also explain why the migrated reflection images generated by the oil company that carried out the seismic survey delineate the NIF as a flower-shaped damaged-zone rather than a single 83 Figure 5.10: Depth slices and cross-sections of our velocity model. (a) Shows a depth slice our velocity model at 50-, 100-, and 600-m depths. (b) Shows four cross-sections of our velocity model plotted on top of the migrated reflection images that were generated by the petroleum company that carried out the Seal Beach survey. To highlight the amount of lateral variation in each of our profiles, we removed the average velocity at each depth. The black dashed lines are used to mark prominent discontinuities in the migrated images. The location of each cross-section is shown in the left-most map in (a). The migrated reflection profiles are courtesy of 3D Seismic Solutions. planar boundary that extends continuously to the deep crust (Figure 5.10b). The identification of this type of structures is relevant not only because it holds the potential of forming profitable hydrocarbon traps [48], but also because it can provide important insights into the type and extent of faulting that we can expect in a particular area [24, 49]. As a final product of the surface-wave analysis, we present an azimuthal anisotropy model for the structure beneath the surveys. This quantity describes how the velocity of a particular wave varies with its direction of propa- 84 gation which, for the very shallow earth, is often caused by the preferred alignment of cracks by the regional tectonic forces [50]. The analysis of anisotropy can thus illuminate the direction of the stress field and potentially detect faults that might be present in the area [e.g., 51, 52]. Figure 5.11 shows how the phase velocities for the 1-s Rayleigh wave vary as a function of azimuth at two different sites of the Long Beach array. These variations display a strong directional dependence that can be well-described by the 2θ components of Smith and Dahlen [53] model for surface wave phase velocity, C, in a weakly anisotropic media: C (T, θ) =C0 (T ) + C1 (T ) cos (2θ) + C2 (T ) sin(2θ)+ C3 (T ) cos (4θ) + C4 (T ) sin(4θ), (5.6) where T represents the period, θ the backazimuth, C0 the isotropic velocity, and C1−4 the azimuthal coefficients [54]. This parametrization allows us to translate the wavefield’s azimuthal dependence into a fast direction term and an amplitude term, which we can then visualize on a horizontal plane to determine whether any spatial patterns exists, and what association they have with the regional geologic structures. Figure 5.12 shows how the anisotropy for the 1-s Rayleigh waves varies across surveys. These measurements are mostly sensitive to the 100 meters depth and exhibit a dominant north-south orientation that is overall consistent with the maximum horizontal compressional direction of southern California [55]. However, with the resolution that is provided by the nodal arrays, we can now begin to illuminate sub-kilometer variations of the stress field, as is evident by the slight rotation of the anisotropy vectors just at the center of the NIF near Signal Hill (yellow star in Figure 5.12). This observation is consistent with the results of Lin et al. [20], who perform a similar analysis with the Long Beach array data and suggested that there might be a small-scale stress change near the fault zone segmentation. With the modest improvements in our processing scheme, which permits the retrieval of velocity estimates with much less spatial averaging, we are able to detect and interpret even smaller changes in the anisotropy pattern that could potentially be related to unmapped faults or other geological boundaries (red dashed lines in Figure 5.12). To validate this claim, Figure 5.13 shows snapshots of a propagating wavefront that is obtained by cross-correlating the ambient noise data of the dense surveys with the broadband station LAF 85 Figure 5.11: Example of the Rayleigh wave azimuthal dependence at 1-s period for two different locations at the Long Beach array. The color of the point cloud represents the density of the velocity measurements and the black curve marks the predicted directionality using the first three terms in Equation 5.6. The best fitting coefficients, along with their 95% confidence intervals, are shown in the bottom right of each panel. from the SCSN (see Figure 5.1 for its location). The path to this particular instrument is oriented parallel to most tectonic faults in the area (relative to the port of Long Beach) and is thus ideal to detect the presence of any geologic structures that has a similar strike. This figure shows that there exist large lateral variations on the propagation properties of the wavefield and that the wavefront appears to be "cut" just along the planes where our anisotropy measurements suggest the presence of a structural barrier (cyan dashed lines in Figure 5.13). This feature is largely evident at propagation time of 45 seconds, where the wavefront is deflected just at the boundary where our anisotropy vectors display a dramatic change in orientation, at the site where the seismic reflection image marks the location of the Garden Grove Fault (see reflection image in Figure 5.8). These results suggest that the analysis of azimuthal anisotropy has the potential of detecting sharp structural boundaries. With the robust velocity estimates that were derived from our ambient noise analysis, it is now possible to investigate the effects of the fine-scale structure 86 Figure 5.12: Anisotropy measurements for the 1-s period Rayleigh wave for the Long Beach array. The bars point to the fast direction of anisotropy and their length is proportional to the magnitude. The colored markers behind each anisotropy vector denote the isotropic velocity (C0 term in equation 5.6). The continuous red line on the right map marks the mapped surface trace of the regional faults, whereas the dashed red lines delineate sharp changes in the anisotropy measurements that we suggest might be indicative of the presence of faults or other geologic boundaries. We associate these structures to the following: (a) NIF zone; (b) Silverado aquifer; (c) Compton-Los Alamitos Fault; (d) Garden Grove Fault; (e) Prograding clinoforms (see reflection image in Figure 5.8). The yellow star marks the location of Signal Hill. on ground motion acceleration. For this endeavor, we propagate realistic 3-D wavefields through our high-resolution velocity model using Salvus, a highperformance spectral-element solver [56], and calculate the amount of shaking that is expected on different parts of the survey. The peak ground acceleration (PGA) that is obtained in this step is then compared to the one that results from propagating an equivalent wavefield through the 1-D model that was obtained from the inversion of the average dispersion curves in (Figure 5.3b). This approach, in principle, allows us to characterize the influence of the local geological structure on the amplitude of the seismic waves and separate the influence of factors such as the source radiation pattern and geometric spreading from the response. 87 Figure 5.13: Wavefield snapshots of the ambient noise correlation functions using station LAF of the SCSN as a virtual source. Clear fundamental mode (FM) and first overtone (FO) surface waves can be observed. The dashed lines are the same as in Figure 5.12 and delineate the zones where the anisotropy measurements exhibit a drastic change in direction. The wavefield was bandpassed filtered between 0.5 and 1.5 Hz. The time-stamp of each frame is shown in the upper right side of the panels. To provide a concrete example of the above, we place a vertical point source at the surface on the north-western end of the Long Beach survey, just a-top of the NIF, and record the entire time evolution of the wavefield (Figure 5.14). Here, because of the choice of source geometry and configuration, most of the energy that is inserted into the system will propagate in the form of surface waves [57]. This exercise allows us to see the largely asymmetric shape of the wavefront and reveals the existence of significant lateral variations on the magnitude of acceleration as the energy propagates away from the source. Such variations are particularly obvious in Figure 5.15, where we present the PGA ratio that is obtained by comparing our estimates with those obtained from the 1-D modeling. From this analysis, we can quantitatively observe how the small-scale variations in the shallow velocities can produce rapid changes in the vertical PGA and, most interestingly, how the structure of the NIF can generate elongated zones of intense amplification that result from the focusing and self-intersection of the propagating wavefront (i.e. from multipathing). These findings are in remarkable agreement with those of Clayton et al. [58], who used data from a nearby earthquake to document that the velocity contrasts underneath the Long Beach array can cause ground accelerations to vary by a factor of five over a horizontal scale length of one kilometer. The results in Figure 5.15 show that near-surface velocity variations can significantly alter the propagation properties of earthquake wavefields. To quantify 88 Figure 5.14: Snapshots of a synthetic wavefield propagated through our velocity model. The location of the source is indicated by the yellow star and corresponds to a vertical force injected by a ricker wavelet with a central frequency of 2 Hz. Note the lateral differences in the wavefront’s amplitude as it propagates away from the source. The time-stamp of each frame is shown in the lower right side of the panels. this phenomenom, we determine the amount of shaking that is expected from a vertically incident shear wave by placing an earthquake source beneath each station of the survey and propagating it to the surface [e.g., 21, 59, 60]. For a source, we use a horizontal force injected by a 2 Hz Ricker wavelet located at 2-km depth (Figure 5.16). We then compare these estimates with the ones that are obtained from an equivalent 1-D run to obtain a measure of the relative amplification that occurs at each site. We record the three-component acceleration and calculate the PGA from the magnitude of acceleration to account for any non-vertical propagation that might occur due to the heterogeneous 3-D structure. The amplification measurements that are obtained from this exercise are presented in Figure 5.17. The amplitude of variation in the amplification factors that we obtain through this analysis is smaller than the one we observe from a distant source. However, it is important to note that these measurements are meant to look at the amount of local amplification at each individual site without taking into con- 89 Figure 5.15: Map of relative amplification for the source configuration shown in Figure 5.14. These values describe the expected PGA relative to the one that is obtained from a 1-D model. The red star marks the location of the source. The right panels show the synthetic waveforms that are recorded at locations A and B of the map for both our 3-D model (black) and a 1-D model (red). Note how the structure of the NIF generates elongated zones of intense amplification that result in ground motion changes that vary for several factors on the sub-kilometer scale. The waveforms are low-passed filtered at 4 Hz, which is the maximum frequency that is resolvable in our simulations. sideration the integrated effects from other lateral heterogeneities. Within this context, we observe that the largest amount of amplification can be expected at the northwest and southeast sections of the survey, and that the area with the least amount of amplification is concentrated around the NIF. This result is theoretically-consistent with what is expected since, as mentioned above, the NIF zone appears to be composed of faster earth materials and, therefore, has stiffer site conditions. The standard quantity that is used to determine near-surface amplification is Vs30, the shear wave velocity in the top 30 m of the crust [61]. In Figure 5.18, we show the velocity variations averaged from the top two layers in our model (top 25 m) and are about a factor of two higher, which is likely due to averaging from deeper layers as indicated by the sensitivity kernel show in the figure. However, the biggest difference is in the contribution from the topographic slope component of the USGS estimate. This is dominated by the small compact edifice of Signal Hill, which rises 100-m above the surrounding area, and is a pop-up structure caused by convergence across two strands of the Newport-Inglewood Fault. This com- 90 Figure 5.16: Experiment set-up for the estimation of the amplification factors for the case of a vertically propagating shear wave. The source wavelet, and its frequency content, are shown in the upper right panel. parison suggests that Vs30 may not be a useful measure of amplification in the Long Beach area. 5.5 Conclusions In this study, we used ambient noise data recorded at three dense petroleum industry surveys that were deployed in Long Beach, California, to map the small-scale heterogeneities of the crust and investigate their impact on the amplification of seismic waves. For this purpose, we developed a fully-automatic neighborhood cross-correlation method for phase arrival picking that allowed us to perform robust phase velocity measurements across a wide frequency range. We then used these dispersion measurements, together with a timebased beamforming approach, to construct high-resolution phase velocity maps and invert for a 3-D shear wave velocity model of the top kilometer of the Long Beach crust. Our velocity estimates compare well with the structure image from the active-source reflection survey that was done with the array, and also show a strong correlation with the main geological features of the area. We also presented evidence that speak to the possibility of using seismic anisotropy observations to detect small-scale changes on the regional stress regime and illuminate faults whose traces are not necessarily visible at the surface. Lastly, 91 Figure 5.17: Map of relative amplification for the case of a vertically propagating shear wave. These values describe the expected PGA relative to the one that is obtained from a 1-D model. we propagated synthetic wavefields through our velocity model to quantify the scale of variability and intensity of amplification that can be expected across the different parts of the experiments. For the case of surface waves, our results revealed that some of the fine-scale structures that are present in this area are capable of causing complicated wave phenomena that can result in rapid lateral changes on the intensity of shaking of over a factor of ten over a horizontal distance of less than one kilometer. For the case of vertically propagating shear-waves, our results suggest that the mildest shaking can be expected in the vicinity of the Newport-Inglewood fault, as the majority of the fault zone is composed of denser earth materials that were exhumed by the transpressional stresses that are acting on the region. In general, our results not only provide insights into the geo-mechanical properties of sedimentary basins but also highlight the enormous potential of dense nodal arrays to improve the reliability of modern hazard assessments. 92 Figure 5.18: Comparison between the Vs30 and the velocity measurements of this investigation. (a) Shows the Vs30 values across the surveys. (b) Shows our shear wave velocity estimates at 15-m depth. (c) Shows the ratio of (a) and (b). (d) Shows the sensitivity kernels of our measurements for the frequency range of analysis. The shallowest kernel peaks at around 30-m depth (red circle). References [1] Harry O Wood. Preliminary report on the long beach earthquake. Bulletin of the Seismological Society of America, 23(2):43–56, 1933. [2] Carl W Stover and Jerry L Coffman. Seismicity of the United States, 1568-1989 (revised). US Government Printing Office, 1993. 93 [3] Robert Bolin and Lois Stanford. The northridge earthquake: communitybased approaches to unmet recovery needs. Disasters, 22(1):21–38, 1998. [4] Giovanni Costa, Giuliano Francesco Panza, Peter Suhadolc, and Franco Vaccari. Zoning of the italian territory in terms of expected peak ground acceleration derived from complete synthetic seismograms. Journal of applied geophysics, 30(1-2):149–160, 1993. [5] David J Wald and Robert W Graves. The seismic response of the los angeles basin, california. Bulletin of the Seismological Society of America, 88(2):337–356, 1998. [6] Arthur J Rodgers, Arben Pitarka, N Anders Petersson, Björn Sjögreen, and David B McCallen. Broadband (0–4 hz) ground motions for a magnitude 7.0 hayward fault earthquake with three-dimensional structure and topography. Geophysical Research Letters, 45(2):739–747, 2018. [7] Robert W Graves. The seismic response of the san bernardino basin region during the 2001 big bear lake earthquake. Bulletin of the Seismological Society of America, 98(1):241–252, 2008. [8] KB Olsen, SM Day, LA Dalguer, J Mayhew, Y Cui, J Zhu, VM CruzAtienza, D Roten, P Maechling, TH Jordan, et al. Shakeout-d: Ground motion estimates using an ensemble of large earthquakes on the southern san andreas fault with spontaneous rupture propagation. Geophysical Research Letters, 36(4), 2009. [9] Robert Graves and Arben Pitarka. Kinematic ground-motion simulations on rough faults including effects of 3d stochastic velocity perturbations. Bulletin of the Seismological Society of America, 106(5):2136–2153, 2016. [10] Ricardo Taborda, Shima Azizzadeh-Roodpish, Naeem Khoshnevis, and Keli Cheng. Evaluation of the southern california seismic velocity models through simulation of recorded events. Geophysical Journal International, 205(3):1342–1364, 2016. [11] Xin Wang and Zhongwen Zhan. Moving from 1-d to 3-d velocity model: automated waveform-based earthquake moment tensor inversion in the los angeles region. Geophysical Journal International, 220(1):218–234, 2020. [12] Christine A Goulet, Curt B Haselton, Judith Mitrani-Reiser, James L Beck, Gregory G Deierlein, Keith A Porter, and Jonathan P Stewart. Evaluation of the seismic performance of a code-conforming reinforcedconcrete frame building—from seismic hazard to collapse safety and economic losses. Earthquake Engineering & Structural Dynamics, 36(13): 1973–1997, 2007. 94 [13] Carl Tape, Qinya Liu, Alessia Maggi, and Jeroen Tromp. Adjoint tomography of the southern california crust. Science, 325(5943):988–992, 2009. [14] Minyan Zhong and Zhongwen Zhan. An array-based receiver function deconvolution method: methodology and application. Geophysical Journal International, 222(1):1–14, 2020. [15] Jack B Muir and Victor C Tsai. Geometric and level set tomography using ensemble kalman inversion. Geophysical Journal International, 220 (2):967–980, 2020. [16] Harold Magistrale, Steven Day, Robert W Clayton, and Robert Graves. The scec southern california reference three-dimensional seismic velocity model version 2. Bulletin of the Seismological Society of America, 90(6B): S65–S76, 2000. [17] MD Kohler, H Magistrale, and RW Clayton. Mantle heterogeneities and the scec reference three-dimensional seismic velocity model version 3. Bulletin of the Seismological Society of America, 93(2):757–774, 2003. [18] En-Jui Lee, Po Chen, Thomas H Jordan, Phillip B Maechling, Marine AM Denolle, and Gregory C Beroza. Full-3-d tomography for crustal structure in southern california based on the scattering-integral and the adjointwavefield methods. Journal of Geophysical Research: Solid Earth, 119(8): 6421–6451, 2014. [19] Guibao Liu, Patricia Persaud, and Robert W Clayton. Structure of the northern los angeles basins revealed in teleseismic receiver functions from short-term nodal seismic arrays. Seismological Research Letters, 89(5): 1680–1689, 2018. [20] Fan-Chi Lin, Dunzhu Li, Robert W Clayton, and Dan Hollis. Highresolution 3d shallow crustal structure in long beach, california: Application of ambient noise tomography on a dense seismic array. Geophysics, 78(4):Q45–Q56, 2013. [21] Zhe Jia and Robert W Clayton. Determination of near surface shear-wave velocities in the central los angeles basin with dense arrays. Journal of Geophysical Research: Solid Earth, page e2020JB021369, 2021. [22] Thomas L Wright. Structural geology and tectonic evolution of the los angeles basin, california: Chapter 3: Part 1. 1991. [23] Zefeng Li, Zhigang Peng, Dan Hollis, Lijun Zhu, and James McClellan. High-resolution seismic event detection using local similarity for large-n arrays. Scientific reports, 8(1):1–10, 2018. 95 [24] Asaf Inbal, Jean Paul Ampuero, and Robert W Clayton. Localized seismic deformation in the upper mantle revealed by dense seismic arrays. Science, 354(6308):88–92, 2016. [25] Nori Nakata, Jason P Chang, Jesse F Lawrence, and Pierre Boué. Body wave extraction and tomography at long beach, california, with ambientnoise interferometry. Journal of Geophysical Research: Solid Earth, 120 (2):1159–1173, 2015. [26] Jorge C Castellanos, Robert W Clayton, and Alan Juarez. Using a timebased subarray method to extract and invert noise-derived body waves at long beach, california. Journal of Geophysical Research: Solid Earth, 125 (5):e2019JB018855, 2020. [27] Robert W Clayton. A detailed image of the continent-borderland transition beneath long beach, california. Geophysical Journal International, 222(3):2102–2107, 2020. [28] Nikolai M Shapiro and Michel Campillo. Emergence of broadband rayleigh waves from correlations of the ambient seismic noise. Geophysical Research Letters, 31(7), 2004. [29] Nikolai M Shapiro, Michel Campillo, Laurent Stehly, and Michael H Ritzwoller. High-resolution surface-wave tomography from ambient seismic noise. Science, 307(5715):1615–1618, 2005. [30] GD Bensen, MH Ritzwoller, MP Barmin, A Lin Levshin, Feifan Lin, MP Moschetti, NM Shapiro, and Yanyan Yang. Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements. Geophysical Journal International, 169(3):1239–1260, 2007. [31] Huajian Yao, Robert D van Der Hilst, and Maarten V De Hoop. Surfacewave array tomography in se tibet from ambient seismic noise and twostation analysis—i. phase velocity maps. Geophysical Journal International, 166(2):732–744, 2006. [32] MP Moschetti, MH Ritzwoller, and NM Shapiro. Surface wave tomography of the western united states from ambient seismic noise: Rayleigh wave group velocity maps. Geochemistry, Geophysics, Geosystems, 8(8), 2007. [33] Fan-Chi Lin, Morgan P Moschetti, and Michael H Ritzwoller. Surface wave tomography of the western united states from ambient seismic noise: Rayleigh and love wave phase velocity maps. Geophysical Journal International, 173(1):281–298, 2008. [34] Michael H Ritzwoller, Fan-Chi Lin, and Weisen Shen. Ambient noise tomography with a large seismic array. Comptes Rendus Geoscience, 343 (8-9):558–570, 2011. 96 [35] Heather Nicolson, Andrew Curtis, Brian Baptie, and Erica Galetti. Seismic interferometry and ambient noise tomography in the british isles. Proceedings of the Geologists’ Association, 123(1):74–86, 2012. [36] George A McMechan and Mathew J Yedlin. Analysis of dispersive waves by wave field transformation. Geophysics, 46(6):869–874, 1981. [37] Matthew M Haney and Victor C Tsai. Perturbational and nonperturbational inversion of rayleigh-wave velocities. Geophysics, 82(3):F15–F28, 2017. [38] John H Shaw, Andreas Plesch, Carl Tape, M Peter Suess, Thomas H Jordan, Geoffrey Ely, Egill Hauksson, Jeroen Tromp, Toshiro Tanimoto, Robert Graves, et al. Unified structural representation of the southern california crust and upper mantle. Earth and Planetary Science Letters, 415:1–15, 2015. [39] Pierre Gouédard, Huajian Yao, Fabian Ernst, and Robert D van der Hilst. Surface wave eikonal tomography in heterogeneous media using exploration data. Geophysical Journal International, 191(2):781–788, 2012. [40] Ge Jin and James B Gaherty. Surface wave phase-velocity tomography based on multichannel cross-correlation. Geophysical Journal International, 201(3):1383–1398, 2015. [41] JC VanDecar and RS Crosson. Determination of teleseismic relative phase arrival times using multi-channel cross-correlation and least squares. Bulletin of the Seismological Society of America, 80(1):150–169, 1990. [42] MP Barmin, MH Ritzwoller, and AL Levshin. A fast and reliable method for surface wave tomography. In Monitoring the comprehensive nucleartest-ban treaty: Surface waves, pages 1351–1375. Springer, 2001. [43] D Gish and S Boljen. Personal Communication, 2021. [44] GHF Gardner, LW Gardner, and AR Gregory. Formation velocity and density—the diagnostic basics for stratigraphic traps. Geophysics, 39(6): 770–780, 1974. [45] Joseph Fairfield Poland. Ground-water geology of the coastal zone, Long Beach-Santa Ana area, California, volume 1109. US Government Printing Office, 1956. [46] Lorraine A Leon, James F Dolan, John H Shaw, and Thomas L Pratt. Evidence for large holocene earthquakes on the compton thrust fault, los angeles, california. Journal of Geophysical Research: Solid Earth, 114 (B12), 2009. 97 [47] Steven G Wesnousky. 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. Journal of Structural Geology, 27(8):1505–1512, 2005. [48] Thomas P Harding. Seismic characteristics and identification of negative flower structures, positive flower structures, and positive structural inversion. AAPG Bulletin, 69(4):582–600, 1985. [49] Lei Yang, Xin Liu, and Gregory C Beroza. Revisiting evidence for widespread seismicity in the upper mantle under los angeles. Science Advances, 7(4):eabf2862, 2021. [50] Stuart Crampin. The fracture criticality of crustal rocks. Geophysical Journal International, 118(2):428–438, 1994. [51] A Arda Ozacar and George Zandt. Crustal structure and seismic anisotropy near the san andreas fault at parkfield, california. Geophysical Journal International, 178(2):1098–1104, 2009. [52] A Licciardi, T Eken, T Taymaz, N Piana Agostinetti, and S YolsalÇevikbilen. Seismic anisotropy in central north anatolian fault zone and its implications on crustal deformation. Physics of the Earth and Planetary Interiors, 277:99–112, 2018. [53] Martin L Smith and FA Dahlen. The azimuthal dependence of love and rayleigh wave propagation in a slightly anisotropic medium. Journal of Geophysical Research, 78(17):3321–3333, 1973. [54] George Backus. A geometrical picture of anisotropic elastic tensors. Reviews of geophysics, 8(3):633–671, 1970. [55] Wenzheng Yang and Egill Hauksson. The tectonic crustal stress field and style of faulting along the pacific north america plate boundary in southern california. Geophysical Journal International, 194(1):100–117, 2013. [56] Michael Afanasiev, Christian Boehm, Martin van Driel, Lion Krischer, Max Rietmann, Dave A May, Matthew G Knepley, and Andreas Fichtner. Modular and flexible spectral-element waveform modelling in two and three dimensions. Geophysical Journal International, 216(3):1675–1692, 2019. [57] David G Harkrider. Surface waves in multilayered elastic media i. rayleigh and love waves from buried sources in a multilayered elastic half-space. Bulletin of the Seismological Society of America, 54(2):627–679, 1964. 98 [58] Robert Clayton, Thomas Heaton, Mani Chandy, Andreas Krause, Monica Kohler, Julian Bunn, Richard Guy, Michael Olson, Matthew Faulkner, MingHei Cheng, et al. Community seismic network. Annals of Geophysics, 54(6):738–747, 2011. [59] CJ Wills, M Petersen, WA Bryant, M Reichle, GJ Saucedo, S Tan, G Taylor, and J Treiman. A site-conditions map for california based on geology and shear-wave velocity. Bulletin of the Seismological Society of America, 90(6B):S187–S208, 2000. [60] Daniel C Bowden and Victor C Tsai. Earthquake ground motion amplification for surface waves. Geophysical Research Letters, 44(1):121–127, 2017. [61] EM Thompson, David J Wald, and CB Worden. Avs 30 map for california with geologic and topographic constraints. Bulletin of the Seismological Society of America, 104(5):2313–2321, 2014. 99 Chapter 6 THE FINE-SCALE STRUCTURE AT LONG BEACH, CALIFORNIA: BODY WAVES. This chapter was adapted from: J. C. Castellanos, R. W. Clayton, and A. Juarez (2020). “Using a timebased subarray method to extract and invert noise-derived body waves at Long Beach, California.” In: Journal of Geophysical Research 125.5. doi: 10.1029/2019JB018855. 6.1 Abstract The reconstruction of body waves from the cross-correlation of random wavefields has recently emerged as a promising approach to probe the fine-scale structure of the Earth. However, because of the nature of the ambient noise field, the retrieval of body waves from seismic noise recordings is highly challenging and has only been successful in a few cases. Here, we use seismic noise data from a 5,200-node oil-company survey to reconstruct body waves and determine the velocity structure beneath Long Beach, California. To isolate the body wave energy from the ambient noise field, we divide the entire survey into small-aperture subarrays and apply a modified double-beamforming scheme to enhance coherent arrivals within the cross-correlated waveforms. The resulting beamed traces allow us to identify clear refracted P-waves traveling between different subarray pairs, which we then use to construct a high-resolution 3D velocity model of the region. The inverted velocity model reveals velocity variations of the order of 3% and strong lateral discontinuities caused by the presence of sharp geologic structures such as the Newport Inglewood fault (NIF). Additionally, we show that the resolution that is achieved through the use of high-frequency body waves allows us to illuminate small geometric variations of the NIF that were previously unresolved with traditional passive imaging methods. 6.2 Introduction Detailed knowledge of the subsurface velocity structure is essential for predicting ground motion and, thus, earthquake-hazard assessment. Over the past 100 few decades, significant advancements in imaging the Earth’s interior have been possible with the advent of ambient noise cross-correlation. This method goes beyond the spatial limitations of classical earthquake seismology and uses the Earth’s background vibrations (i.e. ambient noise field) recorded at a pair of synchronous seismic stations to reconstruct the Green’s function between the two stations [1–3]. Both the robustness and practicality of this technique have allowed the seismological community to investigate the mechanical properties of the shallow earth on regional [e.g., 4, 5] and continental scales [e.g., 6–8]. However, as the number of high-density seismic networks continues to increase, it is being revealed that near-surface velocities can have large lateral contrasts that result in peak ground acceleration variations of a factor of 5 over a horizontal length scale of less than a kilometer [9]. This scale of variations motivates our necessity to push the limits of modern seismic imaging and explore the use of new generation instrumentation and processing techniques to resolve the fine-scale velocity structure of the Earth. During the first half of 2011, an unprecedentedly dense seismic network was deployed in Long Beach, California, area as part of a petroleum industry survey. This network consisted of more than 5,200 vertical velocity sensors that were distributed over a 7x10-km area at the outboard transition from continental southern California to the Inner Borderland, across the Newport Inglewood fault system (Figure 6.1). While the general purpose of this survey was to perform conventional active-source imaging, Lin et al. [10] showed that highfrequency (0.5-4 Hz) ambient noise surface waves recorded by the network could also be used to construct a detailed 3D shear-wave velocity model for the top 800 m of the crust. The clear correlation of the surface wave velocity model with the known geologic structure proved that passive seismic data can constrain the sub-kilometer-scale material properties of the subsurface. More recent work has shown that this type of industry-based instrumentation can further be used to investigate local microseismicity [e.g., 11], the mechanics of active fault zones [e.g., 12] and diverse characteristics of the ambient seismic wavefield [e.g., 13]. As the use of dense seismic networks grew more common, it also became evident that it is possible to extract body waves from the cross-correlation of the ambient noise field [13, 15, 16]. This topic has attracted considerable attention from the geophysical community since body waves penetrate deep into the 101 Figure 6.1: Regional map of southern California and the Continental Borderland. The black rectangle marks the location of the Long Beach array and the inset map shows the distribution of the seismic sensors (as black circles). The mapped faults are from [14]. The main trace of the Newport Inglewood fault is labeled as NIF. The major strands of the Newport Inglewood fault are denoted by red lines in the magnified map. To provide some context on the level of density of the Long Beach array, the stations of the Southern California Seismic Network are plotted as inverted red triangles. Earth and can probe the seismic structure with higher resolution than surface waves [e.g., 17, 18]. However, because ambient noise energy propagates mainly horizontally [3, 4, 19], reconstructing body waves from seismic noise recordings has proven to be a difficult task [20, 21]. In a recent study, Nakata et al. [17] used a 2,500-station array adjacent to the Long Beach survey to extract body waves and performed the first-ever body wave tomography using pure ambient noise recorded at the ground surface. This investigation, however, used a series of selection filters on individual waveforms to isolate the body wave energy, which, based on a set of quality control criteria, only allowed the use of a small portion of the entire dataset (about 35%). In this study, we extend the work of Nakata et al. [17] and use the Long Beach network to extract ambient noise body waves and map the velocity structure 102 beneath the array. First, we validate the existence of coherent body wave arrivals in the cross-correlated waveforms by generating a stacked record section of the Long Beach dataset, inverting for a 1D velocity model, and computing a synthetic wavefield. We then divide the entire array into several small-aperture subarrays and apply an array processing technique to retrieve first-arrival body waves traveling between the different groups of stations. Finally, we make reliable traveltime measurements and perform a 3D traveltime tomography to resolve the P-wave velocity variations of the uppermost part of the crust. 6.3 Ambient Noise Correlation and Body Waves For the Long Beach survey, Lin et al. [10] built a complete set of more than 5,200 virtual sources (i.e. one for each physical receiver) by cross-correlating the receiver at the source location and every other receiver in the array. The processing scheme used to generate such a dataset followed the method of Bensen et al. [22] (without the temporal normalization step) and yielded more than 13.5 million ambient noise cross-correlations. Aside from being able to retrieve clear fundamental Rayleigh waves, Lin et al. [10] showed that weak yet coherent body waves existed in the Long Beach correlograms and that they could be readily discerned above the noise after stacking over all receiver pairs into discrete offset bins (Figure 6.2). These arrivals, however, are not prominent enough to be reliably identified in individual traces. In this study, we use the same gathers and array-processing tools to extract body waves traveling between small groups of stations. A close inspection into the first 5 seconds of the stacked correlated waveforms reveals that body waves start propagating at approximately time zero with velocities larger than 1.7 km/s (Figure 6.3A). As in the work of Nakata et al. [17], we find that the apparent wave speed of the first arrival increases with offset, which suggests that the reconstructed wave corresponds to a refracted or diving P-wave. Moreover, we observe 2 distinct pulses arriving shortly after the first arrival that have slightly different horizontal slownesses and, hence, different propagation paths. To determine the nature of these arrivals, we convert the P-wave portion of the stacked cross-correlations into the T au-p domain and use this representation to invert for a smoothly variant 1D velocity model (Figure 6.4) [23]. We then use the inverted velocity profile to perform a 2D finite difference simulation [24] and resolve whether these arrivals are the result of simple 1D wave propagation or due to large lateral velocity variations 103 Figure 6.2: Stacked record section of the Long Beach virtual shot gathers. The waveforms are band-passed filtered between 2 to 8 Hz. The bin size of the spatial stacking is of 50 m. Clear fundamental mode Rayleigh waves, some combination of higher-mode Rayleigh waves and S-waves, and P-waves are visible at almost all offset ranges. The dashed black lines mark the different time windows in which each of these phases are expected to arrive [from 10]. within the array. Our analysis shows that a simple 1D model is sufficient to replicate all 3 arrivals and that they correspond to a diving P-wave, a PP-wave and a direct P-wave that is propagating close to the surface (Figure 6.3B-C). The strong resemblance between the synthetic and observed waveforms proves the correlation’s convergence to the Green’s function and confirms the presence of coherent body wave energy in the Long Beach correlograms. 104 Figure 6.3: Observed and synthetic wavefields windowed around the expected P-wave time of arrival. (A) Record section of the first 5 seconds of the stacked cross-correlations. The thick black lines mark 3 different reference velocities that approximate the variations in the apparent wave speed of the first-arrival P-wave as a function of offset. (B) Synthetic wavefield computed by 2D finite difference modeling using the inverted 1D velocity model (black profile in Figure 6.4B). The frequency range of the waveforms is the same as Figure 6.2. (C) Snapshots of the 2D velocity wavefield at different times (1.5 s, 2.5 s, and 3.5 s). The explosion at the surface marks the location of the source. Note that attenuation was not considered in the wavefield simulation. 6.4 Double Beamforming and Body Wave Extraction Although P-waves are observable in the stacked correlograms, the generation of a high-resolution velocity model requires a complete characterization of the different wave propagation properties within the seismic survey. To this end, we apply a slightly modified double-beamforming (DBF) scheme to extract first-arrival P-waves propagating between different sections of the array. The DBF method [26] combines classical slant stack processing [27–29] on a source and receiver array to reconstruct a beam or eigenray between the center-point of the two groups of stations. In short, this technique entails finding the appro- 105 Figure 6.4: T au-p representation of the first 5 seconds of the stacked crosscorrelated waveforms (A) and inverted velocity model (B). As a comparison, the model that is obtained from the inversion of the T au-p curve (black line in A) is shown together with a 1D profile of the CVM-S4 model [25] at the center of the Long Beach array. priate slowness, u, and azimuthal direction, θ, of a given phase simultaneously on both sides by applying a systematic delay-and-sum to all the recordings. Following the notation of Nakata et al. [18], the double-beam space, B, at a given time, t, can be constructed by scanning the slowness and azimuth domains at the source and receiver locations through the computation of B(us , θs , ur , θr , t) = 1 XXXX C(xs , ys , xr , yr , t− Ns Nr xs ys xr yr (6.1) τs (xs , ys , us , θs ) + τr (xr , yr , ur , θr )), where the subscripts s and r refer to the sources and receivers, respectively, N is the number of source and receivers, C is the cross-correlation function between the spatial points x and y (e.g. eastward and northward), and τ is the relative time lag from a reference point, which, for a 2D seismic array, can be defined as τ (x, y, u, θ) = u(x − xc ) sin(θ) + u(y − y c ) cos(θ), (6.2) wherein the coordinates xc and y c represent the center of an array. Because of the large number of traces that are generally involved in the construction of a double-beam, the improvement on the signal-to-noise ratio is significant and, 106 therefore, ideal to bring out the weak body waves that are buried within the correlogram’s noise. As in every array-processing problem, it is desirable to explore the complete model space to find the optimal pairs of slownesses and directional azimuths that would result in the most energetic and coherent beam. However, due to the large number of instruments in the Long Beach survey, scanning through the 5D space of the double-beam for every possible source and receiver array configuration is computationally very demanding and basically impractical. To overcome this challenge, we apply two modifications to the traditional DBF method. First, we assume that, because of the relatively short interstation distances of the array, the diving P-waves do not suffer large deviations from their direct raypaths. This assumption allows us to collapse the directional azimuth dimensions, θs and θr , by fixing their optimal values to the azimuth between the geographic center of the source and receiver arrays. Synthetic tests using the real Long Beach station geometry indicate that this simplification is reasonably valid even for lateral velocity contrasts as large as 20%. Second, instead of scanning through the remaining slowness and time space, we compute the time-frequency spectra of each individual record via an S-transform [30] and, from this representation, extract the timing of every energy peak that is ±0.2 seconds within the predicted arrival time of the diving P-wave (Figure 6.5). The length of the tolerance time window is chosen based on the fact that we are focusing on body waves that are mostly higher than 3 Hz [17] and we only want to include a maximum of 3 wavelets to avoid severe cycle skipping. Subsequent analyses of our DBF results suggests that most of the shallow crustal velocity variations of Long Beach can be captured with this window length and that its broadening does not change the first-order distribution of our traveltime observations. As an aside, it is worth noting that here we do not use nor retain any frequency information that is contained in the S-transform spectrograms as we opt to degrade their frequency resolution in order to enhance their temporal resolution and obtain finer traveltime picks [31]. After collecting the time of arrival of every potential refracted P-wave, we cast the entire double-beam calculation as a linear inverse problem and simultaneously solve for the source and receiver side slownesses that best predict the traveltime dataset via weighted least-squares. In other words, we minimize 107 Figure 6.5: Examples of traveltime picking of individual traces. (A) Map of the Long Beach array configuration. The red lines mark the shortest path between the virtual source (station S) and the receivers (stations R1 and R2) involved in the computation of the correlograms shown in (B) and (C). (B) Cross-correlated waveforms between stations S and R1 together with its timefrequency spectra computed via the S-transform. The dashed red line on top of the waveforms marks the predicted arrival time of the P-wave. The star markers in the lower panel depict all of the energy peaks identified in the time-frequency spectra and the red box delimits the time window used to identify the potential P-wave arrivals. The energy peaks in the spectrogram are determined using a watershed algorithm [32]. (C) Same as (B) but for the correlogram between stations S and R2. the prediction error E = eT We e and solve, m = [GT WG]−1 GT Wd, (6.3) where m is a vector containing the source-side and receiver-side slowness, G is a sensitivity matrix that relates the source and receiver array geometries to the traveltime data, W is a weighting matrix, and d is a vector containing the traveltime picks. To weight the inversion, we use the absolute amplitude of the picked energy contours of the S-transformed correlograms. It is important to point out that the amplitudes of the P-wave signals largely depend on the local noise distribution and are only used here to prioritize the alignment of 108 energetic arrivals rather than the one of weak and random oscillations within the traces. This last step turns our processing scheme into something that is equivalent to finding the source-side and receiver-side velocities that best maximizes the energy of the stacked waveforms around the predicted time of arrival of the refracted P-wave. To implement the DBF scheme described above, we first downsample the data to a 0.02 s sampling and band-pass filter every available correlogram between 2-8 Hz. We then loop through every virtual source and construct source arrays by taking their 29 nearest neighbors. Then, for each source array, we loop through every other station that is not part of the source group and take their 29 nearest neighbors to construct the receiver arrays. Each source and receiver array is, therefore, composed by 30 stations each, which results in a total of 900 Green’s functions for each subarray pair. The number of stations that are involved in the construction of the subarrays is determined empirically based on the tradeoff that, to better capture the small-scale velocity variations within the survey, we want to use the least amount of instruments in the DBF and still be able to observe clear body wave arrivals. As a quality control, we only apply the DBF to source and receiver groups in which both array radii were below 900 m and their geographic centers are at least 1 km apart. The first condition is set to ensure strong waveform coherence between the traces while the second condition is set to allow some time separation between the body and surface wave arrivals. Moreover, to assure that the reconstructed beams are statistically significant, we apply a bootstrapping method to calculate 95% confidence limits, and choose to only retain beams in which their confidence region is smaller than 0.3 km2 /s2 using a total of 500 picked traveltime resamples. The size of the threshold for the confidence region is determined based on a visual inspection of the quality of the beams and the overall bootstrapping results. Figure 6.6 shows an example of the modified double-beamforming scheme for a single source-receiver array configuration. A unique aspect of our processing scheme is that the size and geometry of the source and receiver arrays changes every time a beam is constructed at different locations across the survey. This particularity causes the resolution of the DBF, which largely depends on the size of the subarrays, to vary spatially such that its ability to extract wavelets of specific frequencies is not always the same. To prevent spatial aliasing, the general rule in beamforming is that the 109 Figure 6.6: Double-beamforming process for a single source-receiver array configuration. (A) Map showing the entire Long Beach survey together with a source array (red circles) and a receiver array (gray circles). The red radius around the source array marks the 1-km distance threshold imposed in our stacking scheme. For this particular source-receiver configuration, the source radius is of 0.32 km, the receiver radius is of 0.42 km, and the distance between the arrays’ geographic centers is of 2.12 km. (B) Entire source and receiver side velocity space of the DBF. The yellow star indicates the best-fitting velocity pair obtained in the least-squares inversion. The black dots correspond to the best-fitting velocity pairs obtained in the bootstrapping process and the red box marks their 95% confidence limits. For this case, the area of the 95% confidence region is 0.09 km2 /s2 . (C) Empirically derived probability density function of the source-side velocity (SSV) and the receiver-side velocity estimations (RSV) calculated from the bootstrapped results. (D) Beamed trace (black waveforms) plotted behind one of the 900 original correlograms used to build the beam (red waveforms). Note how the P-wave is essentially unidentifiable in the original correlogram and how its signal-to-noise ratio improves dramatically after applying DBF. The retrieved phase corresponds to a diving P-wave propagating between the center of the source and receiver arrays. aperture of the array must be at least >1 wavelength than the longest period of interest and the average interstation spacing must be <0.5 wavelengths [18, 33, 34]. For our case, the Long Beach survey has an average interstation distance of 100 m and, as stated above, all of our source and receiver arrays 110 are composed of 30 stations each. This geometry causes the average aperture of most subarrays to be at around 800 m. Now, assuming that the slowest body wave propagates as a perfectly plane wave with an apparent velocity of 1.5 km/s, the minimum size of the subarrays capable of extracting >3 Hz seismic arrivals is ∼450 m. This result allows us to conclude that our choice of subarray geometries is optimal for extracting the body wave arrivals that are buried within the Long Beach correlograms, and that such configuration would technically allow us to beamform up to 7.5 Hz waves without being spatially aliased. The application of the modified DBF to the entire Long Beach dataset yielded approximately 12 million unique beams that now show enhanced diving Pwaves. As an example, Figure 6.7 presents snapshots of the constructed wavefield before and after the application of DBF. A clear body wave propagating away from the virtual source in all directions can be seen after the waveforms have been coherently stacked. It is interesting to notice how the wavefield is not completely spherical, which indicates the existence of lateral variations in the elastic properties of the medium. To further illustrate this point, we build two orthogonal refraction profiles by stacking all available beams along two corridors that are 0.5 km wide in a 50-m offset bin (Figure 6.8). Aside from showing clear diving P-waves regardless the direction of propagation, both profiles show a prominent step in the traveltime that is typical of a horizontal discontinuity in a buried layer. For the case of the N-S profile, the velocity jump is coincident with the surface trace of the Newport Inglewood fault. This observation supports the premise that we are observing real features of the data and that the constraints imposed in the DBF are flexible enough to capture the variations in the structural properties of the Long Beach crust. 6.5 Traveltime Measurements and Body Wave Tomography After we successfully apply the DBF to every possible source-receiver configuration, we pick the arrival times of the beamed P-waves. For this matter, we first construct a set of reference traces by stacking all available beams in a 50-m offset bin (Figure 6.9A). We then manually identify the diving P-waves in every reference trace and extract their waveform shape so that they can be used as templates to automatically find their time of arrival in individual beams with similar offsets. To extract the waveform shape of these arrivals, we use the seismogram decomposition method of Juarez and Jordan [35]. This 111 Figure 6.7: Wavefield emitted by a virtual source. The location of the virtual source is marked by the yellow star. The lag-times are shown on the upper left corner of each panel. The top panels show snapshots of the 2-8 Hz band-passed wavefield observed at each station. The middle panels show snapshots of the beamed wavefield observed at each station. Strong body waves are now visible propagating away from the virtual source. The lower panels shows snapshot of a synthetic wavefield using the inverted velocity model (black profile in Figure 6.4B). Note how the observed wavefield is not completely spherical as in the synthetic case. technique takes the time-frequency spectra of a given record computed by the S-transform and systematically isolates all apparent arrivals using localized Gaussian filters around every local maxima in the power spectrum. It then 112 Figure 6.8: Examples of P-wave refractions from a N-S line (A) and a E-W line (B). The profiles are formed by stacking all available beams along corridors that are 0.5 km wide (green areas in C) in a 50-m offset bin. The location of the virtual sources are marked by the yellow stars. Note how there exists a clear jump in the N-S profile that is coincident with the Newport Inglewood fault (thick red line in all panels), and shorter jump in the far end of the E-W profile that might be associated with a smaller fault (thick blue line in panel B). integrates every filtered spectrum with time and separately transforms them back into the time domain, thus yielding a finite set of waveforms that are localized in time and frequency. As an example, Figure 6.9B-D presents the waveform decomposition process for a single reference trace, and shows how one can recover the input trace after summing all of its decomposed elements. It is important to note that, for this particular reference trace, the temporal resolution of the S-transform is sufficient to isolate the diving P-wave from the PP- and direct P-waves. Although the latter may not be the case for shorter off-set traces, the waveform decomposition scheme ultimately allows us to re- 113 trieve the characteristic shape of the first-arrival whether it consists of a single diving P-wave or a superposition of different P-wavelets. Figure 6.9: Stacked record section of beams (A) and waveform decomposition process of a single reference trace (B-D). (A) Stacked recordsection showing how the DBF has successfully isolated the body wave energy and suppressed any other arrival that is not arriving with a slowness of a refracted P-wave. The dashed colored lines mark the time window where the higher-mode surface waves were present before the DBF. (B) Reference trace (black waveform) and reconstructed trace (red waveform) that is obtained after linearly adding all of the decomposed waveforms (all traces in D). The reference trace corresponds to the stacked beam in the 7 km offset bin (red arrow in A). (C) Time-frequency spectra of the reference trace computed via the S-transform. The star markers indicate the arrivals identified in the power spectrum and the numbers indicate their order in time. (D) Decomposed traces sorted by time. A total of 27 apparent phases are identified in the reference trace spectra, with the diving P-wave being the most prominent amongst all. The green trace alone is what is used as a template to pick the refracted P-wave arrival in individual beams that have a 7-km offset. Once we build a library of P-wave templates, we apply the same phase decomposition scheme to all individual beams and, for each one of them, retain the timing of the isolated phase that is within the expected time of arrival of the diving P-wave and best correlates with its respective waveform template. To show the performance of our picking method, Figure 6.10 presents the distribution of the traveltime differences between the picked time and the theoretical arrival (from the inverted 1D model) for 2 virtual sources located in opposite sides of the survey. The spatial coherency of the traveltime anomaly maps, together with the general high correlation coefficients (CCs) obtained in 114 the picking process, subjectively attest to the reliability of our measurements and suggest that their smooth fluctuations are likely to be associated with variations in the local velocity structure. Moreover, there is a clear northsouth velocity dichotomy with the transition just on the surface trace of the Newport Inglewood fault that is revealed by the northern virtual source. It is also worth noticing that stations in which a low correlation coefficient were obtained are in agreement with beams that have a low signal-to-noise ratio (SNR) and regions where traveltime measurements have sudden and unrealistic deviations from their neighboring stations. The majority of these random measurements appear to be correspondent to beams that are constructed from either large-offset correlations or from stations that are located in the vicinity of the fault. Furthermore, a general analysis of the CC and SNR distribution of different virtual sources reveals that higher-quality beams are mostly constructed from subarrays that are located on similar sides of the fault. This observation suggests that there must exist some zone of structural complexity along the fault that may be acting as a barrier to wave propagation. Because of the extensively large number of traveltime measurements that are available to estimate the 3D velocities, solving the structural problem using a regularized inversion, which involves calculating and storing the data sensitivity (Fréchet) kernel for each measurement, is an intense computational task. For this reason, we opt to implement an iterative backprojection method [36] that maps, on a ray-by-ray basis, traveltime anomalies into slowness perturbations along the ray paths until the data are satisfied. To parameterize the model, we use a regularly spaced 70-m grid that extends 8 km (easting) x 12 km (northing) x 2.5 km (depth) and the inverted 1D velocity profile as a starting model. We then trace every ray through the model using the fast sweeping method [37] and average the perturbations applied to each model element from all the rays that are influenced by that model element. Once the entire model is updated, we retrace the rays, compute new theoretical traveltimes and remap the new measurements into slowness perturbations. This process is iterated until the updates to the model start to become negligible. It is important to note that the DBF that was applied to the traces will cause the resulting model to be inherently smooth. Figure 6.11 shows the horizontal and vertical slices of the obtained velocity model after 5 iterations of the backprojection method. To ensure robustness in the construction of this model, we only used traveltime measurements that had a correlation coefficient larger 115 Figure 6.10: Traveltime measurements for two virtual sources localized on opposite sides of the survey. The left panels show the distribution of the traveltime differences between the picked time and the theoretical arrival. The middle panels show the normalized correlation coefficient between each picked phase with their respective waveform template. The right panels show the signal-to-noise ratio of the beams. Here, we define the SNR by the ratio between the peak amplitude within a window containing the refracted P-wave arrival to the root-mean-square of the amplitudes before this window. The locations of the virtual sources are marked by the yellow stars. The red lines depict the major strands of the Newport Inglewood Fault and the black dashed lines delimit the regions where traveltime measurements appear to be random. than 0.5 with their respective waveform templates and a signal-to-noise ratio larger than 2. This restriction allowed the use of 11,033,733 rays which, based on the ray coverage maps, we find that can resolve velocities down to a depth of 2 km at the center of the survey. 116 Figure 6.11: Horizontal (A) and vertical (B) slices of the inverted velocity model. Velocities are presented as perturbation (in percent) from the horizontally averaged inverted model. To emphasize the lateral velocity variations, each horizontal slice is plotted with its independent scale bar. The depth of the horizontal slices is given at the upper right corner of each panel (0.14 km, 0.49 km and 0.98 km). The blue lines in the leftmost panel in A show the location of the vertical slices in B (P1, P2, P3 and P4). For this set of plots, the local coordinates have been shifted to the center of the survey. The red lines in every panel mark the major strands of the Newport Inglewood fault. To assess the uncertainty of our estimated velocities, we first generate a set of one-iteration tomograms by bootstrapping through the virtual sources. We then follow a similar approach to the one that is used in traditional Eikonal tomography and, for each point in our model, estimate a mean slowness, s0 , and a standard deviation, σs0 , from the distribution of slownesses, si : n 1X si , s0 = n i=1 (6.4) 117 σs20 = n X 1 (si − s0 ) 2 , n(n − 1) i=1 (6.5) where n is the number of virtual sources [38]. Once these quantities are obtained, we simply compute the uncertainty of the velocity, σv , via: σv = 1 σs . s20 0 (6.6) After resampling the virtual sources 50 times and applying such a process, we find that the uncertainties of our model are below 5 m/s and that they tend to be largest at the edges of the array (Figure 6.12). Note that these uncertainties are related to random noise within our traveltime measurements and not the tomographic inversion itself. Figure 6.12: Depth slices of the velocity uncertainties obtained after bootstrapping the virtual sources 50 times. The depth of the horizontal slices is given at the upper right corner of each panel (0.14 km, 0.49 km, and 0.98 km). 6.6 Results and Discussion The tomographic model shows velocity variations that are, in general, consistent with the surface wave results of Lin et al. [10]. As expected, the most prominent feature within our model is the Newport Inglewood fault zone, which shows up as a high velocity anomaly [probably as a result of strain focusing; 39] that emerges at depths of around 500 m. Although the average velocity structure of this feature has been imaged in previous tomography studies [e.g., 10, 40], the resolution that is now achieved with high-frequency body waves 118 allow us to illuminate some of its geometric variations. Figure 6.11 shows a comparison between different cross-sections of our velocity model cut perpendicular to the main fault trend. From this comparison, we find that there are structural differences amongst the three imaged segments of the fault. For the northwestern strand, the fault trace appears to be mostly vertical and welldefined at shallow depths (Figure 6.13C-D). In contrast, cross-sections across the southeastern strand suggest that the fault extends to deeper depths and might even have a small SW dipping component that grows as it approaches Signal Hill (Figure 6.13F-G). These observations are consistent with reconstructed maps of the internal structures of Los Angeles Basin that suggest that the primary strand of the Newport Inglewood Fault (also known as the Cherry Hill Fault) extends nearly vertical down to depths of 1,100-1,500 m and, at greater depths, may dip as much as 60◦ [41]. Near the central strand, we find that the fault structure is even more complicated and unidentifiable within our velocity model (Figure 6.13E). This complexity may arise since, at this location, the Cherry Hill Fault is known to branch out near the surface through a series smaller-scale fault strands that bound a recently uplifted wedge (Figure 6.13BE). Although the general agreement between our tomographic model with known geologic features attests to the reliability of ambient noise body waves to characterize subsurface velocity structure, we can take advantage of the large station density of the survey to qualitatively validate our results in a different way. Figure 6.14 shows the effects of velocity variations on a Mb2.5 earthquake wavefront propagating through the Long Beach array. From this image, it is clear how the structural variations within the area introduce complexities into the wavefront and can even cause a 180-degree phase change in the S-wave across the Newport Inglewood fault. The partitioning of the wavefront is most extreme in the southeastern strand of the fault, precisely where our velocity model reveals the existence of a prominent and elongated fast velocity anomaly (associated to the fault itself) adjacent to a slow velocity zone that extends to a depth of about 1 km (Figure 6.11A). Moreover, it is important to note that the fast velocities of the fault also causes a slight increase in the wavelength of the surface wave packet, which may have some implications when quantifying the frequency response of the basin at a local scale. The fact that our velocity model captures the structural features that introduce these complexities into the wavefield supports the reliability of our estimations, and 119 Figure 6.13: Fault-perpendicular cross-sections of the inverted velocity model. (A) Map showing the Long Beach survey with the location of the profiles. (B) Cross-sections illustrating the general structural pattern across different segments of the Newport-Inglewood fault zone [modified from 41]. The thick back lines in A show the location of the two profiles (X-X’ and Y-Y’). (C-G) Vertical cross-sections of the inverted velocity model (colored profiles in A). The red arrows on top of each profile indicate the location of the surface trace of the Newport Inglewood fault. The black dashed lines in each profile mark our interpretation of the fault. The inferred fault geometry of Wright [41] across profile X-X’ is overlain on profile E. highlights the relevance of characterizing the fine-scale velocity variations of the shallow crust. The availability of high-resolution velocity models is essential in many geophysical applications and paves the way for more complete interpretation of the Earth’s geologic structures. For this study case, the results are mostly relevant to seismic hazard studies of Los Angeles basin since they can be incorporated into the standard velocity models for the region and, hence, be used in earthquake ground motion predictions. In addition to its role to seismic hazard, our velocity model can be used in an exploration context to derive a set of static corrections that can be applied to seismic reflection analysis and used to improve active-source processing. Static corrections as such have been derived from the surface wave velocity model of Lin et al. [10] and found to be effective in enhancing P-wave reflective signals [42]. Nonetheless, the statics 120 Figure 6.14: Wavefield emitted by the Carson Mb2.5 earthquake at 4 seconds after its origin time (A). Note how the earthquake surface waves clearly break up along the Newport Inglewood fault. The dashed black box in the zoomed panel is the same as in Figure 6.11A and delimits a region where the wavefield is being most advanced and stretched by the high velocities of the fault. The velocity perturbation at 0.98 km depth of this small area is shown in (B). derived from our P-wave velocity model should be even more accurate since they span a larger depth range and are not converted from phase and shear wave velocities. 6.7 Conclusions We used a high-density oil company survey in Long Beach, California, to demonstrate that body waves can be extracted from the ambient noise field and subsequently inverted to produce a high-resolution velocity model of the subsurface. First, we confirmed the existence of coherent body wave energy in the noise correlograms by computing a synthetic wavefield and predicting every prominent arrival that is present in an average record section of the entire dataset. We then retrieved refracted P-waves propagating through different sections of the survey by dividing the Long Beach network into small-aperture subarrays and applying a double-beamforming (DBF) technique that simultaneously solves for the optimal source-side and receiver-side stacking velocities. Profiles of P-wave refractions along particular corridors of the survey after applying DBF revealed that there exist prominent velocity jumps that are likely associated with regional geological features such as the Newport Inglewood fault. After extracting clean refracted P-waves, we measured their absolute traveltime by generating waveform templates and using them to identify their 121 time-of-arrival in individual beams. Once the traveltime measurements were made, we applied an iterative backprojection tomography to map the traveltime measurements into velocity perturbations and solve for the 3D velocity structure beneath the survey. The estimated 3D velocities are in general agreement with previous tomographic results of the area and correlate with the main structural features of the region. A close analysis of the Newport Inglewood fault system reveals that there are structural differences amongst the three imaged segments of the fault. The northwestern strand of the fault seems to extend to depths of around 500 m and appears as a well-imaged vertical velocity discontinuity. On the other hand, the southeastern strand of the fault extends to deeper depths (>1 km) and appears to have a subtle southwest dipping component that grows as the fault approaches Signal Hill. The central strand of the fault is not apparent within our velocity model, which we believe is caused by the structural complexities of the fault zone. In short, the results of this work confirmed the reliability of ambient noise body waves to characterize the seismic structure and showed that, in spite of the processing challenges of retrieving them, they can provide unique constraints that are otherwise unattainable with current passive imaging methods. References [1] Oleg I Lobkis and Richard L Weaver. On the emergence of the green’s function in the correlations of a diffuse field. The Journal of the Acoustical Society of America, 110(6):3011–3017, 2001. [2] Michel Campillo and Anne Paul. Long-range correlations in the diffuse seismic coda. Science, 299(5606):547–549, 2003. [3] Nikolai M Shapiro and Michel Campillo. Emergence of broadband rayleigh waves from correlations of the ambient seismic noise. Geophysical Research Letters, 31(7), 2004. [4] Karim G Sabra, Peter Gerstoft, Philippe Roux, WA Kuperman, and Michael C Fehler. Surface wave tomography from microseisms in southern california. Geophysical Research Letters, 32(14), 2005. [5] Nikolai M Shapiro, Michel Campillo, Laurent Stehly, and Michael H Ritzwoller. High-resolution surface-wave tomography from ambient seismic noise. Science, 307(5715):1615–1618, 2005. 122 [6] Antonio Villasenor, Yingjie Yang, Michael H Ritzwoller, and Josep Gallart. Ambient noise surface wave tomography of the iberian peninsula: Implications for shallow seismic structure. Geophysical Research Letters, 34(11), 2007. [7] Yingjie Yang, Michael H Ritzwoller, Anatoli L Levshin, and Nikolai M Shapiro. Ambient noise rayleigh wave tomography across europe. Geophysical Journal International, 168(1):259–274, 2007. [8] Fan-Chi Lin, Morgan P Moschetti, and Michael H Ritzwoller. Surface wave tomography of the western united states from ambient seismic noise: Rayleigh and love wave phase velocity maps. Geophysical Journal International, 173(1):281–298, 2008. [9] Robert W Clayton, Thomas Heaton, Mani Chandy, Andreas Krause, Monica Kohler, Julian Bunn, Richard Guy, Michael Olson, Mathew Faulkner, MingHei Cheng, et al. Community seismic network. Annals of Geophysics, 54(6):738–747, 2012. [10] Fan-Chi Lin, Dunzhu Li, Robert W Clayton, and Dan Hollis. Highresolution 3d shallow crustal structure in long beach, california: Application of ambient noise tomography on a dense seismic array. Geophysics, 78(4):Q45–Q56, 2013. [11] Zefeng Li, Zhigang Peng, Dan Hollis, Lijun Zhu, and James McClellan. High-resolution seismic event detection using local similarity for large-n arrays. Scientific reports, 8(1):1646, 2018. [12] Asaf Inbal, Jean Paul Ampuero, and Robert W Clayton. Localized seismic deformation in the upper mantle revealed by dense seismic arrays. Science, 354(6308):88–92, 2016. [13] Zack J Spica, Nori Nakata, Xin Liu, Xander Campman, Zijian Tang, and Gregory C Beroza. The ambient seismic field at groningen gas field: An overview from the surface to reservoir depth. Seismological Research Letters, 89(4):1450–1466, 2018. [14] CW Jennings and JS George. Fault activity map of california and adjacent areas, with locations and ages of recent volcanic eruptions, scale 1: 750 000. Divisions of Mines and Geology, Geologic Data Map, (6), 1994. [15] Philippe Roux, Karim G Sabra, Peter Gerstoft, WA Kuperman, and Michael C Fehler. P-waves from cross-correlation of seismic noise. Geophysical Research Letters, 32(19), 2005. [16] Gerrit Olivier, Florent Brenguier, Michel Campillo, Richard Lynch, and Philippe Roux. Body-wave reconstruction from ambient seismic noise correlations in an underground mine. Geophysics, 80(3):KS11–KS25, 2015. 123 [17] Nori Nakata, Jason P Chang, Jesse F Lawrence, and Pierre Boué. Body wave extraction and tomography at long beach, california, with ambientnoise interferometry. Journal of Geophysical Research: Solid Earth, 120 (2):1159–1173, 2015. [18] Nori Nakata, Pierre Boué, Florent Brenguier, Philippe Roux, Valérie Ferrazzini, and Michel Campillo. Body and surface wave reconstruction from seismic noise correlations between arrays at piton de la fournaise volcano. Geophysical Research Letters, 43(3):1047–1054, 2016. [19] Keiiti Aki and Paul G Richards. Quantitative seismology. W.H. Freeman and Company, 2002. [20] Roel Snieder and Eric Larose. Extracting earth’s elastic wave response from noise measurements. Annual Review of Earth and Planetary Sciences, 41, 2013. [21] Nori Nakata and Kiwanu Nishida. Body wave exploration. Seismic Ambient Noise, page 239, 2019. [22] GD Bensen, MH Ritzwoller, MP Barmin, AL Levshin, F Lin, MP Moschetti, NM Shapiro, and Yanyan Yang. Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements. Geophysical Journal International, 169(3):1239–1260, 2007. [23] Peter M Shearer. Introduction to seismology. Cambridge university press, 2019. [24] Dunzhu Li, Don Helmberger, Robert W Clayton, and Daoyuan Sun. Global synthetic seismograms using a 2-d finite-difference method. Geophysical Journal International, 197(2):1166–1183, 2014. [25] John H Shaw, Andreas Plesch, Carl Tape, M Peter Suess, Thomas H Jordan, Geoffrey Ely, Egill Hauksson, Jeroen Tromp, Toshiro Tanimoto, Robert Graves, et al. Unified structural representation of the southern california crust and upper mantle. Earth and Planetary Science Letters, 415:1–15, 2015. [26] Pierre Boué, Philippe Roux, Michel Campillo, and Benoit de Cacqueray. Double beamforming processing in a seismic prospecting context. Geophysics, 78(3):V101–V108, 2013. [27] Jeffrey R Thorson and Jon F Claerbout. Velocity-stack and slant-stack stochastic inversion. Geophysics, 50(12):2727–2741, 1985. [28] M Weber, JP Davis, Ch Thomas, F Krüger, F Scherbaum, J Schlittenhardt, and M Körnig. The structure of the lowermost mantle as determined from using seismic arrays. Seismic modeling of the Earth’s structure, pages 399–442, 1996. 124 [29] Sebastian Rost and Christine Thomas. Array seismology: Methods and applications. Reviews of geophysics, 40(3):2–1, 2002. [30] Robert Glenn Stockwell, Lalu Mansinha, and RP Lowe. Localization of the complex spectrum: the s transform. IEEE transactions on signal processing, 44(4):998–1001, 1996. [31] L Mansinha, RG Stockwell, and RP Lowe. Pattern analysis with two-dimensional spectral localisation: Applications of two-dimensional s transforms. Physica A: Statistical Mechanics and its Applications, 239 (1-3):286–295, 1997. [32] Fernand Meyer. Topographic distance and watershed lines. Signal processing, 38(1):113–125, 1994. [33] Philippe Roux, Bruce D Cornuelle, WA Kuperman, and WS Hodgkiss. The structure of raylike arrivals in a shallow-water waveguide. The Journal of the Acoustical Society of America, 124(6):3430–3439, 2008. [34] Florent Brenguier, Philippe Kowalski, N Ackerley, N Nakata, P Boué, M Campillo, E Larose, S Rambaud, C Pequegnat, T Lecocq, et al. Toward 4d noise-based seismic probing of volcanoes: Perspectives from a largen experiment on piton de la fournaise volcano. Seismological Research Letters, 87(1):15–25, 2015. [35] A Juarez and TH Jordan. Optimization of data functionals for full-3d tomography. In AGU Fall Meeting Abstracts, 2018. [36] JA Hole. Nonlinear high-resolution three-dimensional seismic travel time tomography. Journal of Geophysical Research: Solid Earth, 97(B5):6553– 6562, 1992. [37] Hongkai Zhao. A fast sweeping method for eikonal equations. Mathematics of computation, 74(250):603–627, 2005. [38] Fan-Chi Lin, Michael H Ritzwoller, and Roel Snieder. Eikonal tomography: surface wave tomography by phase front tracking across a regional broad-band seismic array. Geophysical Journal International, 177 (3):1091–1110, 2009. [39] Elvira Papaleo, David G Cornwell, and Nicholas Rawlinson. Seismic tomography of the north anatolian fault: New insights into structural heterogeneity along a continental strike-slip fault. Geophysical Research Letters, 44(5):2186–2193, 2017. [40] Jason P Chang, Sjoerd AL de Ridder, and Biondo L Biondi. Highfrequency rayleigh-wave tomography using traffic noise from long beach, californiatomography using traffic noise. Geophysics, 81(2):B43–B53, 2016. 125 [41] Thomas L Wright. Structural geology and tectonic evolution of the los angeles basin, california: Chapter 3: Part 1. 1991. [42] Dan Hollis, John McBride, David Good, Nicholas Arndt, Florent Brenguier, and Gerrit Olivier. Use of ambient noise surface wave tomography in mineral resource exploration and evaluation. In SEG Technical Program Expanded Abstracts 2018, pages 1937–1940. Society of Exploration Geophysicists, 2018. 126 Chapter 7 A NOVEL FRAMEWORK FOR INVESTIGATING MID-OCEANIC ENVIRONMENTS. This chapter was adapted from: J. C. Castellanos, Z. Zhan, and W. Wu (2020). “Absolute centroid location of submarine earthquakes from 3D waveform modeling of water reverberations.” In: Journal of Geophysical Research 125.5. doi: 10.1029/2019JB018941. 7.1 Abstract Oceanic transform faults (OTFs) represent an attractive tectonic environment to investigate how slip is accommodated within the crust. However, as most of these fault systems grow in the deep ocean, where few local seismic observations are available, characterizing their earthquake behavior is complicated and remains a formidable challenge. Here, we present a novel approach for retrieving precise centroid locations of submarine earthquakes that is based on the modeling of water phases in teleseismic records. Using a hybrid method for simulating far-field body waves with 3D source-side structures, we demonstrate that the scattered energy generated by the continuous bounces of an earthquake’s P-wave trapped in the ocean is modelable and carries information about its source location. As a case study, we use a realistic bathymetry model of the Gofar transform fault on the East Pacific Rise and simulate the seismic wavefield at US Array stations for 4 of its moderate-sized (Mw5.0+) earthquakes. Our modeling results show that water phases are sensitive to a ∼5 km change in the earthquake’s horizontal location and that a remarkable agreement between observed and synthetic water phases is achieved when the location of an event is close to its true one. We then relocate 3 of these events by systematically computing their water phases in candidate locations until a satisfactory waveform fit is achieved. This analysis technique paves a new route for studying earthquake source properties in isolated marine environments and serves as a means to investigate the seismic behavior of OTFs on a global scale. 127 7.2 Introduction Understanding the influence of different tectonic parameters on the earthquake rupture process is key to anticipating fault behavior and assessing seismic hazard. Over the past few decades, significant advancements in this field have been made by studying and monitoring large continental faults such as the San Andreas Fault in California or the Alpine Fault in New Zealand [e.g., 1–5]. However, owing to their obscure origin and complex geologic evolution, the formulation of a mechanical model that can accurately predict seismic behavior along these fault systems has remained one of the most outstanding problems in seismology. Oceanic transform faults (OTFs) on the other hand are among the simplest tectonic environments and, therefore, represent a more advantageous location to study how slip is accommodated within the crust [6–8]. Different from their continental analogues, OTFs have a simple geometric configuration, a close to homogeneous composition and a smooth thermal structure [9–11]. This relative simplicity leads to a somewhat predictable seismic behavior that is less dependent on the geologic history of the plate boundary, allowing us to dissect the different thermal-mechanical factors that control earthquake rupture along transform faults [12]. However, despite the substantial attention that OTFs have received, there are still many of their seismological aspects that are not well understood, especially on a global scale. A general view of the global distribution of oceanic earthquakes reveals a clear difference in seismic behavior between OTFs offsetting slow and fast spreading ridges, with small-to-moderate size earthquakes being common in slow-spreading centers and typically absent in fast-spreading centers [13, 14]. This commonality between ridge systems suggests that the thermal structure of the lithosphere controls the seismicity along OTFs [15–17]. However, recent observations suggest that temperature cannot be the only control on seismicity since even the largest OTF earthquakes rupture only small segments of the thermally expected seismogenic zone [18, 19]. This unusual rupture style results in a dramatic slip deficit in which the majority of the lateral motion on OTFs is accommodated aseismically by mechanisms such as slow creeping events, and, on a smaller portion, seismically on fully coupled patches that are thought to develop above the 600 ◦ C isotherm [7, 18, 20, 21]. To explain the nature of the slow transients, laboratory experiments suggest the presence of large amounts serpentinite along the faults, which under shallow crustal conditions have the potential to result in aseismic fault creep [e.g., 22, 23]. 128 This argument implies that there must exist some level of interaction between seawater and cold mantle rocks that is allowing serpentinization to occur, and that fluid circulation along the fault may, therefore, be as important of a control on OTF seismicity as temperature [e.g., 7]. An unresolved issue, though, is the mechanism that determines the size and location of the coupled patches. Mapped surface traces of historical continental strike-slip earthquakes indicate that geometrical discontinuities in the form of fault step-overs and fault bends hold the potential to act as physical barriers to rupture propagation [24, 25]. Conversely, high-resolution studies along a few OTFs suggest that along-strike variations in fault frictional properties dictate whether or not a rupture can propagate through a particular segment of the fault [7, 26, 27]. Recent studies have also suggested that OTF seismicity has a high level of predictability over different spatial and temporal scales. On short time scales, OTFs show a clear increase in foreshock activity just hours before any mainshock occurs [6]. These foreshock sequences are believed to result from the increased stressing rates caused by preceding slow transients [28, 29]. On long time scales, OTFs exhibit a quasi-periodic seismic behavior with interevent times that are dependent on their slip rates [30]. Furthermore, OTF mainshocks appear to be rupturing on the same overlapping fault patches, which in turn seem to be acting in phase with each other [30]. This type of rupture style suggests that there must exist some type of along-strike variations in the rheology of OTFs that is allowing the synchronization of large earthquakes [31]. However, it is still obscure if both seismic and aseismic slip can occur on the same fault patches or if seismic slip persistently occurs on the coupled patches over many earthquake cycles, while the rest of the fault always slips aseismically [18]. Elucidating any of the enigmas listed above will tremendously improve our understanding of the behavior of OTFs, which could then be extended to their more complicated and hazardous continental counterparts. Such a task requires a detailed characterization of their seismic behavior on a global scale. However, owing to their remoteness, OTF earthquakes are generally located with an uncertainty that is a few times larger than their actual rupture dimensions, thus limiting our capability of probing their fine-scale structure and investigating their spatial variations in seismic properties [26, 30]. This obstacle has resulted in substantial efforts towards developing methodologies 129 that will allow us to accurately constrain the location of mid-oceanic earthquakes in general [e.g., 32]. However, the majority of these techniques can only solve for relative locations and require a small number of events located with insignificant uncertainty, which in most cases are not available [30, 33, 34]. Attempts to determine the absolute locations of OTF earthquakes have been made by associating earthquake focal mechanisms with seafloor bathymetry [35]. Nonetheless, as the resolution of bathymetric data continues to improve, it is being revealed that the morphology of OTFs can be complicated to interpret and that constraining earthquake locations by seafloor bathymetry alone is not always possible. More direct attempts to investigate OTF behavior have been made by deploying ocean bottom seismometers (OBS) [e.g., 7, 36, 37]. However, these studies have been limited to a small number of faults around the world due to the high cost of placement, maintenance, and retrieval of the seismic equipment. This limitation calls for a more efficient yet comparably robust scheme to study OTFs worldwide. In this chapter, we present a novel approach for determining accurate centroid locations of oceanic earthquakes at unprecedented resolutions using land-based seismic stations alone. First, we demonstrate that the waveform complexity that is introduced by the continuous water reverberations of any submarine earthquake is modelable and carries valuable information about the absolute location of the rupture within the fault. We then simulate 4 moderate-sized (Mw5.0+) earthquakes along the Gofar transform fault (GTF) on the equatorial East Pacific Rise (EPR) and show that we can improve their reported locations by modeling their water phases in a set of candidate locations. Finally, we test the sensitivity of water phases to bathymetry changes and show that the current resolution of publicly available bathymetric data is sufficient to model them in any oceanic environment in the world. 7.3 Waveform Modeling of Water-Reverberations Teleseismic recordings of earthquakes occurring directly beneath the seafloor show prominent arrivals caused by the reverberations of P-wave energy trapped in the water column [Figure 7.1; 38, 39]. These arrivals, however, do not often appear as the systematic set of reflections that are expected to result from an acoustic wave propagating in a conventional flat seafloor model (i.e. with decreasing amplitude and alternating polarity). In a series of investiga- 130 tions, [40, 41] examined the effects of a planar dipping water-crust interface on teleseismic waveforms and determined that the high variability of their P-wave coda amplitude results from changes in the take-off angle of water phases. Moreover, Okamoto and Miyatake [42] used 2D numerical simulations to demonstrate that the multipathing effects on far-field body waveforms due to realistic seafloor topography are remarkable and cannot be approximated by a planar interface model. Taken together, these studies suggest that most of the water phase complexity derives from the wavefield’s interaction with the seafloor’s irregular structure. To exemplify this connection, Figure 7.2 shows a set of 2D finite-difference seismograms [43] from a strike-slip source calculated with a two-domain (i.e. fluid and solid) velocity model that includes realistic bathymetry with different levels of smoothness. The synthetic waveform coda shows a systematic increase in amplitude and complexity as the ocean bottom roughness increases that results from the intrinsic relationship between the source radiation pattern and the irregular near-source structure. Figure 7.1: Schematic diagram showing the ray paths of the direct P-wave, the depth phases (pP, sP) and the first four water multiples (pwn P, swn P) (left). Synthetic velocity waveforms for a 20 km deep strike-slip earthquake under a homogeneous water layer (right). For simplicity, the synthetic seismogram is shown without water phases (top) and with water phases (bottom). The water multiples alone are presented in the middle. The waveforms were calculated for a station at a distance of 84◦ using the Preliminary Reference Earth Model [PREM; 44] with a 3-km-thick water layer on the source side. Note how the ray parameter of the water phases is the same as the one of the direct P-wave. Even though water phases in realistic geological settings appear to be overwhelmingly complicated, the causality between seafloor geometry and waveform intricacy presents the unique opportunity to associate particular P-wave coda patterns to both earthquake source properties and the Earth’s structure. 131 Figure 7.2: Example illustrating the effect of seafloor roughness on teleseismic waveforms. The synthetic waveforms show a noticeable increase in amplitude and complexity as the ocean bottom roughness increases (A). The dashed green line on top of the elevation profiles depict the event’s epicenter. The high-resolution bathymetry profile used in the simulations is extracted from ETOPO1 global relief model along profile P-P’ in (B). The low-resolution bathymetry profiles are builtt from applying running averages of different lengths to the high-resolution one. Note that we only show the first 1000 km of each profile to accentuate the near-source structure but the models extend all the way to the inland station. Focal mechanism and takeoff angle of the direct P wave are shown in (C). Such a task, however, heavily relies on our ability to model water reverberations adequately which, in turn, depends on the level of accuracy to which we can capture the bathymetry effect on teleseismic waveforms. To this end, we employ a hybrid method for calculating teleseismic waveforms with 3D source side structures [45]. This method uses the Spectral Element Method [SEM; 46, 47] to compute a complex 3D wavefield in the source region and then propagates that wavefield to large epicentral distances using the Direct Solution Method [DSM; 48, 49]. The coupling of these two methods allows us to capture the full 3D effects of complex source side structures while reducing the heavy computational cost of calculating high-frequency waveforms on a global scale. To show the relevance of considering 3D wave propagation effects in oceanic environments, Figure 7.3 shows a direct comparison of waveform modeling results for the 2013 Mw6.1 GTF earthquake using different structural models (1D, 2D and 3D) and the reported Global Centroid Moment Tensor solution [GCMT; 50, 51]. From this comparison, four critical points should be noted. First, the majority of the P-wave coda of a submarine earthquake is composed 132 of water reverberations that cannot be modeled by synthetics computed for a 1D Earth. Second, for a shallow strike-slip earthquake, a flat seafloor model is insufficient to match the true amplitude of water phases; as the take-off angle of the water multiples is nearly vertical and therefore close to the P and SV nodal planes. Third, an irregular ocean bottom allows larger regions of the focal sphere to be sampled by different rays and causes a larger energy distribution in the P-wave coda. However, 2D modeling is insufficient to match the later part of the waveforms at a desirable level. Fourth, the inclusion of the 3D structure in the epicentral region at an adequate resolution allows us to match the full complexity of the observed waveforms and is, therefore, required to model water reverberations. Note that an even better waveform fit might not have been achieved due to the inherent location error of the GCMT catalog or the assumption of a 1D crust and mantle structure. 7.4 Earthquake Relocation by Waveform Modeling of Water Reverberations: A Case Study in the Gofar Transform Fault Once the bathymetry effects on teleseismic P-waveforms are properly accounted for, water phases should be modeled with high levels of accuracy. Any significant waveform disagreement between observed and synthetic water reverberation phases can, therefore, be attributed to an error on the earthquake’s location, focal mechanism, and/or to the presence of any near-source heterogeneity that is not included in our velocity model. However, given our current knowledge of plate-motion direction and the structural simplicity of OTFs, it is less likely that the two latter factors are the main causes of poor waveform fitting. Moreover, because the impedance contrast is largest at the fluid-solidus interface, the teleseismic P-wave coda of these types of events will be dominated by acoustic waves scattered at the ocean bottom rather than by other type of complex waves scattered at the fault’s and lithosphere’s irregular structure. Receiver-side structures can also introduce some complexity into the observed waveforms, but here we base our analyses on data from dense seismic networks to minimize their effect. Under this set of assumptions, and for a given faulting geometry, we can then generate a library of water phase templates along and around any mid-oceanic fault system to relocate earthquakes whose current catalog location results in poorly modeled coda waveforms. To demonstrate the effectiveness of this approach, we use the fast-slipping (∼ 14cmyr−1 ) GTF as a test case and simulate 4 of its Mw5.0+ events at all available US Array 133 Figure 7.3: Waveform modeling results for the 2013 Mw6.1 Gofar earthquake using different structural models (see Figure 7.4 for the event location). Black traces are observed data and red traces are synthetics (0.01-0.2 Hz). Station network and name is shown in panel (A) only but is the same for all 4 panels. Synthetics in (A) are computed using a 1D Earth (Model 1). Synthetics in (B) are computed using a 2D Earth with a flat ocean bottom on the source side (Model 2). Synthetics in (C) are computed using a 2D Earth with realistic topography (Model 3). Synthetics in (D) are computed using the hybrid method with 3D source side structures (Model 4). The average waveform misfit of each set of synthetics is shown in the upper part of the panels. Examples of the velocity models (for the 1D and 2D simulations) and the bathymetry mesh used in the 3D simulations are shown on the right. The elevation data for Model 4 is extracted from the Global Multi-Resolution Topography synthesis [GMRT; 52]. 134 stations, including the 2008 Mw6.0 Gofar earthquake which was captured by a nearby OBS array (Figure 7.4). The location of the 2008 Mw6.0 event was previously determined relative to its Mw5.1 aftershock and assumed that the aftershock’s centroid location is equal to its own epicenter location [7]. The earthquake depth, on the other hand, is not well-constrained but suggested to be above the 600◦ C isotherm, which given the warm thermal structure of the GTF is ∼5 km for the center of the transform fault [19]. Using the reported GCMT moment tensor solutions and centroid locations, our initial modeling results show that a remarkable waveform fit is achieved for the 2008 Mw6.0 earthquake (most-likely due to its proximity to its revised and probably true location; Figure 7.5C) whereas the water reverberation phases for the remainding 3 earthquakes are poorly matched by the synthetic seismograms (Figure 7.5A,B,D). The large waveform discrepancy between the observed and synthetic traces for the 2002 Mw5.1, 2014 Mw5.5, and 2016 Mw5.7 events is therefore likely to be caused by severe earthquake mislocation. Before relocalizing the poorly-modeled events, it is an interesting exercise to analyze the regional bathymetric effect of the GTF on the teleseismic wavefield. To do so, we simulate two orthogonal lines of equi-spaced strike-slip sources placed at a constant 5-km depth along and across the G3 segment of the fault (locations 7-27 in Figure 7.4) all the way to a 207-station array uniformly spread across the United States. We then use the full synthetic dataset to analyze the wavefield’s sensitivity to source location and perform a source-side beamforming to detect the direction and apparent wave speed at which the multiple arrivals in the P-wave coda are traveling outward from the source array (Figure 7.6). Our analysis reveals that, for the 0.01-0.2 Hz frequency band, water phases are sensitive to a ∼5 km change in the earthquake’s horizontal location and that a large portion of the water multiples are composed of arrivals that are generated from multiple azimuths. However, our most relevant finding is that the fault topographic relief appears to behave as a strong scatterer, as is confirmed by the large concentration of slow arrivals that are coincident with the fault’s strike (Figure 7.6C). This observation supports our hypothesis that water phases can be used to retrieve valuable information about the earthquake’s centroid location as the timing and amplitude of these scattered waves can be used to constrain the event’s position along and across the fault. 135 Figure 7.4: Topographic map of the Gofar transform fault system. G1, G2, G3 mark the different segments of the fault that are separated by intratransform spreading centers. Beachballs show the focal mechanisms and centroid locations from the Harvard GCMT catalogue. The circled numbers represent candidate locations for the 2002 Mw5.1, 2007 Mw6.1, and 2016 Mw5.7 earthquakes. The average intersource spacing between neighboring candidate locations is 5 km. The small inset map is a zoom of the 2016 Mw 5.7 event epicentral region. The dashed purple line marks the extension of a potential sub-parallel fault (from surface bathymetry). The red circled numbers mark the centroid locations determined by this study. Elevation data is extracted from the Global Multi-Resolution Topography synthesis [GMRT; 52]. Main-Fault Events In the simplest of cases, and for earthquakes whose initial reported location is the vicinity of prominent bathymetric features, we can correlate the water reverberation signals generated by several equi-spaced strike-slip sources along any obvious candidate locations until an optimal solution is found. For this case, we limit the depth of all simulated sources to 5 km due to the warm temperature of the GTF, and separate all test locations by 5 km on account of the frequency band of analysis. An example of such a search process is presented in Figure 7.7, which shows the water phase variations along the G2 segment of the fault and their comparison with the 2002 Mw5.1 seismic recordings. Apart from exhibiting significant water phase changes as the source location is shifted along the fault, our computed seismograms are, at a particular location, able to fit almost every single arrival that is present in the observed P-wave codas. 136 Figure 7.5: Waveform modeling results for the 2002 Mw5.1 (A), 2007 Mw6.1 (B), 2008 Mw6.0 (C), and 2016 Mw6.1 (D) Gofar earthquakes using the reported GCMT solutions and a 5-km source depth. Black traces are observed data and red traces are synthetics (0.01-0.2 Hz). Station name and network are shown above each trace pair. By analyzing the goodness-of-fit between the observed and synthetic traces in a total of 99 US Array stations via a shape-fitting criterion, which is defined as follows: [Max(RXY )]2 Misfit = 1 − P 2 P 2 , Xi Yi (7.1) where Xi is the observed waveform, Yi is the synthetic waveform, and RXY is the cross-correlation of the traces, we are able to decrease the average waveform misfit by 22% by taking location 1 from Figure 7.4 as the earthquake’s true centroid location. The difference between our solution and the reported GCMT for this event is over 19 km, which is approximately twice the average rupture length of an Mw5.0 earthquake. Our preferred location displaces the current catalog location away from the extensional pull-apart basin that is between the G1 and G2 fault segments, which is tectonically more consistent with the 137 Figure 7.6: Record section of 3D synthetic waveforms generated for one teleseismic station using source locations 7-27 in Figure 7.4 (A). All waveforms have been band-passed between 0.01 to 0.2 Hz. The blue shaded region delimits the traces for the fault-perpendicular sources whereas the yellow shaded region delimits the traces for the fault-parallel sources. The mean azimuth and distance to the receiver is of 0.70◦ and 51.76◦ , respectively. Beamforming output for the 3 time windows that are enclosed by the red dashed lines in the record section (B). Slow and prominent arrivals can be observed being generated from multiple azimuths. Collection of all maximums in every beamforming output for the full station synthetic array (C). The size of the markers represents the beam power whereas the color represents the apparent wave speed. The presence of slow arrivals concentrated in the azimuth along the fault’s strike suggests that the fault itself is behaving as a strong scatterer. expected faulting geometry of an oceanic strike-slip event. It is important to note that our best-fitting location is at the western end of the source array, and that the true hypocenter of the event may still be further west from our proposed site. However, given the small room for improvement in terms of waveform fit, it is highly unlikely that this location will deviate substantially from our current solution. Low-SNR Events Given the moderate magnitudes of OTF events, as well as their large distance from any seismic station, one difficulty that commonly arises when analyzing their teleseismic waveforms is the low signal-to-noise ratio (SNR) of the earthquake recordings. This complication can bias the interpretation of wa- 138 Figure 7.7: Waveform modeling results for the 2002 Mw5.1 using source locations 1-6 in Figure 7.4. Black traces are observed data and red traces are synthetics (0.01-0.2 Hz). Station name and network are shown above each trace pair in the first panel only. The green contoured panel shows the waveform comparison of the observed and synthetic traces for the location that results in the best fit (location 1) whereas the the blue contoured panel shows the waveform comparison for the location that is closest to the GCMT solution (location 5). ter phases and limit their use in relocating events. Nevertheless, with the availability of dense seismic networks, it is possible to apply beamforming to enhance the coherent signals buried within the traces and suppress the incoherent noise [53]. Here, given that the ray-parameter of the pwn P and swn P phases is substantially close to the one of the direct P-wave arrival (Figure 7.1), we can apply a systematic delay-and-stack using local apparent P-wave slownesses at different groups of stations to extract coherent sequences of water phase arrivals. This process increases the overall SNR of the seismic array at the elevated cost of resolution and, therefore, should only be applied when the seismic traces are uninterpretable. An example of a set of records that can benefit from this SNR enhancement is the one from the 2007 Mw6.1 event. To beamform the signals, we first divide the entire US Array into multiple smallaperture 20-station subarrays such that, for each individual seismic station, we find the optimal stacking velocity using the 19 closest stations and generate a single-beamed trace (Figure 7.8A). To ensure waveform coherence among stations, and maintain the incident plane wave approximation for teleseismic arrivals [53], we only apply beamforming if the radius of a given subarray is smaller than 400 km. Although the choice of stacking parameters is somewhat 139 subjective, one should attempt to use the least amount of stations to allow the subtle waveform changes of the water phases as a function of distance and azimuth to be preserved and, at the same time, still be below the alias limit of the subarray [54, 55]. The outcome of implementing this type of spatio-temporal beamforming results in a remarkable increase in the SNR of the traces such that coherent water phases can now easily be observed (Figure 7.8BC). Figure 7.8: Illustration exemplifying the beamforming process applied to the 2007 Mw6.1 recordings (A). The map shows the stations used to produce a beamed trace at station CI.ARV (red inverted triangle). The blue inverted triangles represent the 19 closest stations used to form the subarray. The inset plot shows the beamforming output for a time window around P-wave arrival of the 20 traces. The subarray and beam information is shown in the bottom white box. The record section in (B) shows the raw waveforms of the 2007 Mw6.1 event at all available US Array stations (gray inverted triangles in A). The record section in (C) shows the waveforms at the same stations after applying the spatio-temporal beamforming. Both sets of records have been aligned with the P-wave onset time and band-pass filtered between 0.01 and 0.2 Hz. The red dashed lines in (B) and (C) mark the P-wave and PP-wave arrival time. The green arrows in (C) mark the arrival time of two prominent water phase arrivals. Once the records have been beamed, and the water reverberation signals are discernible, we take the same approach as the one in the 2002 Mw5.1 event relocation and place multiple strike-slip sources at 5-km depth along the G3 segment of the fault. The waveform modeling results for the 2007 Mw6.1 event 140 are presented in Figure 7.9. A simple visual inspection of the synthetic waveforms reveals that not all candidate locations are able to replicate the water phase patterns that are observed in the recorded seismograms. In particular, there exist two prominent and persistent water phase arrivals at 35 s and 70 s after the direct P-wave (green arrows in Figure 7.8C and Figure 7.9) that are only generated when the earthquake source is placed within a small section of the fault. Based on the source-side beamforming results, we believe that these arrivals are produced somewhere at the edges of the fault probably as a result of seismic guided waves trapped in the fault’s low topographic relief. For obvious reasons, this type of waveform features is useful to investigate near-source seismic structures and, in this case, allow us to tie down the absolute location of the source to a particular point in space. From a goodness-of-fit analysis in over 460 beamed traces, our preferred solution displaces the earthquake’s location 11.6 km east from the reported GCMT location by taking location 15 in Figure 7.4 as the event’s true centroid location. Different from the other events analyzed in this study, the direct P-wave of the 2007 Mw6.1 earthquake has a double-pulse feature that could easily be interpreted as complexity in the rupture process. However, through our waveform analysis, we find that a single point source is sufficient to explain most of the information present in the seismograms, and that, to fully match the direct P-wave and the entire water phase packet simultaneously, the earthquake depth must be close to 15 km (Figure 7.9). This observation suggests that the secondary pulse in the direct P-wave corresponds to a depth phase reflecting from the bottom of the seafloor rather than the result of a complicated source time function. If the latter is correct, then this would imply that either not all earthquakes along the GTF are thermally controlled, or that the thermal structure of the fault cannot be explained by a simple half-space cooling model. This claim is reasonably consistent with the OBS observations of McGuire et al. [7], which revealed large along-strike variations in the material properties of the GTF such that there might be regions of enhanced fluid circulation that are affecting the thermal state of the fault. Although investigating the mechanics of this particular earthquake is out of the scope of this study, our analysis highlights the potential of water phases to place valuable constraints on both the earthquake rupture and the fault structure from teleseismic observations alone. 141 Figure 7.9: Waveform modeling results for the 2007 Mw6.1 earthquake. Black traces are observed data and red traces are synthetics (0.01-0.2 Hz). The leftmost panel shows the synthetic waveforms of earthquake sources 7-17 in Figure 7.4 for station CI.SWS. The black waveform plotted on top corresponds to the observed beamed trace around the same station. The top green arrows mark the arrival time of the two prominent water phases arrivals. Note how source locations 13-16 are the only ones capable of producing the first prominent water phase arrival (blue-circled numbers). The three rightmost panels correspond to the waveform simulation results for a strike-slip source placed at the GCMT reported location at 5 km depth, at source location 15 at 5 km depth, and at source location 15 but at 15 km. Station name and network are shown above each trace pair in the first panel only. The green contour marks the waveform fit of our preferred solution. Off-Main-Fault Events Off-main-fault events such as the 2016 Mw5.7 earthquake present a greater challenge to relocate since there is no primary surface feature to which we can initially tie the earthquake’s epicenter. This requires a broader exploration of the solution space and can therefore become computationally demanding, but not unfeasible. An example of such a search process is shown in Figure 7.10, where we first place a total of 10 strike-slip sources in a 5 km radius around 142 the 2016 Mw5.7 GCMT catalog location and quantify the waveform similarity at each point. We then use the position with the minimum average misfit as our new starting location and simulate 10 new strike-slip-sources around it. This process is applied iteratively until a well-established waveform misfit minimum is found. By quantifying the goodness-of-fit between observed and synthetic records in over 114 US Array stations, we are able to decrease the average waveform misfit by 11% by taking position 42 from Figure 7.4 as the event’s true centroid location. Our new solution is over 10 km east from the current catalog location and is, once again, able to fit almost every single arrival present in the observed P-wave codas. Interestingly enough, our solution places the 2016 Mw5.7 event in close proximity to the 2013 Mw6.1 earthquake, whose current GCMT location also results in a remarkable waveform fit (Figure 7.3D). Both of these events are in close alignment with a linear surface feature that appears to be composed of a series of subparallel faults that appear to be capable of hosting Mw5.0+ earthquakes (Figure 7.4). It is important to mention, however, that since this type of events are not in the immediate vicinity of a prominent bathymetric feature, the water phase sensitivity to horizontal location drops dramatically (as evidenced by the smooth misfit surface exhibited by the candidate locations in Figure 7.10). This limitation increases the uncertainty with which this type of events are relocated and may, therefore, require higher-frequency simulations to determine their accurate location. 7.5 Discussion The modern-day availability of high-resolution multibeam bathymetric data, well-documented slip rates, and comprehensive thermomechanical models of mid-ocean ridges shed light on some of the most relevant factors that control seismicity in OTFs. However, the lack of well-determined earthquake source parameters along these fault systems hinders the accurate characterization of the degree of influence of different tectonic parameters (e.g. offset length, spreading rates and lithosphere temperature) on the distribution of earthquake slip. A thorough characterization of OTF seismicity at a global scale is, therefore, necessary to place strict constraints on the primary controls on earthquake behavior in oceanic plate boundaries. To evaluate the limits of our analysis technique, and whether it can be applied to any oceanic setting on Earth, we generate a series of 2D synthetic seismograms computed with 143 Figure 7.10: Search space for the 2016 Mw5.7 event optimal location. The circles depict the earthquake candidate locations (see magnified map in Figure 7.4) and the traces show the waveform fit between observed (black) and synthetics (red) at station TA.L02F for 4 different locations (0.01-0.2 Hz). The color of the markers represent the waveform misfit for the same station. The waveform modeling results for this earthquake best-fitting location are shown on the right. Station name and network are shown next to each trace pair. The dashed purple line marks the extension of a potential sub-parallel fault (from surface bathymetry). different levels of bathymetric resolution and quantify the waveform similarity as a function of frequency (Figure 7.11). More in detail, this exercise entails computing numerous synthetic waveforms with the same bathymetry profile but sampled at different resolutions and meshed with variable grid sizes so that the wavefields are calculated at different frequencies. Then, for each frequency, we take the synthetic recording that was simulated with the highest resolution velocity model and use it a "true" waveform to be compared with the lower resolution simulations. Our analysis reveals that significant waveform changes begin to occur when bathymetric features with scales larger than the P-wave wavelength in the water are modified. Now, given that water phases propagate in a medium bounded by an elastic half-space and a free surface, the superpo- 144 sition of opposed reverberations leads to a standing wave with one-quarter of a period. Thus, for an average 3 km deep ocean, the fundamental mode of a sequence of water phases will have a central frequency of 0.125 Hz (assuming a constant water velocity of 1.5 km/s). This result suggests that numerical simulations accurate down to ∼8 s period are required to model the simplest form of water phases in an typical mid-oceanic setting, and that bathymetric features of a size comparable to the ocean depth need to be adequately represented in our velocity models. Current publicly available bathymetry models constructed from both satellite data and ship soundings have a 1-arc minute base global resolution and can go as high as 15-arc seconds in regions where more modern ship sounding is available [56, 57]. This type of global resolution technically enable us to accurately model ∼6-s water reverberations in any OTF system in the world. To this point, our methodology has only been applied to relocate events occurring in the equatorial EPR, which has a high spreading rate and, hence, a relatively smooth seafloor geometry. In such a tectonic setting, and because of the frequency band of analysis, our numerical simulations provide us with enough resolution to capture waveform changes resulting from a 5-km location difference. However, in a marine environment with a more chaotic seafloor, such as the mid-Indian ridge, the waveform changes that result from a 5-km source shift are likely to be more notable. As a result, the centroid locations obtained from modeling the water reverberations of earthquakes occurring in slow-spreading ridges should be of higher resolution, and the uncertainties associated with those solutions should be lower. To illustrate this concept, Figure 7.12 shows a comparison between the recordings of two pairs of similarsized transform earthquakes occurring in two different marine environments. Both pairs of events have the same intersource distances, as determined from the cross-correlation of their Rayleigh wave packets [34], and, for each pair, their waveforms are shown for the same seismic stations. From this analysis, it is easily observed that, despite the fact that both pairs have almost the same intersource distances (∼25 km), the water phases from the events occurring in the fast-spreading EPR are more similar to each other than the water phases from the events occurring in the slow-spreading Aden Ridge. Consequently, the location resolution that one may achieve by modeling the water phases of the events in the Aden Ridge is technically higher than the one that is achievable in the EPR. This particularity causes that both the sensitivity and the 145 Figure 7.11: Heatmap showing the waveform sensitivity to bathymetry resolution as a function of frequency. The thick lines depict the P-wave wavelength in the water (λ=1 equals to one wavelength assuming a constant water velocity of 1.5 km/s). The white dashed lines mark the 1-arc minute base resolution of current global bathymetry models, the resolution limit of the resampled SRTM30 global bathymetry model [56], and the the resolution limit of the resampled SRTM15+ global bathymetry model [57]. The thick black arrow marks the bathymetry resolution for the Gofar transform fault (61 m) [52]. The topographic profile used in these calculations is cut from the EPR to the center of the United states (from -4.75◦ /-105.40◦ to 41.38◦ /-119.17◦ lat/lon). The misfit is quantified via the shape-fitting criteria (Equation 7.1) in which the record generated in the highest resolution run is taken as the "true" waveform. Synthetics are calculated for a station at a distance of 48◦ using a 15-km-deep strike-slip source. uncertainties of our technique to vary between regions and, therefore, a careful evaluation of the wavefield’s sensitivity to location changes should always be performed before relocalizing submarine events in different study areas. Here, we have also restricted the analysis of water phases to mid-oceanic settings due to their relatively simple velocity structure and lack of sediments, which have been shown to have a direct effect on the amplitude of teleseismic water reverberations [41]. Nonetheless, our method of analysis can also be applied to extract source characteristics of earthquakes occurring in subduction zones where offshore seismic data has been collected, and high-resolution 3D 146 Figure 7.12: Waveform comparison of two similar-sized pairs of earthquakes occurring in a fast spreading oceanic environment (top) and a slow spreading environment (bottom). Both pairs of events are separated by the same distances, as determined by the time shifts of their fundamental-model Rayleigh waves for a full azimuth range. Event names, and focal mechanisms are shown on the left size of each panel. Station name, distance and azimuth are shown above each pair of waveforms. Distances and azimuths are approximate and calculated from the center of the fault. Fault name and slip rate are shown in the top part of each panel together with the estimated intersource distance. The normalized cross-correlation coefficient of each pair of waveforms is shown in blue next to the traces. Note that, despite having the same intersource distances, the water phases for the events occurring in the fast spreading ridge are more similar than the ones occurring in the slow spreading ridge. All waveforms have been band-passed between 0.01 to 0.2 Hz. For each panel, the focal mechanism and waveforms are color-associated. 147 velocity structure models are available [58]. The analysis of subduction-related water phases may help us to constrain important source characteristics such as the up-dip extent of the earthquake rupture and whether any breakage of the seafloor occurred, which would generate unusually large amplitude water reverberations [59, 60]. Understanding the water phase signature of complicated rupture processes in subduction zones has important implications for probabilistic tsunami hazard assessments and may one day be used to identifying tsunamigenic events in real time. Although so far we have solely focused in improving earthquake locations along the Gofar transform fault, we have shown that our technique can be optimally applied to investigate the seismic behavior of different mid-oceanic environments. An analysis as the one presented here but at a global scale should allow us to answer some of the fundamental open questions regarding the mechanics of OTFs and earthquake faulting in general. For instance, the distribution of seismically active patches along multiple oceanic fault strands around the globe will provide us with sufficient observational information to resolve whether structural discontinuities affect the size and location of the seismogenic zone, or if along-strike variations in fault frictional properties are more likely to limit rupture propagation. Moreover, with an accurate global catalog that spans over a few earthquake cycles, we will be able to determine where the transition from aseismic to seismic takes place along a given fault and whether the frictional properties of OTF segments evolve temporally between stable and unstable slip. Understanding the relationship between seismic and aseismic processes is essential to hazard analyses given that major continental fault systems have been observed to accommodate part of their cumulative strain by aseismic slip [e.g., 61, 62] and because slow transients have been observed to occur minutes before devastating earthquakes [e.g., 63]. The formulation of a reliable global OTF seismic catalog will also serve as an important foundation for numerous future investigations that go well-beyond earthquake source studies. For example, once accurate absolute centroid locations are available, it will become possible to use wavefield methods to probe the fine-scale structure of OTFs and explore the physical conditions within oceanic seismogenic zones. Consequently, with better constraints on the architecture of OTFs, we will be able to investigate how these fault systems serve as pathways for the hydration of oceanic plates and how they interact with the underlying mantle [e.g., 64]. This last topic is of great interest not only to communities 148 studying subduction processes but also to biogeochemical groups investigating the chemical reactions associated with aqueous alteration of peridotite to serpentine minerals, as the interaction between seawater and the cold mantle lithosphere is an important mechanism in the formation of hydrogen, methane and small organic molecules [65, 66]. 7.6 Conclusion We showed that the water reverberation signals of submarine earthquakes are both interpretable and modelable, and can be used to constrain the absolute location of the rupture within a given fault. By computing synthetic water phases using a hybrid method for modeling body wave arrivals with 3D sourceside structures, we relocated 3 OTF events with a ∼5-km uncertainty along the GTF in the EPR using land-based stations alone. The solutions we obtain are within the area of the average location error of the GCMT catalog for OTF earthquakes (∼20 km) and are consistent with the tectonics of the region. We also demonstrated that current global bathymetry maps have enough resolution to allow our methodology to be applied in any oceanic setting in the world. Our approach of using water reverberation phases as a fingerprint of earthquake location is a first step towards a reliable global OTF earthquake catalog and opens the door to a new and unexplored way of studying OTFs in general. References [1] CH Scholz, Max Wyss, and SW Smith. Seismic and aseismic slip on the san andreas fault. Journal of geophysical research, 74(8):2049–2069, 1969. [2] JC Savage and RO Burford. Geodetic determination of relative plate motion in central california. Journal of Geophysical Research, 78(5):832– 845, 1973. [3] Christopher H Scholz, John MW Rynn, Robert W Weed, and Cliff Frohlich. Detailed seismicity of the alpine fault zone and fiordland region, new zealand. Geological Society of America Bulletin, 84(10):3297–3316, 1973. [4] R Jolivet, M Simons, PS Agram, Z Duputel, and Z-K Shen. Aseismic slip and seismogenic coupling along the central san andreas fault. Geophysical Research Letters, 42(2):297–306, 2015. 149 [5] Ian J Hamling, Sigrún Hreinsdóttir, Kate Clark, John Elliott, Cunren Liang, Eric Fielding, Nicola Litchfield, Pilar Villamor, Laura Wallace, Tim J Wright, et al. Complex multifault rupture during the 2016 mw 7.8 kaikōura earthquake, new zealand. Science, 356(6334):eaam7194, 2017. [6] Jeffrey J McGuire, Margaret S Boettcher, and Thomas H Jordan. Foreshock sequences and short-term earthquake predictability on east pacific rise transform faults. Nature, 434(7032):457, 2005. [7] Jeffrey J McGuire, John A Collins, Pierre Gouédard, Emily Roland, Dan Lizarralde, Margaret S Boettcher, Mark D Behn, and Robert D Van Der Hilst. Variations in earthquake rupture properties along the gofar transform fault, east pacific rise. Nature Geoscience, 5(5):336, 2012. [8] B Froment, Jeffrey J McGuire, RD Hilst, P Gouédard, Emily C Roland, H Zhang, and John A Collins. Imaging along-strike variations in mechanical properties of the gofar transform fault, east pacific rise. Journal of Geophysical Research: Solid Earth, 119(9):7175–7194, 2014. [9] J Tuzo Wilson. Transform faults, oceanic ridges, and magnetic anomalies southwest of vancouver island. Science, 150(3695):482–485, 1965. [10] Paul J Fox and David G Gallo. A tectonic model for ridge-transform-ridge plate boundaries: Implications for the structure of oceanic lithosphere. Tectonophysics, 104(3-4):205–242, 1984. [11] Donald Turcotte and Gerald Schubert. Geodynamics. Cambridge University Press, 2014. [12] Margaret S Boettcher. Slip on ridge transform faults: insights from earthquakes and laboratory experiments. PhD thesis, Massachusetts Institute of Technology, 2005. [13] Douglas A Wiens and Seth Stein. Age dependence of oceanic intraplate seismicity and implications for lithospheric evolution. Journal of Geophysical Research: Solid Earth, 88(B8):6455–6468, 1983. [14] DR Bohnenstiehl, RP Dziak, J Steele, S Thorpe, and K Turekian. Midocean ridge seismicity. Encyclopedia of Ocean Sciences, page 15, 2008. [15] Rachel E Abercrombie and Göran Ekström. Earthquake slip on oceanic transform faults. Nature, 410(6824):74, 2001. [16] DV Rundquist and PO Sobolev. Seismicity of mid-oceanic ridges and its geodynamic implications: a review. Earth-Science Reviews, 58(1-2): 143–161, 2002. [17] Margaret S Boettcher, Greg Hirth, and Brian Evans. Olivine friction at the base of oceanic seismogenic zones. Journal of Geophysical Research: Solid Earth, 112(B1), 2007. 150 [18] Margaret S Boettcher and TH Jordan. Earthquake scaling relations for mid-ocean ridge transform faults. Journal of Geophysical Research: Solid Earth, 109(B12), 2004. [19] Emily Roland, Mark D Behn, and Greg Hirth. Thermal-mechanical behavior of oceanic transform faults: Implications for the spatial distribution of seismicity. Geochemistry, Geophysics, Geosystems, 11(7), 2010. [20] Peter Bird, Yan Y Kagan, David D Jackson, S Stein, and JT Freymueller. Plate tectonics and earthquake potential of spreading ridges and oceanic transform faults. 2002. [21] Lynn R Sykes and Göran Ekström. Earthquakes along eltanin transform system, se pacific ocean: Fault segments characterized by strong and poor seismic coupling and implications for long-term earthquake prediction. Geophysical Journal International, 188(2):421–434, 2012. [22] Linda A Reinen, John D Weeks, and Terry E Tullis. The frictional behavior of lizardite and antigorite serpentinites: Experiments, constitutive models, and implications for natural faults. Pure and Applied Geophysics, 143(1-3):317–358, 1994. [23] Linda A Reinen. Slip styles in a spring-slider model with a laboratoryderived constitutive law for serpentinite. Geophysical Research Letters, 27(14):2037–2040, 2000. [24] Stefan B Nielsen and Leon Knopoff. The equivalent strength of geometrical barriers to earthquakes. Journal of Geophysical Research: Solid Earth, 103(B5):9953–9965, 1998. [25] Steven G Wesnousky. Predicting the endpoints of earthquake ruptures. Nature, 444(7117):358, 2006. [26] Monica Wolfson-Schwehr, Margaret S Boettcher, Jeffrey J McGuire, and John A Collins. The relationship between seismicity and fault structure on the discovery transform fault, east pacific rise. Geochemistry, Geophysics, Geosystems, 15(9):3698–3712, 2014. [27] Vaclav Matej Kuna, J Nabelek, and J Braunmiller. Mode of seismic slip at an ocean transform fault: Results from the blanco transform fault obs experiment. In AGU Fall Meeting Abstracts, 2017. [28] Rowena B Lohman and Jeffrey J McGuire. Earthquake swarms driven by aseismic creep in the salton trough, california. Journal of Geophysical Research: Solid Earth, 112(B4), 2007. [29] Emily Roland and Jeffrey J McGuire. Earthquake swarms on transform faults. Geophysical Journal International, 178(3):1677–1690, 2009. 151 [30] Jeffrey J McGuire. Seismic cycles and earthquake predictability on east pacific rise transform faults. Bulletin of the Seismological Society of America, 98(3):1067–1084, 2008. [31] JC Lynch, R Bürgmann, MA Richards, and RM Ferencz. When faults communicate: Viscoelastic coupling and earthquake clustering in a simple two-fault system. Geophysical research letters, 30(6), 2003. [32] Thomas H Jordan and Keith A Sverdrup. Teleseismic location techniques and their application to earthquake clusters in the south-central pacific. Bulletin of the Seismological Society of America, 71(4):1105–1130, 1981. [33] Felix Waldhauser and William L Ellsworth. A double-difference earthquake location algorithm: Method and application to the northern hayward fault, california. Bulletin of the Seismological Society of America, 90(6):1353–1368, 2000. [34] K Michael Cleveland and Charles J Ammon. Precise relative earthquake location using surface waves. Journal of Geophysical Research: Solid Earth, 118(6):2893–2904, 2013. [35] Jianfeng Pan, Michael Antolik, and Adam M Dziewonski. Locations of mid-oceanic earthquakes constrained by seafloor bathymetry. Journal of Geophysical Research: Solid Earth, 107(B11), 2002. [36] Hugh Rowlett. Seismicity at intersections of spreading centers and transform faults. Journal of Geophysical Research: Solid Earth, 86(B5):3815– 3820, 1981. [37] G Horning, Robert A Sohn, J Pablo Canales, and Robert A Dunn. Local seismicity of the rainbow massif on the mid-atlantic ridge. Journal of Geophysical Research: Solid Earth, 123(2):1615–1630, 2018. [38] Robert B Herrmann. Focal depth determination from the signal character of long-period p waves. Bulletin of the Seismological Society of America, 66(4):1221–1232, 1976. [39] Steven N Ward. Ringing p waves and submarine faulting. Journal of Geophysical Research: Solid Earth, 84(B6):3057–3062, 1979. [40] Douglas A Wiens. Effects of near source bathymetry on teleseismic p waveforms. Geophysical Research Letters, 14(7):761–764, 1987. [41] Douglas A Wiens. Bathymetric effects on body waveforms from shallow subduction zone earthquakes and application to seismic processes in the kurile trench. Journal of Geophysical Research: Solid Earth, 94(B3):2955– 2972, 1989. 152 [42] Taro Okamoto and Takashi Miyatake. Effects of near source seafloor topography on long-period teleseismic p waveforms. Geophysical Research Letters, 16(11):1309–1312, 1989. [43] Dunzhu Li, Don Helmberger, Robert W Clayton, and Daoyuan Sun. Global synthetic seismograms using a 2-d finite-difference method. Geophysical Journal International, 197(2):1166–1183, 2014. [44] Adam M Dziewonski and Don L Anderson. Preliminary reference earth model. Physics of the earth and planetary interiors, 25(4):297–356, 1981. [45] Wenbo Wu, Sidao Ni, Zhongwen Zhan, and Shengji Wei. An sem-dsm three-dimensional hybrid method for modelling teleseismic waves with complicated source-side structures. Geophysical Journal International, 215(1):133–154, 2018. [46] Dimitri Komatitsch and Jeroen Tromp. Introduction to the spectral element method for three-dimensional seismic wave propagation. Geophysical journal international, 139(3):806–822, 1999. [47] Dimitri Komatitsch and Jeroen Tromp. Spectral-element simulations of global seismic wave propagation—ii. three-dimensional models, oceans, rotation and self-gravitation. Geophysical Journal International, 150(1): 303–318, 2002. [48] Robert J Geller and Takao Ohminato. Computation of synthetic seismograms and their partial derivatives for heterogeneous media with arbitrary natural boundary conditions using the direct solution method. Geophysical Journal International, 116(2):421–446, 1994. [49] Robert J Geller and Nozomu Takeuchi. A new method for computing highly accurate dsm synthetic seismograms. Geophysical Journal International, 123(2):449–470, 1995. [50] AM Dziewonski, T-A Chou, and JH Woodhouse. Determination of earthquake source parameters from waveform data for studies of global and regional seismicity. Journal of Geophysical Research: Solid Earth, 86 (B4):2825–2852, 1981. [51] Göran Ekström, Meredith Nettles, and AM Dziewoński. The global cmt project 2004–2010: Centroid-moment tensors for 13,017 earthquakes. Physics of the Earth and Planetary Interiors, 200:1–9, 2012. [52] William BF Ryan, Suzanne M Carbotte, Justin O Coplan, Suzanne O’Hara, Andrew Melkonian, Robert Arko, Rose Anne Weissel, Vicki Ferrini, Andrew Goodwillie, Frank Nitsche, et al. Global multi-resolution topography synthesis. Geochemistry, Geophysics, Geosystems, 10(3), 2009. 153 [53] Sebastian Rost and Christine Thomas. Array seismology: Methods and applications. Reviews of geophysics, 40(3), 2002. [54] Philippe Roux, Bruce D Cornuelle, WA Kuperman, and WS Hodgkiss. The structure of raylike arrivals in a shallow-water waveguide. The Journal of the Acoustical Society of America, 124(6):3430–3439, 2008. [55] Florent Brenguier, Philippe Kowalski, N Ackerley, N Nakata, P Boué, M Campillo, E Larose, S Rambaud, C Pequegnat, T Lecocq, et al. Toward 4d noise-based seismic probing of volcanoes: Perspectives from a largen experiment on piton de la fournaise volcano. Seismological Research Letters, 87(1):15–25, 2015. [56] JJ Becker, DT Sandwell, WHF Smith, J Braud, B Binder, Jl Depner, D Fabre, J Factor, S Ingalls, SH Kim, et al. Global bathymetry and elevation data at 30 arc seconds resolution: Srtm30_plus. Marine Geodesy, 32(4):355–371, 2009. [57] B Tozer, DT Sandwell, WHF Smith, C Olson, JR Beale, and P Wessel. Global bathymetry and topography at 15 arc sec: Srtm15+. Earth and Space Science, 6(10):1847–1864, 2019. [58] Yunyi Qian, Shengji Wei, Wenbo Wu, Hongyu Zeng, Aurélie CoudurierCurveur, and Sidao Ni. Teleseismic waveform complexities caused by near trench structures and their impacts on earthquake source study: application to the 2015 illapel aftershocks (central chile). Journal of Geophysical Research: Solid Earth, 124(1):870–889, 2019. [59] Thorne Lay and Andrea Rhode. Evaluating the updip extent of large megathrust ruptures using pcoda levels. Geophysical Research Letters, 46 (10):5198–5206, 2019. [60] Thorne Lay, Chengli Liu, and Hiroo Kanamori. Enhancing tsunami warning using p wave coda. Journal of Geophysical Research: Solid Earth, 2019. [61] William P Irwin and Ivan Barnes. Effect of geologic structure and metamorphic fluids on seismic behavior of the san andreas fault system in central and northern california. Geology, 3(12):713–716, 1975. [62] Frédérique Rolandone, R Bürgmann, DC Agnew, IA Johanson, DC Templeton, MA d’Alessio, SJ Titus, C DeMets, and B Tikoff. Aseismic slip and fault-normal strain along the central creeping section of the san andreas fault. Geophysical Research Letters, 35(14), 2008. [63] Hiroo Kanamori and John J Cipar. Focal process of the great chilean earthquake may 22, 1960. Physics of the Earth and Planetary Interiors, 9(2):128–136, 1974. 154 [64] Caroline M Eakin, Catherine A Rychert, and Nicholas Harmon. The role of oceanic transform faults in seafloor spreading: A global perspective from seismic anisotropy. Journal of Geophysical Research: Solid Earth, 123(2):1736–1751, 2018. [65] TJG Francis. Serpentinization faults and their role in the tectonics of slow spreading ridges. Journal of Geophysical Research: Solid Earth, 86 (B12):11616–11622, 1981. [66] Mathilde Cannat, Fabrice Fontaine, and Javier Escartin. Serpentinization and associated hydrogen and methane fluxes at slow spreading ridges. In Diversity of hydrothermal systems on slow spreading ocean ridges, volume 188, pages 241–264. American Geophysical Union Washington, 2010. 155 Chapter 8 CONCLUSIONS AND FUTURE DIRECTION In this thesis, I presented a set of independent case studies where I used wavefield methods to infer useful properties of the Earth’s structure and dynamics. Chapter 2 represents a case where I use ambient noise surface wave analysis in its most classical way to place constraints on an enigmatic feature of the Mexico subduction system: the transition anatomy of the Cocos slab and its relation with the abrupt termination of the Trans-Mexican Volcanic Belt at its eastern end. This work relied solely on travel-time measurements of noise-derived surface waves to illuminate the large-scale lithospheric structure of south-central Mexico. I also used the velocity discrepancy between orthogonally polarized surface waves to detect the presence of coherently oriented structures within the lithosphere such as fluid melt-filled cracks, lenses of partial melts, and other intrusive bodies such as dykes and sill complexes. The joint analysis of this set of observations pointed to the fact that the transition from flat to steeper subduction is more likely to be accommodated by a slab tear rather than a sharp flexure. With this model, I suggested that the difference in the rollback rates of the subducting plate introduces a suction force that increases the strength of the toroidal flow through the tear and drags the mantle materials under the central portion of the Cocos slab laterally toward southern Mexico. Consequently, any material that is flowing through the tear is unable to reach the surface directly at the site of the slab window, explaining the abrupt termination of the volcanic arc. Chapters 3 and 4 described a set of studies in which I investigate the role that upper mantle buoyancy anomalies play in determining the behavior of the continental crust. In both chapters, I used noise-derived surface waves to build crustal anisotropy models of three tectonically intriguing regions: the Wallowa Mountains in northeastern Oregon, the Rocky Mountain complex near Wyoming, and the Sierra Nevada in California. From the analysis of these three regions, I was able to make a few key observations: (i) In the western US, crustal seismic anisotropy is spatially coherent but unrelated to mantle anisotropy. (ii) There are two main mechanisms that underlie the creation of anisotropy in the lower crust: poloidal flow driven by mantle loading 156 and toroidal flow due to plate interaction. (iii) A mantle lithospheric load depresses the Moho, driving in crustal flow, thickening the crust, and creating an isostasy without topographic expression. (iv) Crustal flow due to mantle vertical loading is especially effective at straining large volumes of lower crust. In general, these results provide important insights into the interplay between mantle processes and crustal dynamics, and allows some broad predictions about continental crust during orogenies. Chapters 5 and 6 explored the use of oil-industry instrumentation and seismic interferometry methods to map small-scale elastic properties of the shallow crust at Long Beach, California. In Chapter 5 I focused on the surface wave portion of the Green’s function to produce a shear wave velocity model of the top kilometer of the crust. In the process of doing so, I developed an automatic neighborhood-based cross-correlation method for phase arrival picking that allows us to deal with complicated waveforms and make reliable phase velocity measurements across multiple frequencies. With the remarkable resolution that was allowed with this type of instrumentation, I illuminated several wellknown geologic features of the area such as as the Newport-Inglewood Fault, the Silverado Aquifer, and the San Gabriel River. I also presented evidence that speaks to the possibility of using seismic anisotropy observations to detect small-scale changes on the regional stress regime and illuminate buried faults. Lastly, I used numerical methods to show that the presence of subsurface smallscale structures have a clear impact on the intensity of the expected shaking, and can cause ground-motion acceleration to change by several factors over a sub-kilometer horizontal scale. In Chapter 6, on the other hand, I focused on the body wave portion of the Green’s function and use array-processing tools to produce a compressional wave velocity model of the top two kilometers of the crust. The results of this work confirmed the reliability of ambient noise body waves to characterize the seismic structure and showed that, in spite of the processing challenges of retrieving them, they can provide unique constraints that are otherwise unattainable with current passive imaging methods. Different from other chapters, Chapter 7 introduced a novel method of analysis that holds great potential for investigating both the seismic behavior and the physical condition of remote mid-oceanic environments. This work used numerical modeling tools to calculate synthetic 3-D wavefields of submarine earthquakes and, with them, associate complicated coda patterns on teleseis- 157 mic signals to specific source properties, particularly location. While the generalization of this approach is still at an early stage, the capability of sparsely investigating oceanic transform seismicity around the world will unlock the door to a broad range of new geophysical research. For instance, with a few single-event location estimations on a given fault, relative relocation methods can be applied to obtain an accurate characterization of its full seismic behavior. Once accurate earthquake catalogs are available, it will not only be possible to study the seismo-tectonics of mid-oceanic environments, but it will also become possible to apply wavefield methods (like the ones presented in this thesis) to explore the physical condition within oceanic seismogenic zones.