SEISMIC STUDIES OF THE ALASKA AND TONGA SUBDUCTION ZONES WITH BODY- WAVE TOMOGRAPHY By Fan Wang A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Geological Sciences – Doctor of Philosophy 2024 ABSTRACT This dissertation focuses on the seismic constraints on the Alaska Peninsula section of the Alaska subduction zone and the central part of the Tonga subduction zone, leveraging the latest onshore and offshore seismic station deployments. I utilized multiple seismic body-wave tomography methods to build 3D P-wave (VP), S-wave (VS), and compressional to shear wave ratio (VP/VS) models at different scales, providing comprehensive insights into key geological processes in subduction zones. Chapter 1 provides an overview of the geologic background and outlines the motivations for investigating the subduction zone structures at the Alaska Peninsula and Tonga. This chapter highlights the distinct regional tectonics and earthquake distributions that characterize these areas, setting the stage for the detailed analyses in the following chapters. Chapter 2 focuses on along-strike variations of the megathrust slip behaviors and structural heterogeneities at the Alaska Peninsula. I present 3D VP and VP/VS models of the Alaska Peninsula to 50 km depth, derived from body-wave double-difference tomography using a newly acquired arrival-time dataset from the Alaska Amphibious Seismic Experiment (AACSE). This chapter also introduces a new slip deficit model obtained from geodetic modeling that is compared with seismic tomography images. The inferred fluid content from the VP/VS model reveals distinct differences at the plate interface and in the overriding plate across the Shumagin, Chignik, Chirikof, and Kodiak segments. High fluid content is inferred at both the plate interface and in the overriding plate at the Shumagin segment, where a low slip deficit is observed. In contrast, low fluid content is inferred at the plate interface and in the overriding plate at the Chignik and Chirikof segments, which display moderate-to-high slip deficit. In the Kodiak segment characterized by the highest slip deficit, the plate interface is fluid-rich while the overriding plate is relatively dry. The notable anti-correlation of the fluid content in the overriding plate and slip deficit distribution highlights the critical role of the overriding plate in controlling along-strike variations of slip behaviors in the Alaska Peninsula. Chapter 3 focuses on the deep part of the Alaska subduction from ~50 to 700 km depths, including the mantle wedge and the subducted slab in the upper mantle. I introduce two sets of 3D velocity models of the Alaska Peninsula at different scales: large-scale VP and VS models extending to 700 km, and local-scale VP and VP/VS models focusing on structures within the top 200 km. The large-scale models reveals slab gaps in the mantle transition zone (MTZ) beneath the Chirikof segment, while the local-scale models exhibit strong along-strike variations, correlating with seismic activity and sub-arc melting controlled by slab dehydration. The multi-scale velocity models presented in Chapter 3, together with those from Chapter 2, provide a comprehensive investigation of structural variations in the Alaska Peninsula. They offer valuable insights into a variety of subduction-related processes, ranging from shallow ones including megathrust and interseismic slip to deeper processes such as sub-arc melting and slab deformation in the MTZ. Chapter 4 investigates the deep slab morphology and mantle wedge structures in the Tonga subduction zone. Two sets of VP and VS models were developed using arrival-times from manual picks and machine learning detection, showing notable consistency. The subducted Pacific Plate flattens out in the MTZ north of 19° S, while it penetrates the 660-km discontinuity south of 19° S. A slab tear at ~19° S implies strong slab deformation in the MTZ. Sub-slab low-velocity anomalies below 300 km suggest the presence of subducted asthenosphere or the influence of the Samoan plume from the north. The deep slab morphology and mantle structures revealed by body- wave tomography provide crucial insights into the complex tectonic processes in the Tonga subduction zone, where significant back-arc spreading and plume interactions are involved. This thesis is dedicated to my parents. Thank you for always being supportive of me iv ACKNOWLEDGEMENTS First and foremost, I would like to express my deepest gratitude to my Ph.D. advisor, Dr. S. Shawn Wei, for your unwavering support, inspiration, and belief in me throughout my doctoral journey. Thank you for guiding me to develop scientific ideas, improve my writing and presentation skills, and grow into a mature independent scientist. I have learned so much from you, particularly in how you carefully deal with details and create solid plans to achieve goals step by step. Your consistent belief in my abilities has meant the world to me, even during the challenging moments of my Ph.D. I am grateful for your help in securing a postdoctoral position and for your care about my future plans when I felt uncertain. Your encouragement to join the research cruise and your offer to keep me in the group during a time of great anxiety about my future have provided immense comfort and reassurance. I want to thank my Ph.D. committee: Dr. Allen McNamara, Dr. Jeffrey Freymueller, Dr. Susannah Dorfman, for your invaluable support and guidance throughout my research and in shaping my future career. From you, I learned how to develop a broad perspective on my research, the importance of gaining knowledge from other disciplines, and how to communicate complex concepts in a way that is accessible to those outside the field. I want to thank my lab colleagues and friends—Connor Drooff, Dongdong Tian, Guanyuan Shuai, Jiaqi Li, Jiaxin Zhang, Katarina Vance, Keyi Chen, Luis Martinetti, Mingda Lv, Ruxin Shi, Shuqi Wang, Xinming Deng, Yaqi Jie, Yuqian Zhang, Yurong Zhang, Ziyi Xi, Zechao Zhuo, Zhuohui Xiao, and Zhuoran Zhang. Your friendship and support have brought so much joy and meaning to my life. I feel truly fortunate to have had you by my side when I needed it. I want to deeply thank my parents for your constant support and for always standing by me, especially during times when I felt uncertain about my career. Your love and belief in me have v given me the strength and foundation to live a fulfilling life. vi TABLE OF CONTENTS CHAPTER 1: INTRODUCTION ................................................................................................... 1 REFERENCES ........................................................................................................................... 6 CHAPTER 2: FLUIDS CONTROL ALONG-STRIKE VARIATIONS IN THE ALASKA MEGATHRUST SLIP ............................................................................................ 9 REFERENCES ......................................................................................................................... 62 APPENDIX ............................................................................................................................... 69 CHAPTER 3: SLAB MORPHOLOGY, DEHYDRATION, AND SUB-ARC MELTING BENEATH THE ALASKA PENINSULA REVEALED BY BODY-WAVE TOMOGRAPHY .................................................................................................. 75 REFERENCES ....................................................................................................................... 127 CHAPTER 4: DEEP SLAB MORPHOLOGY CONSTRAINED BY BODY-WAVE TOMOGRAPHY IN THE TONGA SUBDUCTION ZONE ............................. 135 REFERENCES ....................................................................................................................... 170 vii CHAPTER 1: INTRODUCTION The subduction of negatively buoyant oceanic lithosphere is a fundamental process in plate tectonics, where dense oceanic plates sink into the Earth's mantle along convergent boundaries. This process is usually accompanied by abundant seismic activity, including some of the most powerful ones on Earth, known as megathrust earthquakes. These earthquakes typically occur at the plate interface and can generate devastating tsunamis, posing severe risks to coastal regions. In addition to these shallow (0–50 km) megathrust earthquakes, subduction zones also host most intermediate-depth (70-300 km) and deep-focus earthquakes (>300 km) (Hasegawa and Nakajima, 2017; Zhan, 2020). The mechanisms of these deeper earthquakes remain enigmatic, as they occur well beyond the brittle-ductile transition and exhibit significant regional variability (Zhan, 2020). Beyond their seismological significance, subduction zones are crucial for material exchange between the Earth's surface and interior, playing a key role in geochemical and geodynamic processes. The subduction of the oceanic lithosphere facilitates the transport of surface materials into the mantle, which is eventually recycled back to the surface through arc volcanism (Stern, 2002; van Keken et al., 2011). A thorough investigation of multiple subduction zones is essential for advancing our understanding of earthquakes, plate tectonics, and the Earth's evolution over geological time. Seismological studies in subduction zones have utilized various methods to investigate subduction zone related processes by analyzing different seismic phases in seismograms, such as body-wave and surface-wave phases (van der Hilst et al., 1997; Martin-Short et al., 2018). Multi- scale seismic tomography with various seismic phases is essential for imaging structural heterogeneities in subduction zones. Tomographic studies using body waves offer higher lateral resolution, while surface-wave tomography provides improved depth resolution, though its 1 imaging depth is limited by low-frequency observations. Subducted slabs typically appear as high- velocity anomalies, contrasting with low-velocity structures associated with mantle plumes and sub-arc magmatism. Global tomography has reported diverse large-scale slab morphologies in the mantle transition zone, showing some slabs flattening at its base while others penetrating into the lower mantle (Fukao & Obayashi, 2013). This global seismic tomography captures velocity structures over hundreds of kilometers but often lacks finer details. In comparison, regional seismic tomography achieves a more detailed view of slabs and mantle wedge structures in subduction zones, revealing features such as slab tears, buckling, and sub-arc magmatic paths with higher resolution (Martin-Short et al., 2018; Yang & Gao, 2020; Zhang et al., 2019). However, the resolution in subduction zones remains constrained by sparse seismic station coverage as major areas lie beneath the ocean where instrumentation is limited. Outstanding scientific questions are emerging related to detailed subduction zone processes as seismic instrumentation and imaging techniques advance. At the seismogenic zone depth (<50 km), devastating megathrust earthquakes happen at the plate interface that displays substantial variations globally (Bassett & Watts, 2015; Wallace, 2020; Wang & Bilek, 2014). Structural heterogeneities at the plate interface and in overriding plate—such as seamounts, sediment thickness, fluid content, and lithology—are implied to influence megathrust and slipping behaviors (Bassett et al., 2016; Shillington et al., 2015; Wei et al., 2021; Wang et al., 2024). However, the primary factors governing slab slipping behaviors remain debated. At upper mantle depths, the distribution of intermediate-depth earthquakes shows significant variability across global subduction zones, often attributed to thermal control (Syracuse et al., 2010; Wei et al., 2017). Yet, detailed regional studies reveal complexities in intermediate-depth seismicity, slab dehydration, and associated arc volcanism, which cannot be solely explained by slab thermal structures (Wei et 2 al., 2021). In some subduction zones, along-strike slab dehydration shows large variations, influenced by factors such as sediment, pre-existing seafloor fabrics, and the lithology of the overriding plate (Martin-Short et al., 2018; Wei et al., 2021). These factors likely interact in complex ways, highlighting the need for further investigation. At the mantle transition zone depths, morphological features such as slab gaps/tears suggest interactions between subducting slabs and mantle flow, sometimes complicated by mantle plumes rising from greater depths (Chang et al., 2016; Tang et al., 2014). The dynamics of slabs within a convecting mantle remain uncertain, especially when regional tectonics are unclear. Therefore, detailed 3D imaging of slab and mantle structures is essential for advancing our understanding of megathrust slip, slab dehydration, sub- arc volcanism, and mantle dynamics. This dissertation focuses on seismological studies of the Alaska Peninsula section of the Alaska subduction zone in the North Pacific and the Tonga subduction zone in the Southwest Pacific—two regions with distinct seismicity distributions and subduction dynamics. The structures of these two subduction zones were not well-constrained until recently due to a lack of offshore seismic stations. The Pacific Plate subducts northward beneath the Alaska Peninsula at approximately 65 mm/yr, whereas the subduction rate in the Tonga subduction zone is significantly higher, ranging from 160 to 240 mm/yr (Bevis et al., 1995; Bird, 2003). This faster subduction rate and older plate age in the Tonga trench make the Tonga slab thermally distinct from the Alaskan slab (Syracuse et al., 2010). The Tonga subduction zone hosts around two-thirds of all the global deep-focus earthquakes, while seismic activity extends to only about 200 km depth beneath the Alaska Peninsula. Previous seismic tomography studies reveal a complex deep slab morphology in the MTZ beneath Tonga, while beneath the Alaska Peninsula, the slab exhibits a relatively simple low-angle dip above 200 km depth (Nayak et al., 2020; van der Hilst, 1995). 3 Additionally, a double seismic zone, double layers of seismic activity, extends to ~120 km depth in the Alaska Peninsula and reaches much deeper at ~300 km depth in Tonga, implying a thermal controlling mechanism on intermediate-depth earthquakes (Florez and Prieto, 2019; Wei et al., 2021; Wei et al., 2017; Yamasaki and Seno, 2003). Besides, the back-arc and regional tectonics exhibit distinct differences in these two subduction zones. Multiple active back-arc spreading centers have been developed in the Lau back-arc basin west of the Tofua volcanic arc and the Tonga Trench (Martinez and Taylor, 2002), whereas the Alaska Peninsula has not experienced widespread back-arc spreading recently (Tibaldi and Bonali, 2017). The subduction of the Pacific Plate beneath the Alaska Peninsula has remained stable since ~40 Myr ago (Madsen et al., 2006), while the subduction of the Tonga slab is affected by the ancient slab from the fossil Vitiaz trench and Samona plume to the north (Chen and Brudzinski, 2001; Smith et al., 2001). The distinct differences in the regional tectonics in the Alaska and Tonga subduction zones make them ideal case studies to investigate a variety of factors that influence the slab structure, sub-arc/back-arc melting, and earthquakes. In this dissertation, I utilize multiple body-wave seismic tomography techniques to image the slab and mantle wedge structures of the Alaska Peninsula section of the Alaska subduction zone and the central part of the Tonga subduction zone at different scales. High-resolution imaging of subduction zone structures benefits from offshore data recorded by recent ocean bottom seismographs (OBSs), as the majority of the regions are undersea. Besides, recent developments in the double-difference tomography technique (tomoDD), utilized in this work, provide improved constraints on structural heterogeneity near earthquake source regions (Zhang & Thurber, 2003). This tomoDD technique leverages event pairs in close spatial proximity—a common feature in subduction zones. Because closely spaced event pairs share similar ray paths outside the 4 earthquake source region, the effects of along-path velocity heterogeneities in the far field will cancel out if taking differential travel time between these event pairs to the same station. Therefore, the use of differential travel times between closely spaced event pairs can increase resolution in near-source regions. Furthermore, I use an improved version of tomoDD (teleTomoDD; Pesicek et al., 2014) that incorporates teleseismic arrival times to enhance the resolution of the slab structure down to the base of the mantle transition zone, as presented in Chapters 3 and 4. This large-scale velocity model obtained from teleTomoDD serves as a good 3D initial model for further regional body-wave tomography in Chapters 2 and 3, where a refined version of tomoDD is applied to invert for VP and VP/VS structures (Guo et al., 2018). This improved tomoDD method for VP and VP/VS structures minimizes biases from different ray paths of P and S waves by selecting only closely aligned P and S wave ray paths, thereby improving inverted VP/VS structures. Multiple seismic velocity models developed in this dissertation shed light upon key subduction processes ranging from megathrust behavior at seismogenic depths to deep slab morphology at the mantle transition zone. The newly constructed high-resolution VP, VS, and VP/VS models of the Alaska Peninsula reveal strong structural controls of the overriding plate on megathrust slip behaviors at shallow depths (Chapter 2). These models also highlight the impact of slab dehydration on sub-arc melting in the forearc and the presence of slab gaps in MTZ (Chapter 3). In the Tonga subduction zone, complex slab morphology and deep sub-slab low-velocity observed in our VP and VS models suggest a strong influence of regional tectonics on slab structure (Chapter 4). The 3D velocity models presented in this dissertation will serve as reference models for sophisticated analyses of megathrust behavior, arc volcanism, and slab deformation in the MTZ. This work will enhance our understanding of seismic hazards, deep earthquake mechanisms, and subduction dynamics. 5 REFERENCES Bassett, D., Sandwell, D.T., Fialko, Y., Watts, A.B., 2016. Upper-plate controls on co-seismic slip in the 2011 magnitude 9.0 Tohoku-oki earthquake. Nature 531, 92–96. https://doi.org/10.1038/nature16945 Bassett, D., Watts, A.B., 2015. Gravity anomalies, crustal structure, and seismicity at subduction zones: 2. Interrelationships between fore-arc structure and seismogenic behavior. Geochemistry, Geophys. Geosystems 16, 1541–1576. https://doi.org/10.1002/2014gc005685 Bevis, M.G., Taylor, F.W., Schutz, B.E., Recy, J., Isacks, B.L., Helu, S., Singh, R., Kendrick, E., Stowell, J., Taylor, B., Calmant, S., 1995. Geodetic observations of very rapid convergence and back-arc extension at the Tonga arc. Nature 374 (6519), 249-251. DOI: https://doi.org/10.1038/374249a0 Bird, P., 2003. An updated digital model of plate boundaries. Geochemistry, Geophys. Geosystems 4 (3). DOI: https://doi.org/10.1029/2001gc000252 Chang, S.-J., Ferreira, A.M.G., Faccenda, M., 2016. Upper- and mid-mantle interaction between the Samoan plume and the Tonga–Kermadec slabs. Nat. Commun. 7, 10799. https://doi.org/10.1038/ncomms10799 Chen, W.P., Brudzinski, M.R., 2001. Evidence for a large-scale remnant of subducted lithosphere beneath Fiji. Science 292 (5526), 2475-2479. DOI: https://doi.org/10.1126/science.292.5526.2475 Florez, M.A., Prieto, G.A., 2019. Controlling Factors of Seismicity and Geometry in Double Seismic Zones. Geophys. Res. Lett 46 (8), 4174-4181. DOI: https://doi.org/10.1029/2018gl081168 Fukao, Y., Obayashi, M. 2013. Subducted slabs stagnant above, penetrating through, and trapped below the 660 km discontinuity, J. Geophys. Res. Solid Earth, 118(11), 5920-5938, DOI: https://10.1002/2013jb010466. Guo, H., Zhang, H., Froment, B., 2018. Structural control on earthquake behaviors revealed by high-resolution Vp/Vs imaging along the Gofar transform fault, East Pacific Rise. Earth Planet. Sci. Lett. 499, 243–255. https://doi.org/10.1016/j.epsl.2018.07.037 Hasegawa, A., Nakajima, J., 2017. Seismic imaging of slab metamorphism and genesis of intermediate-depth intraslab earthquakes. Phys. Earth Planet. Inter 4 (1). DOI: https://doi.org/10.1186/s40645-017-0126-9 Madsen, J.K., Thorkelson, D.J., Friedman, R.M., Marshall, D.D., 2006. Cenozoic to Recent plate configurations in the Pacific Basin: Ridge subduction and slab window magmatism in western North America. Geosphere 2 (1). DOI: https://doi.org/10.1130/ges00020.1 6 Martinez, F., Taylor, B., 2002. Mantle wedge control on back-arc crustal accretion. Nature 416, 417-420. DOI: https://doi.org/10.1038/416417a Martin-Short, R., Allen, R., Bastow, I.D., Porritt, R.W., Miller, M.S., 2018. Seismic imaging of the Alaska subduction zone: Implications for slab geometry and volcanism. Geochemistry, Geophys. Geosystems. 19, 4541–4560. DOI: https://doi.org/10.1029/2018gc007962 Nayak, A., Eberhart-Phillips, D., Ruppert, N.A., Fang, H., Moore, M.M., Tape, C., Christensen, D.H., Abers, G.A., Thurber, C.H., 2020. 3D seismic velocity models for Alaska from joint tomographic inversion of body-wave and surface-wave data. Seismol. Res. Lett 91 (6), 3106-3119. DOI: https://doi.org/10.1785/0220200214 Pesicek, J.D., Zhang, H., Thurber, C.H., 2014. Multiscale seismic tomography and earthquake relocation incorporating differential time data: Application to the maule subduction zone, Chile. Bull. Seismol. Soc. Am. 104, 1037–1044. https://doi.org/10.1785/0120130121 Shillington, D.J., Bécel, A., Nedimović, M.R., Kuehn, H., Webb, S.C., Abers, G.A., Keranen, K.M., Li, J., Delescluse, M., Mattei-Salicrup, G.A., 2015. Link between plate fabric, hydration and subduction zone seismicity in Alaska. Nat. Geosci. 8, 961–964. https://doi.org/10.1038/ngeo2586 Smith, G.P., Wiens, D.A., Fischer, K.M., Dorman, L.M., Webb, S.C., Hildebrand, J.A., 2001. A complex pattern of mantle flow in the Lau backarc. Science 292 (5517), 713-716. DOI: https://doi.org/10.1126/science.1058763 Stern, R.J., 2002. Subduction Zones. Rev. Geophys 40 (4). DOI: https://doi.org/10.1029/2001rg000108 Syracuse, E.M., van Keken, P.E., Abers, G.A., 2010. The global range of subduction zone thermal models. Physics of the Earth and Planetary Interiors 183 (1-2), 73-90. DOI: https://doi.org/10.1016/j.pepi.2010.02.004 Tang, Y., Obayashi, M., Niu, F., Grand, S.P., Chen, Y.J., Kawakatsu, H., Tanaka, S., Ning, J., Ni, J.F., 2014. Changbaishan volcanism in northeast China linked to subduction-induced mantle upwelling. Nat. Geosci 7, 470–475. https://doi.org/10.1038/ngeo2166 Tibaldi, A., Bonali, F.L., 2017. Intra-arc and back-arc volcano-tectonics: Magma pathways at Holocene Alaska-Aleutian volcanoes. Earth Sci Rev 167, 1-26. DOI: https://doi.org/10.1016/j.earscirev.2017.02.004 van der Hilst, R.D., 1995. Complex morphology of subducted lithosphere in the mantle beneath the Tonga trench. Nature 374, 154-157. DOI: https://doi.org/10.1038/374154a0 van der Hilst, R.D., Widiyantoro, S., Engdahl, E.R., 1997. Evidence for deep mantle circulation 7 from global tomography. Nature 386, 578–584. DOI: https://doi.org/10.1038/386578a0 van Keken, P.E., Hacker, B.R., Syracuse, E.M., Abers, G.A., 2011. Subduction factory: 4. Depth- dependent flux of H2O from subducting slabs worldwide. J. Geophys 116 (B1). DOI: https://doi.org/10.1029/2010jb007922 Wallace, L.M., 2020. Slow slip events in New Zealand. Annu Rev Earth Planet Sci. 48, 175–203. https://doi.org/10.1146/annurev-earth-071719-055104 Wang, K., Bilek, S.L., 2014. Invited review paper: Fault creep caused by subduction of rough seafloor relief. Tectonophysics 610, 1–24. https://doi.org/10.1016/j.tecto.2013.11.024 Wang, F., Wei, S.S., Drooff, C., Elliott, J.L., Freymueller, J.T., Ruppert, N.A., Zhang, H., 2024. Fluids control along-strike variations in the Alaska megathrust slip. Earth Planet. Sci. Lett. 633, 118655. https://doi.org/10.1016/j.epsl.2024.118655 Wei, S.S., Ruprecht, P., Gable, S.L., Huggins, E.G., Ruppert, N., Gao, L., Zhang, H., 2021. Along-strike variations in intermediate-depth seismicity and arc magmatism along the Alaska Peninsula. Earth Planet. Sci. Lett 563, 116878. DOI: https://doi.org/10.1016/j.epsl.2021.116878 Wei, S.S., Wiens, D.A., van Keken, P.E., Cai, C., 2017. Slab temperature controls on the Tonga double seismic zone and slab mantle dehydration. Sci. Adv 3 (1), e1601755. DOI: https://doi.org/10.1126/sciadv.1601755 Yamasaki, T., Seno, T., 2003. Double seismic zone and dehydration embrittlement of the subducting slab. J. Geophys. Res. Solid Earth 108 (B4). DOI: https://doi.org/10.1029/2002jb001918 Yang, X., Gao, H., 2020. Segmentation of the Aleutian-Alaska subduction zone revealed by full- wave ambient noise tomography: Implications for the along-strike variation of volcanism. J. Geophys. Res.: Solid Earth 125. DOI: https://doi.org/10.1029/2020jb019677 Zhan, Z., 2020. Mechanisms and implications of deep earthquakes. Annu Rev Earth Planet Sci 48 (1), 147-174. DOI: https://doi.org/10.1146/annurev-earth-053018-060314 Zhang, H., Thurber, C.H., 2003. Double-difference tomography: The method and its application to the Hayward fault, California. Bull. Seismol. Soc. Am. 93, 1875–1889. https://doi.org/10.1785/0120020190 Zhang, H., Wang, F., Myhill, R., Guo, H., 2019. Slab morphology and deformation beneath Izu- Bonin. Nat Commun 10, 1310. https://doi.org/10.1038/s41467-019-09279-7 8 CHAPTER 2: FLUIDS CONTROL ALONG-STRIKE VARIATIONS IN THE ALASKA MEGATHRUST SLIP An edited version of this chapter was published in Earth and Planetary Science Letters: Wang, F., Wei, S. S., Drooff, C., Elliott, J. L., Freymueller, J. T., Ruppert, N. A., et al. (2024), Fluids control along-strike variations in the Alaska megathrust slip, Earth Planet. Sci. Lett., 633, 118655, doi: 10.1016/j.epsl.2024.118655. 9 Abstract Interseismic and coseismic slip on the subduction zone megathrust show complex along-strike variations and the controlling factors remain debated. Here we image seismic velocity anomalies to infer fluid distribution in the Alaska subduction zone (Alaska Peninsula section) with seismic tomography. The weakly locked Shumagin segment is characterized by abundant fluids on the plate interface and in the overriding plate, whereas the moderately-to-highly locked Chignik and Chirikof segments that hosted M>8 earthquakes are relatively dry. The Kodiak segment, which was ruptured by the 1964 M9.2 Great Alaska earthquake and is presently fully locked near the trench, is characterized by a fluid-rich plate interface but a fluid-poor and inferred low-porosity overriding plate. Multiple large earthquakes have nucleated in fluid-rich regions at the plate interface. Our study highlights the important roles of fluids, particularly in the overriding plate, in controlling interseismic slip deficit and earthquake rupture. We propose that a fluid-rich overriding plate means that there has been a high sustained flux of fluids across the plate interface, which enhances the formation of clay minerals that promote fault creep. 10 2.1. Introduction Almost all great earthquakes (M>8) occur on megathrust faults in subduction zones. Previous observations and simulations focus on the plate interface and subducted plate and suggest that varying roughness, sediment thickness, and fluids at the plate interface may cause different slip behaviors, ranging from aseismic slip to megathrust earthquakes (Audet et al., 2009; Guo et al., 2021; Heise et al., 2017; Kodaira et al., 2004; Lay et al., 2012; Li et al., 2018a; Moreno et al., 2014; Shillington et al., 2015; Sun et al., 2020). Specifically, fluids released from sediment compaction and metamorphic dehydration reactions of hydrous minerals in the subducted slab play a critical role in controlling the stress state and frictional properties along the megathrust (Saffer and Tobin, 2011; Tobin and Saffer, 2009). These fluids can be either trapped at the plate interface or released into the overriding plate. The porosity of the overriding plate thus depends on not only the fluid supply from slab dehydration but also the permeability and evolution history of the materials overlying the megathrust. Although a classic depth-dependent model focusing on fluid conditions at the plate interface (Lay et al., 2012; Saffer and Tobin, 2011) can well explain seismic and aseismic activities in the Nankai subduction zone, more complexities have been observed in other subduction margins (Nishikawa et al., 2019). In recent years, properties of the overriding plate in the forearc region, including its lithology, porosity, and rigidity, have received increasing attention as controlling factors of megathrust slip behavior and the earthquake rupture process (Arnulf et al., 2022; Bassett et al., 2016; Bassett et al., 2014; Chesley et al., 2021; Egbert et al., 2022; Sallares and Ranero, 2019; Shillington et al., 2022). However, questions remain on which factor or factors dominate the interseismic and coseismic slip behaviors. A comprehensive investigation of a single subduction zone with earthquake and slip features that vary along strike is needed to understand the controlling factors of slip 11 behaviors and their influence on large megathrust earthquakes. The Alaska Peninsula section of the Alaska-Aleutian subduction zone displays complex along-strike variations in interseismic slip behavior and crustal and mantle structures, providing an ideal natural laboratory to study slip behaviors on the megathrust (Cordell et al., 2023; Drooff and Freymueller, 2021; Li et al., 2018a; Liu et al., 2022a; Lynner, 2021; Shillington et al., 2015; Wei et al., 2021; Xiao et al., 2021). We identify four along-strike segments from southwest to northeast based on the seismic and geodetic observations in this study (Figure 2.1): Shumagin, Chignik, Chirikof, and Kodiak. Although the region between the Shumagin and Kodiak islands was conventionally referred to as the “Semidi” segment (e.g., Shillington et al., 2022), we divide it into the Chignik and Chirikof segments along a trench-normal boundary crossing the Semidi Islands, because these two segments exhibit distinct characteristics from our seismic tomography images and geodetic modeling. The most recent megathrust earthquakes, the 2020 M7.8 Simeonof and 2021 M8.2 Chignik events, ruptured two close-spaced segments that show distinct differences in interseismic coupling (Davies et al., 1981; Drooff and Freymueller, 2021; Elliott et al., 2022, Xiao et al., 2021) (Figure 2.1), which shows that interseismic coupling alone cannot predict seismic hazard. In this region, the Pacific Plate subducts beneath the North American Plate at a convergence rate of ~6.3 cm/yr in a nearly trench-normal direction (Argus et al., 2010). Interseismic geodetic observations indicate that the subducted slab and the overriding plate are fully locked together outboard of Kodiak Island, whereas the megathrust slip transitions to weak locking outboard of the Shumagin Islands (Drooff and Freymueller, 2021; Xiao et al., 2021). It is worth noting that a region of low or no frictional locking adjacent to a fully frictionally locked region may still exhibit a high interseismic slip deficit due to the stress shadow effect (Lindsey et al., 2021; Zhao et al., 2022). Therefore, we use slip deficit instead of plate locking to describe the 12 geodetic constraints on interseismic slip behavior in the rest of this study. In addition to the 1964 M9.2 Great Alaska earthquake that ruptured the high-slip-deficit Kodiak segment, large megathrust earthquakes have occurred every ~50-75 years in the nearby Chignik segment where a high slip deficit is also observed (Davies et al., 1981; Li et al., 2018a; Ye et al., 2021) (Figure 2.1). The most recent 2021 M8.2 Chignik earthquake mostly ruptures the 20- to 40-km depth of the plate interface within the Chignik segment. In contrast, the Shumagin segment with a low slip deficit did not rupture in a great megathrust earthquake for at least 200 years prior to the 2020 M7.8 Simeonof earthquake (Elliott et al., 2022; Witter et al., 2014; Ye et al., 2021). What causes these along-strike changes despite the uniform convergence rate and direction? One possible answer is that the incoming plate exhibits along-strike changes. The pre-existing seafloor fabric and normal faults in the Pacific Plate have a low angle with the trench to the southwest in the Shumagin segment, while their orientations are nearly perpendicular to the trench to the northeast in the Chirikof and Kodiak segments (Shillington et al., 2015) (Figures 2.1, 2.4). Consequently, more bending faults are re-activated at the outer rise, and thus the incoming plate is likely to be more hydrated in the Shumagin segment than in the neighboring Chignik segment. The Shumagin segment has a thin sediment layer being subducted, whereas thick sediments are present in the Chignik, Chirikof, and Kodiak segments due to the subduction of the Zodiac and Surveyor sedimentary fans (Li et al., 2018a). Given these systematic along- strike changes, Shillington et al. (2015) hypothesized that changes in the hydration state of the incoming lithosphere may control interseismic plate coupling and megathrust earthquake potential in the Alaska Peninsula through the volume of subducted fluids. On the other hand, the changes in the overriding plate crust properties may also influence megathrust earthquakes in this region (Shillington et al., 2022). In addition, a recent electromagnetic study suggests that the 13 weak plate coupling in the Shumagin segment does not solely result from fluids (Cordell et al., 2023). Yet the details of the controlling factors remain unclear due to the lack of 3D images of the megathrust and overriding plate in this region. In this study, we utilize a seismic body-wave double-difference tomography technique to construct a 3D high-resolution velocity structure of the Alaska Peninsula. We also develop a new slip deficit model based on continuous and campaign GNSS measurements, independent from our seismic imaging. Furthermore, we investigate the structural variations of the overriding plate and the plate interface in different segments along the strike, and their relationships with interseismic slip deficit. 2.2. Data and Methods 2.2.1. Seismic travel-time tomography In this study, we image the high-resolution 3D seismic velocity structure of the Alaska Peninsula and the adjacent regions using body-wave travel-time data recorded by the Alaska Amphibious Community Seismic Experiment (AACSE) and the adjacent stations of the EarthScope Transportable Array and the Alaska Earthquake Center (Ruppert et al., 2022) (Figure 2.2). AACSE ocean-bottom seismographs (OBSs) in the outer rise and fore-arc regions provide critical data for high-resolution images of the megathrust, which is mainly offshore. We conduct a body-wave travel-time tomography imaging on multiple scales. First, we develop a large-scale velocity model using regional and teleseismic data. Second, we use this model as a starting model to image the VP and VP/VS structure down to 200 km depth with only regional seismic data. With this strategy, our final results are largely independent of the starting model (Figure 2.3). In the first step, we assemble the P- and S-wave arrival times from May 2018 through 14 August 2019 recorded by onshore and offshore seismic stations of the AACSE, EarthScope USArray, and the Alaska regional network (Figure 2.2c), as well as teleseismic traveltime data of Alaskan earthquakes recorded by global stations and global events recorded by Alaskan stations. We use an improved version of teletomoDD (Pesicek et al., 2014), a teleseismic double- difference tomography package, to invert the regional and teleseismic arrival times for 3-D VP and VS models down to 700 km depth. The starting velocity model is constructed by combining the top 50 km part of a recent regional VS model constrained by ambient noise and teleseismic surface waves (Li et al., submitted) and a 3D global model TX2019slab below 50 km depth. We perform four iterations of joint inversion for velocity and relocation and two iterations of inversion for only relocation. In the second step, we invert for higher-resolution VP and VP/VS models from 0 to 200 km depths using only regional seismic data, using two sets of starting models, and a new version of tomoDD (Guo et al., 2018; Zhang and Thurber, 2003) that can simultaneously determine VP, VS, and VP/VS models using P, S, and S-P arrival times. The first starting VP and VP/VS models directly come from the first step. However, because the VP and VS models from the first step have different resolutions, using them to construct a starting VP/VS model will produce some unrealistic values and potentially bias inversion in the second step. Therefore, for the second set of starting models, the VP model still comes from the first step, whereas the starting VP/VS model is assumed with a uniform value of 1.75. The tomographic inversion includes 9 iterations of joint inversion for velocity and relocation as well as 6 iterations of inversion for only relocation. The VP and VP/VS models are always inverted simultaneously to avoid mapping errors in one model into the other. As shown in Figure 2.3, the output VP/VS models given two different starting models are similar. Therefore, we prefer the inverted VP and VP/VS models with a starting VP/VS 15 value of 1.75 (Figures 2.3a, 2.3c, 2.4-2.6). Details of the data preparation and inversion are described in Appendix Text 2.1. Although our tomography method does not require a slab surface (i.e., the top of the plate interface), we also construct a geometry model of the slab surface for interpretation. This model is based on a model constrained by active-source seismic experiments (Kuehn, 2019) in addition to the Slab2 model (Hayes et al., 2018). We start with the Slab2 model as the slab surface geometry (Hayes et al., 2018) (Figure 2.20a). In the regions with active-source imaging (Kuehn, 2019), we update the slab surface geometry in the Shumagin, Chignik, and Chirikof segments because these results are presumably more accurate than the Slab2 model. The slab surface becomes about 7 km deeper, particularly in the Chignik and Chirikof segments. We then add a 50-km-wide buffer zone to transit from the active-source-imaged region to the background Slab2 model. 2.2.2. Geodetic modelling We estimate the slip deficit along the subduction interface using the TDEFNODE software package (McCaffrey, 2002) following Drooff and Freymueller (2021). This approach predicts velocities at the GPS sites on the Peninsula block by combining the motions due to tectonic block rotation with the elastic deformation associated with the accumulation of interseismic slip deficit. Compared to the slip deficit model of Drooff and Freymueller (2021), we extended the slip deficit model further to the east to include the entire Kodiak Island and updated the slip deficit distribution in the eastern part of the model with more data. We approximate the Aleutian megathrust as a series of planar dislocations in an elastic halfspace (Okada, 1992) and create a grid of nodes with 10 km spacing along the strike and 5 km spacing along the dip, assuming the Slab2 (Hayes et al., 2018) slab surface geometry from 10 to 16 75 km depth. Then we estimate the slip deficit at each grid node and interpolate it between the node grids. Using the updated slab surface model from Section 2.2.1 will result in negligible changes on the estimated slip deficit along the strike and small changes along the dip, because the maximum change in the slab surface depth is only ~7 km. We define block boundaries, angular velocities, and rotational poles according to the published data from the Peninsula (Li et al., 2016), Southern Alaska (Fletcher and Freymueller, 2003), and the Bering (Cross and Freymueller, 2008). The Pacific and North American plate boundaries and rotational poles are defined by Argus et al. (2010). All sites in this study lie on the Peninsula block, and the relative motion of the Peninsula block and Bering Plate is small (Cross and Freymueller, 2008; Elliott and Freymueller, 2020; Li et al., 2016). Therefore, we assume that the Alaska-Aleutian subduction zone is the only source of substantial elastic strain accumulation. We invert for slip deficit using only the horizontal components of site velocities. The vertical site velocity residuals are ~3 mm/yr in our results. Thus, long-wavelength patterns are likely due to the inaccuracies of regional Glacial Isostatic Adjustment models for the collapse of the forebulge of the Laurentide ice sheet as well as inaccurate loading models of the historic glaciers in the Alaska Peninsula (Li and Freymueller, 2018). The velocities of the sites in the Alaska Peninsula and west of Kodiak are taken from Drooff and Freymueller (2021), including their corrections to remove the volcanic inflation signal from Veniaminof volcano. We use velocities from Elliott and Freymueller (2020) for sites in the Kodiak archipelago. They are derived from the same set of daily GPS solutions, and simply represent different spatial subsets. We divide up the subduction interface into discrete segments, using the same segment boundary locations as Drooff and Freymueller (2021) (5-segment model) outboard of the Shumagin and Semidi Islands. We test the effect of an additional coupling boundary inserted 17 near Tugidak Island in the southwestern Kodiak archipelago (6-segment model), corresponding to the boundary between the Chirikof and Kodiak segments (Figure 2.7). We assume a Gaussian distribution along the dip for the slip deficit in each segment. The slip deficit is uniform within each segment, but there is no enforced along-strike smoothing between segments (Li and Freymueller, 2018; Drooff and Freymueller, 2021). The Gaussian model predicts a much lower slip deficit at shallower depths than the model that assumes the slip deficit monotonically decreases with depth (Drooff and Freymueller, 2021; Xiao et al., 2021). Compared to the monotonic model, the Gaussian model leads to a low slip deficit near the trench, which may help explain the lack of a significant tsunami. We are unable to distinguish between the Gaussian and monotonic coupling models based on the onshore data alone, but both are useful in providing endmember scenarios for minimum and maximum extents of the coupling patches. We choose to regularize our coupling model using a Gaussian method in order to provide a minimum extent of the regions of coupling as required to fit the onshore geodetic data. 2.3. Results 2.3.1. Seismic tomography images Our results reveal the seismic velocity structure shallower than 50 km depth, extending from the Alaska Trench to arc volcanoes on the Alaska Peninsula (Figures 2.4-2.6). Compared to a recent body-wave tomography study using a similar dataset (Gou et al., 2022), our results show similar large-scale patterns but provide a much higher resolution. This is because our method is more sensitive to small-scale anomalies in the earthquake source region, and we minimize any systematic biases by simultaneously inverting for VP, VP/VS, and earthquake hypocenters. A gently dipping high-VP structure (~6.6-8.5 km/s) is imaged beneath the slab surface (Figure 2.5). The overriding plate at 20 km depth shows higher VP values (~7.2-7.5 km/s) in the Shumigan 18 segment compared to the segments to the east (Figure 2.4a). But we do not observe dramatic changes in VP on the slab surface (Figure 2.4b). The VP/VS structure shows strong variations along the strike in the overriding plate and on the slab surface (Figures 2.4c, 2.4d, 2.6). We do not observe VP/VS anomalies higher than 2.0, such as was imaged by studies of receiver functions in the Cascadia subduction zone (Audet et al., 2009), possibly because our tomography models have relatively lower spatial resolution. In the Shumagin segment to the southwest, the overriding plate and slab surface show a high VP/VS signature from near the trench to ~40 km depth (Figures 2.4, 2.6). In the Chignik and Chirikof segments to the east, low VP/VS anomalies are the dominant features in the overriding plate and on the slab surface at depths greater than 20 km. At depths shallower than 20 km towards the trench, the slab surface shows a high VP/VS in the Chignik segment but low VP/VS in the Chirikof segment. In the Kodiak segment further to the northeast, the overriding plate displays low VP/VS anomalies, but the plate interface shows high VP/VS values (Figure 2.6). These strong along- strike changes of VP/VS structures in each segment are moderate-scale features with sizes larger than 40 km horizontally and 10 km vertically. 2.3.2. Seismic tomography resolution tests To evaluate the tomography model resolution, we perform 4 synthetic tests with different input models (details in Appendix Text 2.1.2). The first and second tests are checkerboard tests with different sizes of input velocity anomalies to evaluate the spatial resolution (Figures 2.10- 2.12). The third and fourth restoration tests are designed to validating our capability of imaging high VP/VS anomarlies at the plate interface with along-strike variations (Figures 2.13-2.16). The recovered checkerboard models suggest that the slab structure is well imaged at 20–40 km depths, particularly near earthquake sources (Figures 2.10-2.12). In the first checkerboard test 19 (Figure 2.10), we calculate the semblance between the recovered and input models at each node, and use it as a proxy for model resolvability (Zelt, 1998). By visually comparing the input and recovered models, we choose the semblance value of 0.54 as the threshold of model resolvability and use this contour to mask out the low-resolution regions in the final VP and VP/VS velocity models. As shown in Figures 2.10-2.12, we achieve a resolution of 25 and 10 km in the horizontal and vertical directions, respectively, in the unmasked region. Therefore, we focus on discussing large-scale (≥30 km horizontally and ≥10 km vertically) structures in our models because some small-scale (<20 km horizontally) anomalies may reflect model resolution changes along the strike and dip directions due to ununiform distributed stations and earthquakes. The third resolution test shows that the high VP/VS layer at the plate interface, representing hypothetically hydrated slab sediment and crust, can be well recovered at 15–50 km depths throughout the region (Figure 2.13). In addition, the high VP/VS anomaly in the overriding plate, representing a possible fore-arc sedimentary basin, will not be mistakenly mapped to greater depths. Therefore, we are confident with the imaging results near the plate interface. In addition, the similar patterns between the input and recovered models in the fourth resolution test indicate the major along-strike changes are robust in our tomography images (Figures 2.14-2.16). In summary, all resolution tests suggest that our VP/VS model has a horizontal resolution of 25 km in the unmasked regions. Thus, we are confident with all moderate-scale anomalies that are mentioned in Section 2.3.1 and will be discussed in Section 2.4. Our resolution limit is larger than the typical width of the plate interface deformation zone (<10 km) or the scale of fluid-rich layers (<1 km) (Li et al., 2015). Therefore, the VP and VP/VS structures imaged at the slab surface represent the average velocity structure across the plate interface. 20 2.3.3. Plate coupling model We find a ~14% decrease in the total weighted residual sum of square velocities for sites fit in the 6-segment model than in the 5-segment model which does not include the new boundary near Tugidak Island, and assumes the Chirikof and Kodiak segments have the same slip deficit distribution. The best-fit coupling distributions and residual site velocities are shown in Figure 2.7. The main new feature of this model is the higher slip deficit outboard of Kodiak than in the Chignik segment and the high slip deficit outboard of Chirikof Island. In addition, the new slip deficit model has a narrower highly locked patch in the Kodiak segment than the previous model (Drooff and Freymueller, 2021), and this highly locked patch is at a shallower depth compared with the one in the Chignik segment. The nearly 100% slip deficit outboard of Kodiak is similar to that found in models that were optimized for the 1964 rupture zone (i.e., Elliott and Freymueller, 2020; Suito and Freymueller, 2009). We note that our model includes a zone of partial coupling by the trench outboard of Kodiak where the model resolution is minimal. The estimated slip deficit at shallow depth is controlled mainly by the model regularization rather than the data. We observe a distinct trend of trench-parallel residual velocities for sites on Kodiak Island (Figure. 2.7), which we interpret as unaccounted-for localized forearc translation on the order of ~5 mm/yr at most sites. This is a common phenomenon observed in subduction settings globally (Jarrard, 1986; Wang et al., 2003). In the case of Kodiak, this effect was accounted for in Elliott and Freymueller (2020) by adding an upper crustal strike-slip fault aligned approximately with the 25 km depth contour underneath Kodiak Island which bounded two upper plate blocks. We note that the transition from partial coupling to creeping beneath Kodiak in our model is in the same place as the along-strike upper crustal fault included by Elliott and Freymueller (2020). 21 However, it will dramatically increase the model misfit in the Chirikof segment if we add a creeping boundary west of the Kodiak segment in the same location as in Elliott and Freymueller (2020). In order to test whether this crustal sliver affects the estimated slip deficit distribution, we projected the Kodiak site velocities into the trench-normal direction and re-ran the inversion for slip deficit without the trench-parallel motion. We find no difference in the coupling model from such a model run. Therefore we conclude that the effect of a Kodiak crustal sliver on the coupling model is negligible and the misfit caused by the projection will be passed on into the residuals rather than biasing the model parameters. Besides, adding an additional block to the model would require us to adopt the estimated block velocity from Elliott and Freymueller (2020). This additional block velocity was not purely trench parallel and it is provided to achieve a better fit to the GPS data further to the east, which is outside of our study area. Geodetic modeling of this additional block will require to extend the model domain even further east by adding more data to help constrain the block motion, and re-considering block boundaries in a complex region outside of our study area. 2.4. Discussions 2.4.1. Geological implications of VP/VS anomalies and effects of seismic anisotropy and slab geometry In this study, we focus on the VP/VS images while also using the VP model to help with interpretation. We estimate the VP and VP/VS of potential compositions on the Alaska slab surface using a MATLAB toolbox developed by Abers and Hacker (2016), and compare them with the VP and VP/VS values in our tomography models (details in Appendix Text 2.2.1; Figure 2.17). The comparisons suggest that the extremely high VP/VS (>1.83) anomalies on the slab surface indicate the existence of fluids. In addition, we evaluate the effects of the azimuthal anisotropy, 22 titled transversely isotropy, and sub-slab anisotropy in the asthenosphere on the along-strike variations of the VP/VS anomalies, and conclude that these effects are minor and will not change our interpretations (details in Appendix Text 2.2.2; Figure 2.19). It is important to emphasize that our tomography images do not have a sufficient resolution to capture details of fluid-rick layers at the plate interface that are on a scale of hundreds of meters to a few kilometers (Li et al., 2015). The VP/VS values of particular features are likely underestimated due to the limited resolution. Therefore, we do not quantitatively translate the imaged VP/VS values to rock porosity. Instead, we qualitatively interpret high VP/VS anomalies as a proxy for fluids and hydrous minerals in highly deformed rocks and cracks. Besides the geometry model of the slab surface described in Section 2.2.1, we also try other slab geometry models, including the original Slab2 model and a model purely constrained by our relocated seismicity distribution (Figures 2.20a, 2.20b), for interpretation. In addition, we average VP/VS from 10 km above to 10 km below the slab surface of the chosen geometry model (Figures 2.20c, 2.20d, 2.20e). We find no obvious changes in the along-strike VP/VS variations among all the slab surface models (Figure 2.20). Therefore, we conclude that the choice of the slab surface geometry has minor effects on our interpretations, partly because our tomography resolution (~10 km vertically) is larger than the slab surface model discrepancies. 2.4.2. Seismic imaging of fluids in highly deformed rocks The subducted slab is imaged as a gently dipping (< 30°) region of high-VP anomalies (~6.6- 8.5 km/s) beneath the slab surface (Figure 2.5). A thin (< 10 km) relatively lower VP layer (~6.6- 6.9 km/s) overlays the high-VP anomalies (~ 7.0-8.5 km/s) right beneath the slab surface. These layers indicate the slab crust and mantle, respectively. The slab VP structure is in good agreement with previous active-source and receiver-function seismic imaging in this region (Li et al., 23 2018a; Onyango et al., 2022; Shillington et al., 2022; Shillington et al., 2015) (Figure 2.18). The VP/VS structure shows more complexities, with values ranging from 1.65 to 1.90 (Figures 2.4, 2.6). This is because VP is dominantly controlled by bulk composition, temperature, and pressure, whereas VP/VS is highly sensitive to the existence of fluids (Takei, 2002). VP/VS varies with depth on the slab surface, indicating the presence of hydrous minerals and fluids from different sources (Figure 2.4d). Near the trench at 10–20 km depths, we observe high VP/VS anomalies in almost all segments, in agreement with the fluid overpressure observed in other subduction zones from rock experiments and seismic studies (Bassett et al., 2014; Tobin and Saffer, 2009). Furthermore, a recent study using absolute pressure gauge data from the AACSE detected a possible slow-slip event in late 2018 near the trench (He et al., 2023), also implying high fluid content at shallow depths. At greater depths, the distribution of high VP/VS anomaly varies. The fluids at depths shallower than 20 km presumably come from sediment compaction, whereas the sources for fluids at 20–50 km depths remain under debate. Some numerical simulations suggest that the Alaska slab is cold at these depths and that slab dehydration primarily happens below 70 km (Abers et al., 2017), although a recent study demonstrated that the metamorphic dehydration reactions may take place at shallower depths if non-steady-state thermal evolution of the slab is considered (Holt and Condit, 2021). A recent electromagnetic study in the Shumagin segment suggests that slab mantle dehydration may provide substantial fluids in the fore-arc (Cordell et al., 2023). Nevertheless, we attribute the high VP/VS anomalies in the subducted slab to the locally hydrated slab crust and fluids either from in-situ dehydration or from migration along the plate interface up-dip from deeper sources. The high VP/VS anomalies imaged in the overriding plate crust indicate fluids migrating from the plate interface 24 because we do not expect any dehydration reactions in the fore-arc crust. We observe systematic along-strike variations in VP/VS and thus the inferred fluid content (Figure 2.6). We acknowledge that the model resolution changes along strike as the station and earthquake distributions are not uniform. However, we focus on velocity anomalies on a scale of ≥30 km horizontally and ≥10 km vertically, which are well resolved across the entire study region (Figures 2.10-2.16). These anomalies are likely smoothed images of geological features, such as the interplate deformation zone or localized fluid-rick pockets, which are on a scale smaller than our tomography resolution. Thus, the imaged velocity values should not be over- interpreted. The Shumagin segment (represented by cross-sections A-A’ and B-B’ in Figure 2.6) is characterized by continuous high VP/VS anomalies (1.86-1.90) and high fluid content both at the plate interface and in the overriding plate crust. Interestingly, the fluid distribution inferred from a recent 2D electromagnetic study (Cordell et al., 2023) is nearly identical to that inferred from our cross-section B-B’, showing a fluid-rich plate interface at ~30 km depth and a fluid- starved plate interface shallower than 25 km. However, this electromagnetic profile is located near our cross-section A-A’, where we image a fluid-rich plate interface from 30 km to at least 10 km depth. We suspect the discrepancies between the 2D electromagnetic profile with our 3D seismic images result from the along-strike variations within the Shumagin segment that cannot be captured by a 2D profile. Our 3D images show that the megathrust at shallow (<25 km) depths is more hydrated than that revealed by Cordell et al. (2023), although the fluid distribution is not uniform throughout the entire Shumagin segment. In contrast to the Shumagin segment, moderately high VP/VS anomalies (1.80-1.84) are sporadically observed at the plate interface in the Chignik and Chirikof segments (represented by cross-sections C-C’ and D-D’ in Figure 2.6, respectively), whereas low VP/VS anomalies 25 dominate the overriding plate there, implying dry conditions in these segments. In the Kodiak segment (represented by cross-section E-E’ in Figure 2.6), the plate interface has high VP/VS anomalies (1.82-1.86), indicating high fluid content, whereas the overriding plate shows low VP/VS and low inferred porosity. Cross-section F-F’ in Figure 2.6 shows the VP/VS structure along the 30-km depth contour of the slab surface, suggesting that the plate interface and overriding plate crust are more hydrated in the Shumagin segment than the segments to the east. This trend is in good agreement with active-source seismic studies which suggest more water subducted in the southwest than in the northeast (Shillington et al., 2015). The changes along cross-section F- F’ can also explain the along-strike variations in azimuthal anisotropy revealed by teleseismic shear-wave splitting (Lynner, 2021), as the highly hydrated slab at 20–40 km depths in the Shumagin and Kodiak segments produce trench-parallel-fast-direction signals. 2.4.3. Fluids enhance interseismic slip along the megathrust The VP/VS images at the slab surface and in the overriding plate show first-order similarities with our updated slip deficit model, with VP/VS anti-correlating with slip deficit along the dip (Figure 2.6) and along the strike (Figure 2.8), except in the near-trench area where the seismic and geodetic model resolutions are both limited. In the Shumagin segment to the southwest, high VP/VS anomalies at the plate interface and in the overring plate coincide with low slip deficit and high seismic activity, suggesting abundant fluids enhance interseismic creep along the megathrust (Figures 2.4, 2.6, 2.7). In contrast, in the Chignik segment to the northeast, low VP/VS anomalies at the plate interface and in the overriding plate, a higher slip deficit, and low seismic activity imply that the low fluid content inhibits interseismic creep and promotes frictional locking. Further to the northeast, the Chirikof segment, which has the widest region of high slip deficit has moderate-to-low VP/VS, implying a dry plate interface. At the north-eastern end of this 26 region, the Kodiak segment has a high slip deficit and high VP/VS anomalies at the plate interface but low VP/VS in the overring plate. In Figure 2.22, we also compare our VP/VS model with the slip deficit model from Zhao et al. (2022), which models interseismic plate coupling with asperities based on historical and recent rupture zones. This comparison shows a similar anti- correlation between VP/VS and slip deficit. The contrast in the Kodiak segment is surprising, given that many studies suggest that fluids at the plate interface should enhance the interseismic slip along the megathrust (Audet et al., 2009; Guo et al., 2021; Heise et al., 2017; Kodaira et al., 2004; Li et al., 2018a; Moreno et al., 2014). But our observation is not the only case that shows an anti-correlation between VP/VS in the overriding plate and slip deficit. For example, a recent tomography study of the Honshu subduction zone (Zhao et al., 2022) using offshore seismic data shows that VP/VS changes little at the megathrust, but reduces dramatically in the overriding plate from northern Honshu to central Honshu, coinciding with the increase in interseismic slip deficit and the 2011 M9.1 Tohoku-Oki earthquake rupture area. In the Hikurangi subduction zone, Henrys et al. (2020) also image an increase in VP/VS in the forearc overriding plate from the southern North Island to Cook Strait, anti-correlating with the slip deficit change. In the southern Cascadia subduction zone with a high slip deficit, high VP/VS anomalies are confined within the slab crust and low VP/VS anomalies lie above the plate interface (Delph et al., 2018; Guo et al., 2021). Given all observations above, we propose that the along-strike variations in interseismic megathrust slip behaviors at 20-50 km depths are dominantly controlled by the changes in fluid content in the subduction zone, particularly in the overriding plate crust that directly contacts the plate interface at 15–35 km depths (Figure 2.9). We cannot rule out the other factors, such as the plate interface roughness and changes in the slab surface geometry, that may also influence 27 interseismic slip (Cordell et al., 2023; Jiang et al., 2022). However, the fluid content in the overriding plate is the most evident factor that varies along the strike and well anti-correlates with the interseismic slip deficit. Systematic changes in the hydration state of the plate interface are presumably caused by the varying input of the incoming plate and the slab dehydration after subduction. The hydration state of the overriding plate is fundamentally affected by the permeability of the overriding materials and the fluid supply from the plate interface. It is worth noting that some previous studies suggest that changes in the overriding plate lithology, such as the deposit of quartz into the overriding forearc crust, can affect the inter-plate slip behaviors (Audet and Burgmann, 2014). Here at the Alaska Peninsula, two major Mesozoic and Paleogene tectonic boundaries cut across all the boundaries between segments with different slip deficits (Horowitz et al., 1989). The Border Ranges fault system (BRFS) formed when the Chugach terrane was juxtaposed against and beneath the Peninsula terrane to the north starting from the early Jurassic (Plafker et al., 1994), and the Prince William terrane was underthrust by the Chugach terrane along the Contact fault system (CFS) further south (Figure 2.21). Since both the BRFS and CFS are parallel to the trench, we do not expect significant changes in the overriding plate lithology along the strike (cross-section F-F’ is in the same Chugach terrane from Figure 2.21). Therefore, we suspect that the along-strike variations in VP/VS are predominantly influenced by the fluid supply from slab dehydration rather than the overriding plate lithology. We propose that the flux of fluids across the plate interface has a greater influence on the slip behavior than simply the presence of fluids along the interface. The high fluid pore pressure that decreases the effective normal stress is important but not the most critical factor in controlling the megathrust interseismic slip at the Alaska Peninsula. A high fluid flux from the subducting plate through the 28 plate interface would lead to the formation of clay minerals that promote fault creep, through processes of alteration, dissolution, and re-precipitation of mineral grains (Ikari et al., 2009; Katayama et al., 2015). In the Shumagin segment, the incoming plate crust and uppermost mantle are highly hydrated through extensive outer-rise bending faults (Li et al., 2023; Liu et al., 2022b; Shillington et al., 2015) (Figure 2.9). As this hydrated plate subducts, hydrous minerals in the slab crust and mantle break down, releasing water, and triggering abundant intermediate-depth earthquakes (Cordell et al., 2023; Wei et al., 2021). The released fluids migrate up-dip along the strongly deformed plate interface (Kuehn, 2019; Shillington et al., 2022), and hydrate the overriding plate. In addition, a thin, heavily faulted sediment layer subducts in this segment (Li et al., 2018a). Consequently, the combination of more fluids and high deformation at the plate interface, high porosity in the overriding plate, and the lack of sediment facilitate interseismic megathrust slip, resulting in more small earthquakes but a lack of great earthquakes. In contrast, the incoming plate is less faulted and hydrated with more sediments subducted in the Chignik segment (Li et al., 2018a; Li et al., 2023; Shillington et al., 2015), and thus less water is subducted to great depths to trigger intermediate-depth earthquakes (Wei et al., 2021). The low fluid supply from deep dehydration, combined with the weakly deformed plate interface, a thick sediment layer, and a low-porosity overriding plate, result in a lower aseismic fault slip, and thus fewer small earthquakes but several great events in the past. In the Chirikof segment further to the east, the plate interface from the trench to ~30 km depth and the overriding plate crust are so dry that the two plates are fully frictionally locked. In the Kodiak segment at the eastern end, the incoming plate is more hydrated compared to that in the Chirikof segment, because more normal faults that are sub-parallel to the trench are re-activated at the outer rise (Reece et al., 2013) and 29 several seamounts enter into the trench (von Huene et al., 2021). Consequently, more water is subducted to great depths and then released to form a sporadically fluid-rich plate interface, which is consistent with petrological studies of exhumed fore-arc rock (Vrolijk, 1987) and a slow-slip event in 2009 at 20–40 km depths (Holtkamp, 2017). However, the overriding plate crust remains relatively dry, possibly because the less deformed lower crust of Kodiak Island acts as a permeability barrier above the plate interface (Saffer and Tobin, 2011), resulting in a lack of fluid flux across the interface and through the rocks along the interface. Therefore, the interseismic megathrust slip deficit is high in this region. 2.4.4. Complicated roles of fluids in generating megathrust earthquakes Our seismic images and interseismic slip model help explain the diverse mechanisms of historical and recent megathrust earthquakes at the Alaska Peninsula (Davies et al., 1981; Elliott et al., 2022). Note that the velocity anomalies and the inferred high-porosity regions are possibly more concentrated than those appearing in our tomography images. The 2020 M7.8 Simeonof earthquake ruptured down-dip from 25 to 40 km depths and propagated westward within the Shumagin segment (Xiao et al., 2021; Ye et al., 2021) (Figures 2.4d, 2.6). This segment is characterized by a low interseismic slip deficit (~25% of the plate motion rate), extremely high VP/VS anomalies and thus high fluid content at the plate interface and overriding plate, and a lack of great megathrust earthquakes (M > 8) for at least 200 years. However, the Simeonof earthquake rupture initiated near the boundary between the Shumagin and Chignik segments where VP/VS changes from high to low and the slip deficit changes from low to moderate (Figures 2.4d, 2.7). We suspect that shear stress had accumulated due to the stress shadow effect of the adjacent Chignik segment, and high fluid content lowered the plate interface strength to facilitate rupturing. Northeast of the Shumagin segment, the 2021 M8.2 Chignik earthquake 30 nucleated in a high VP/VS region and ruptured into low VP/VS regions in both the down-dip and eastward directions (Figures 2.4d, 2.6). The rupture area occurred from 20 to 40 km depths within the moderately high slip deficit region of the Chignik segment (Elliott et al., 2022; Liu et al., 2022a). We propose that this earthquake started in a fluid-rich region and mainly ruptured the dry plate interface. So the relatively dry, partially locked Chignik segment acted as a barrier for the 2020 Simeonof earthquake rupturing, and then was ruptured by the 2021 Chignik earthquake. These diverse slip behaviors of fault barriers have been seen in numerical simulations (Molina- Ormazabal et al., 2023). Although the rupture area of the 1938 M8.3 Alaska Peninsula earthquake remains under debate, it likely ruptured from near the trench to 30 km depth and overlaps with the eastern part of the 2021 Chignik earthquake rupture area (Freymueller et al., 2021). More than 80 years after this earthquake, our results show low fluid content and high slip deficit in this region. The southwestern end of the 1964 M9.2 Great Alaska earthquake exhibits a high slip deficit and low VP/VS in the overriding plate crust. The low VP/VS anomalies reflect either a low-permeability overriding plate or a low fluid supply from slab dehydration in this region, and nevertheless indicate stronger materials above the plate interface. We speculate that it is the low-porosity overriding plate rather than the fluid-rich plate interface that increases interseismic slip deficit, as the high rigidity and density of the overriding plate favor the accumulation of elastic strain at the plate interface (Sallares and Ranero, 2019; Sun et al., 2020). More sophisticated laboratory and numerical experiments focusing on the overriding plate rigidity are needed to verify our speculation. Similar conditions have been observed beneath the southern North Island, New Zealand (Henrys et al., 2020) and southern Cascadia (Guo et al., 2021). The low-porosity overriding plate may result from limited fluid flux across the plate interface because the slab 31 dehydration appears to be weaker in the Kodiak segment compared to the Shumagin segment (Wei et al., 2021). 2.5. Conclusions We image the 3D high-resolution VP and VP/VS structure of the megathrust and the fore-arc overriding plate at the Alaska Peninsula. We also construct a new slip deficit model at seismogenic depths and find anti-correlation between VP/VS and slip deficit along the trench. Fluids indicated by high VP/VS anomalies play an important role in controlling the along-strike variations in interseismic megathrust slip behavior at the Alaska Peninsula. In the Shumagin segment to the west, abundant fluids at the plate interface and in the overriding plate lead to a low slip deficit, less frequent large megathrust earthquakes, and higher background seismic activity. In the Chignik, Chirikof, and Kodiak segment to the east, low fluid content in the overriding plate contributes to a higher slip deficit and more frequent occurrence of great megathrust earthquakes regardless of whether the plate interface is rich or starved in fluids. Since the overriding plate shows little along-strike changes in lithology, there is no evidence of permeability variations. Therefore, we suspect that the inferred porosity variations reflect the changes in the fluid supply from slab dehydration, and a sustained fluid flux across the plate interface is critical to enhancing interseismic slip. In light of this study, we suggest future geophysical and geological observations need to pay more attention to the fore-arc overriding plate, and numerical simulations of megathrust slip need to consider the complex along-strike heterogeneities of fluid content in the overriding plate. 32 FIGURES Figure 2.1. Topography and bathymetry map of the Alaska Peninsula and adjacent regions. The new along-strike segmentation of the megathrust is determined based on the seismic tomography and plate coupling models in this study. Note that we avoid using “Semidi segment”, a term commonly used in previous studies (Li et al., 2018a; Lynner, 2021; Shillington et al., 2015; Wei et al., 2021), because the Semidi Islands lie on the boundary between the Chignik and Chirikof segments. Dark black circles are seismic stations used in this study (Figure 2.2c). The magenta curves illustrate the slab surface (top of the plate interface) at 10, 20, 30, 40, and 50 km depths. The slab surface geometry is based on a model constrained by active-source seismic experiments (Kuehn, 2019) in addition to the Slab2 model (Hayes et al., 2018). Black dashed contours outline the rupture zones of historical megathrust earthquakes determined from 33 Figure 2.1 (cont’d) their aftershock distributions (Davies et al., 1981). Black solid contours show 1-m slip areas of previous megathrust earthquakes whose finite slip models have been recently determined by Freymueller et al. (2021), Xiao et al. (2021) and Elliott et al. (2022). Cyan-shaded ellipse shows the preferred area of a slow slip event in 2018 reported by He et al. (2023). Red-filled beachballs indicate the Global CMT focal mechanisms (Ekström et al., 2012) and hypocenters of the 2020 M7.8 Simeonof, 2021 M8.2 Chignik, and 2023 July M7.2 Sand Point megathrust earthquakes. 34 Figure 2.2. Event and station distributions used in this study. Distribution of global stations (a) and teleseismic events (b) are used in the first step of tomography for a large-scale velocity model. (c) Distribution of regional stations (blue triangles) and local events (red dots) used in the second step of tomography for high-resolution VP and VP/VS models. Other features are the same as in Figure 2.1. 35 Figure 2.3. Maps of VP and VP/VS models using different starting VP/VS values at constant depths of 10, 20, 30, and 40 km. (a) Final VP model. (b) VP/VS model inverted using the VP/VS values from the first step as the starting model. (c) Final VP/VS model inverted using a uniform value of 1.75 as the starting model. Bold black lines indicate the cross-sections in Figures 2.5 and 2.6. Other features are the same as in Figure 2.1. 36 Figure 2.4. Along-strike variations of VP and VP/VS structure in the plate interface and overriding plate. (a) VP structure at 20 km depth, representing the overriding plate. Regions with low resolutions are masked. The area southeast of the 20 km depth contour of the slab surface is semi-transparently masked because this area is within the subducted slab rather than the overriding plate. Dark yellow contours on the Pacific Plate outline seafloor magnetic anomalies (Meyer et al., 2017). Yellow-shaded areas show significant sediment input to the trench (von Huene et al., 2012). The thick brown contour up-dip of the 2021 M8.2 Chignik 37 Figure 2.4 (cont’d) rupture zone indicates the area of a slow slip event in 2018 (He et al., 2023). Bold black lines indicate the cross-sections in Figures 2.5 and 2.6. Other features are the same as Figures 2.1. (b) VP/VS structure on the slab surface. The slab surface geometry is based on a model constrained by active-source seismic experiments (Kuehn, 2019) in addition to the Slab2 model (Hayes et al., 2018). Other features are the same as in (a). (c) and (d) show the VP/VS structure in the overriding plate and at the plate interface, respectively. Other features are the same in (a) and (b). 38 Figure 2.5. VP structure along 6 cross-sections shown in Figure 2.4. Cross-sections A-A’ to E-E’ are normal to the trench, whereas F-F’ is along the 30-km depth contour of the slab surface and is perpendicular to all other cross-sections. In each panel, the grey curve shows the slab surface based on the Slab2 (Hayes et al., 2018) model and with the geometry changes constrained by previous active-source seismic studies (Kuehn, 2019), and the black parts indicate the rupture areas of previous megathrust earthquakes. The dark-green and red dashed curves show the bottom of the sediment layer and the overriding plate Moho, respectively, based on an independent surface-wave tomography study (Li et al., submitted). Vertical black dashed lines in F-F’ show the segment boundaries. Red bars above F-F’ indicate the rupture areas of previous megathrust earthquakes. Regions with low resolutions are masked. Topography/bathymetry is plotted on the top of each panel with vertical exaggeration. The relocated earthquake distribution (red dots) in each cross-section is right below the velocity model. 39 Figure 2.6. VP/VS structure along 6 cross-sections shown in Figure 2.4. Other features are the same as in Figure 2.5 except that the relocated events are plotted with the velocity structures together. Cross-section F-F’ is in the Chugach terrane from northeast to southwest despite the along-strike changes in VP/VS. The slip deficit from Figure 2.7 along each cross-section is shown as red curves atop the velocity structure. 40 Figure 2.7. New slip deficit model constrained by geodetic data from 1992 to 2017. Compared to the previous model by Drooff and Freymueller (2021), we add an additional segment boundary near Tugidak Island (See light blue arrow), where the new boundary is added between the Kodiak and Chirikof segments. The magenta ellipses and vectors indicate the data uncertainties and the misfit of GNSS stations using the new slip deficit model, respectively. 41 Figure 2.8. The comparison of the new slip deficit model with the VP/VS structure along the strike. The red, grey, blue, and dotted green curves indicate the average VP/VS on the slab surface, average VP/VS in the overriding plate, average slip deficit, and seismic activity, respectively. All the average values are derived by taking the average of the structure between 15 and 35 km depths along the slab surface or in the overriding plate. The number of earthquakes within a 30-km-wide along-dip cross-section is normalized by the largest number of earthquakes along the strike. 42 Figure 2.9. Schematic cartoons showing fluids at the plate interface and in the overriding plate enhancing interseismic slip on the megathrust. Semi-transparent blue shades mark the fluid-rich regions revealed by high VP/VS anomalies. Wavy arrows indicate fluid migration, and water drops represent the possible fluid migration from deep sources. Black contours outline the rupture areas of previous megathrust earthquakes shown in Figure 2.1, whereas the beachballs and stars represent the hypocenters of the 2020 M7.8 Simeonof and 2021 M8.2 Chignik 43 Figure 2.9 (cont’d) earthquakes. The dashed black contour up-dip of the 2021 M8.2 Chignik earthquake rupture zone shows the area of the 2018 slow slip event in the Chignik segment. The color-coded curves on each plate interface illustrate the slip deficit in each segment. Black lines on the surface represent the terrane boundaries between the Prince William terrane (PWT), Chugach terrane (CT), and Peninsula terrane (PT). Note that the inferred properties of the overriding plate crust do not correlate with geological terranes, suggesting that the overriding plate lithology is not the predominant factor in controlling the observed along-strike variations. 44 Figure 2.10. Checkerboard tests for VP and VP/VS models. The input models are checkerboard-shaped velocity anomalies in 3D with VP varying from -5% to 5% and VP/VS varying from -8% to 8%. (a) First checkerboard test. The size of the input velocity anomalies is 78, 77, and 15 km in the east-west, north-south, and vertical directions, respectively. (b) Second 45 Figure 2.10 (cont’d) checkerboard test. The anomaly size increases to 26, 22, and 10 km in the east-west, north-south, and vertical directions, respectively. Model resolvability contours of 0.54 (purple) in the first checkerboard test are chosen to mask out low-resolution regions. Grey dashed lines in the top panels indicate the vertical cross-sections shown in Figure 2.11. Other features are the same as in Figure 2.1. 46 Figure 2.11. Vertical cross-sections of the checkerboard tests of VP (a, b, e, f) and VP/VS (c, d, g, h) models in Figure 2.10. Red bars on the top of each panel indicate the rupture areas of previous megathrust earthquakes. Model resolvability contours of 0.54 in the first checkerboard test (purple) are used to mask out low-resolution regions. 47 Figure 2.12. Vertical cross-sections of the input and recovered model along A-A’ to F-F’ in Figures 2.5 and 2.6 from the first checkerboard test. The size of input velocity anomalies is 78, 77, and 15 km in the east-west, north-south, and vertical directions, respectively. (a) The 48 Figure 2.12 (cont’d) input velocity anomalies along cross-sections from A-A’ to E-E’. (b) The recovered velocity anomalies along cross-sections from A-A’ to E-E’. (c) The input velocity anomalies along cross- sections F-F’. (d) The recovered velocity anomalies along cross-sections F-F’. Other features are the same as in Figure 2.6. 49 Figure 2.13. A synthetic test of high VP/VS in the slab and high VP/VS sediment layer in the overriding crust. The plate interface is represented by a 10-km-thick high VP/VS layer immediately beneath the slab surface, and a 10-km thick high VP/VS layer on top of the overriding plate represents a hypothetically thick sediment layer. (a) Input model along cross- 50 Figure 2.13 (cont’d) sections A-A’ to E-E’ shown in Figure 2.4. (b) Tomographic inversion recovery model. Other features are the same as in Figures 2.5 and 2.6. 51 Figure 2.14. A restoration test for an input structure that mimics Figure 2.4d. (a) Input model with 10-km-thick high or low VP/VS anomalies. (b) Tomographic inversion recovery model. Other features are the same as in Figure 2.4. The short-wavelength structure near the edge of the structure boundary is related to the plotting artifact. 52 Figure 2.15. Vertical cross-sections of the input model along A-A’ to F-F’ in Figure 2.4 from the restoration test in Figure 2.14. The input 10-km thick VP/VS anomalies are put beneath the slab surface based on Slab2 with geometry changes constrained by the active-source imaging. Other features are the same as in Figure 2.6. 53 Figure 2.16. Vertical cross-sections of the recovered model along A-A’ to F-F’ in Figure 2.4 from the restoration test in Figure 2.14. Other features are the same as in Figure 2.6. 54 Figure 2.17. Comparison between observed seismic wave speed and theoretical values. Color-coded curves show VP and VP/VS along the slab surface of each cross-section in Figure 2.4. Dashed lines with numbers indicate theoretical VP and VP/VS of common rocks and mineral assemblies calculated following Abers and Hacker (2016) and the numbers on lines show the corresponding mineral assemblies. High VP/VS (> 1.8) anomalies indicate the presence of fluids. More details of theoretical calculations are described in Appendix Text 2.2.1. 55 Figure 2.18. Comparisons with previous studies in this region. (a) Cross-sections of our VP model along ALEUT Lines 3 and 5 from Shillington et al. (2022). A higher VP structure is imaged above the slab surface at ~35 km depth in ALEUT Line 5, which is consistent with the active source imaging results. (b) Cross-sections of our VP, VS and VP/VS models in the Kodiak segment, along the receiver function profile from Onyango et al. (2022). The solid grey and dashed magenta lines show the slab surface from the Slab2 model (Hayes et al., 2018) and the inferred Moho 8 km beneath the slab surface, respectively. A high VP/VS structure is imaged below the slab surface near the center of the profile, which is consistent with the receiver function image. 56 Figure 2.19. The analysis of seismic anisotropy. (a) and (b) indicate azimuth distributions of the seismic ray paths that sample the plate interface and the overriding plate, respectively. The rose diagrams from top to bottom correspond to the Shumagin, Chignik, Chirikof and Kodiak segments, respectively. The red line in each panel shows the fast direction of the azimuth 57 Figure 2.19 (cont’d) direction in each segment based on an independent study of surface waves (Liu et al., 2022b). (c) Incident angle distribution of the seismic ray paths that sample the plate interface. 90˚ indicates the vertical direction. The blue fan indicates the varying slab dip angle along the strike at 10–40 depths at the Alaska Peninsula. (d) The directions of the ray paths that sample beneath the bottom of the slab. We assume the slab is ~100 km thick. The red lines point at the directions of the ray paths traveling in the asthenosphere. Other features are the same as in Figure 2.1. 58 Figure 2.20. The VP/VS structure in different set up of the slab surface. (a) The VP/VS structure at the Slab2 (Hayes et al., 2018) surface. (b) The VP/VS structure at the slab surface constrained by the relocated seismicity distribution. The seismicity distribution is fitted to the 3rd-order polynomial with all the events right beneath the slab surface. (c), (d) and (e) show the VP/VS structure which takes the average of the velocity structure 10 km below and above the slab surface constrained by the relocated seismicity distribution, Slab2, and active-source imaging results, respectively. Other features are the same as in Figure 2.4. 59 Figure 2.21. The VP/VS structure in the plate interface with terrane features atop. The thick trench-parallel curves show the locations of the terrane boundary (Plafker et al., 1994). BRFS and CFS stand for Border Ranges fault system and Contact fault system, respectively. The cross- section F-F’ is located in the Chugach terrane. 60 Figure 2.22. VP/VS structure (a-e) along 5 cross-sections shown in (f). (f) shows the VP/VS structure at the slab surface. The locations of cross-sections B-B’ to E-E’ are slightly different from those in Figures 2.5-2.6, so that the new cross-sections cross the locked patches from Zhao et al. (2022) (purple contours in (f)). Along each cross-section, the black part of the slab surface (gray curve) represents the locked patches in (f) from Zhao et al. (2022). The plate coupling distribution from Zhao et al. (2022) along each cross-section is shown as red curves atop the velocity structure. Other features in (a)-(e) and (f) are the same as Figure 2.6 and Figure 2.4d, respectively. 61 REFERENCES Abers, G.A., van Keken, P.E., Hacker, B.R., 2017. The cold and relatively dry nature of mantle forearcs in subduction zones. Nat. Geosci. 10 (5), 333-337. https://doi.org/10.1038/ngeo2922 Argus, D.F., Gordon, R.G., Heflin, M.B., Ma, C., Eanes, R.J., Willis, P., Peltier, W.R., Owen, S.E., 2010. The angular velocities of the plates and the velocity of Earth's centre from space geodesy. Geophys. J. Int. 180 (3), 913-960. https://doi.org/10.1111/j.1365- 246X.2009.04463.x Arnulf, A.F., Bassett, D., Harding, A.J., Kodaira, S., Nakanishi, A., Moore, G., 2022. Upper- plate controls on subduction zone geometry, hydration and earthquake behaviour. Nat. Geosci. 15 (2), 143-148. https://doi.org/10.1038/s41561-021-00879-x Audet, P., Bostock, M.G., Christensen, N.I., Peacock, S.M., 2009. Seismic evidence for overpressured subducted oceanic crust and megathrust fault sealing. Nature 457 (7225), 76-78. https://doi.org/10.1038/nature07650 Audet, P., Burgmann, R., 2014. Possible control of subduction zone slow-earthquake periodicity by silica enrichment. Nature 510 (7505), 389-392. https://doi.org/10.1038/nature13391 Bassett, D., Sandwell, D.T., Fialko, Y., Watts, A.B., 2016. Upper-plate controls on co-seismic slip in the 2011 magnitude 9.0 Tohoku-oki earthquake. Nature 531 (7592), 92-96. https://doi.org/10.1038/nature16945 Bassett, D., Sutherland, R., Henrys, S., 2014. Slow wavespeeds and fluid overpressure in a region of shallow geodetic locking and slow slip, Hikurangi subduction margin, New Zealand. Earth Planet. Sci. Lett. 389, 1-13. https://doi.org/10.1016/j.epsl.2013.12.021 Chesley, C., Naif, S., Key, K., Bassett, D., 2021. Fluid-rich subducting topography generates anomalous forearc porosity. Nature 595 (7866), 255-260. https://doi.org/10.1038/s41586- 021-03619-8 Cordell, D., Naif, S., Evans, R., Key, K., Constable, S., Shillington, D., Bécel, A., 2023. Forearc seismogenesis in a weakly coupled subduction zone influenced by slab mantle fluids. Nat. Geosci. https://doi.org/10.1038/s41561-023-01260-w Cross, R.S., Freymueller, J.T., 2008. Evidence for and implications of a Bering plate based on geodetic measurements from the Aleutians and western Alaska. J. Geophys. Res. 113 (B7), B07405. https://doi.org/10.1029/2007jb005136 Davies, J., Sykes, L., House, L., Jacob, K., 1981. Shumagin seismic gap, Alaska Peninsula: History of great earthquakes, tectonic setting, and evidence for high seismic potential. J. Geophys. Res. 86 (B5), 3821-3855. https://doi.org/10.1029/JB086iB05p03821 62 Delph, J.R., Levander, A., Niu, F., 2018. Fluid controls on the heterogeneous seismic characteristics of the Cascadia margin. Geophys. Res. Lett. 45 (20), 11,021-011,029. https://doi.org/10.1029/2018gl079518 Drooff, C., Freymueller, J.T., 2021. New Constraints on Slip Deficit on the Aleutian Megathrust and Inflation at Mt. Veniaminof, Alaska From Repeat GPS Measurements. Geophys. Res. Lett. 48 (4), e2020GL091787. https://doi.org/10.1029/2020gl091787 Egbert, G.D., Yang, B., Bedrosian, P.A., Key, K., Livelybrooks, D.W., Schultz, A., Kelbert, A., Parris, B., 2022. Fluid transport and storage in the Cascadia forearc influenced by overriding plate lithology. Nat. Geosci. 15 (8), 677-682. https://doi.org/10.1038/s41561- 022-00981-8 Ekström, G., Nettles, M., Dziewoński, A.M., 2012. The global CMT project 2004–2010: Centroid-moment tensors for 13,017 earthquakes. Phys. Earth Planet. Inter. 200-201, 1- 9. https://doi.org/10.1016/j.pepi.2012.04.002 Elliott, J., Freymueller, J.T., 2020. A block model of present‐day kinematics of Alaska and western Canada. J. Geophys. Res. Solid Earth 125 (7), e2019JB018378. https://doi.org/10.1029/2019jb018378 Elliott, J.L., Grapenthin, R., Parameswaran, R.M., Xiao, Z., Freymueller, J.T., Fusso, L., 2022. Cascading rupture of a megathrust. Sci. Adv. 8 (18), eabm4131. https://doi.org/10.1126/sciadv.abm4131 Fletcher, H.J., Freymueller, J.T., 2003. New constraints on the motion of the Fairweather fault, Alaska, from GPS observations. Geophys. Res. Lett. 30 (3), 1139. https://doi.org/10.1029/2002gl016476 Freymueller, J.T., Suleimani, E.N., Nicolsky, D.J., 2021. Constraints on the slip distribution of the 1938 Mw 8.3 Alaska Peninsula earthquake from tsunami modeling. Geophys. Res. Lett. 48 (9), e2021GL092812. https://doi.org/10.1029/2021gl092812 Gou, T., Xia, S., Huang, Z., Zhao, D., 2022. Structural heterogeneity of the Alaska‐Aleutian forearc: Implications for interplate coupling and seismogenic behaviors. J. Geophys. Res. Solid Earth 127 (11), e2022JB024621. https://doi.org/10.1029/2022jb024621 Guo, H., McGuire, J.J., Zhang, H., 2021. Correlation of porosity variations and rheological transitions on the southern Cascadia megathrust. Nat. Geosci. 14 (5), 341-348. https://doi.org/10.1038/s41561-021-00740-1 Hayes, G.P., Moore, G.L.P., Daniel E., Hearne, M.F., Hanna., Furtney, M.S., Gregory M., 2018. Slab2, a comprehensive subduction zone geometry model. Science 362 (6410), 58-61. https://doi.org/10.1126/science.aat4723 He, B., Wei, X., Wei, M., Shen, Y., Alvarez, M., Schwartz, S.Y., 2023. A shallow slow slip 63 event in 2018 in the Semidi segment of the Alaska subduction zone detected by machine learning. Earth Planet. Sci. Lett. 612, 118154. https://doi.org/10.1016/j.epsl.2023.118154 Heise, W., Caldwell, T.G., Bannister, S., Bertrand, E.A., Ogawa, Y., Bennie, S.L., Ichihara, H., 2017. Mapping subduction interface coupling using magnetotellurics: Hikurangi margin, New Zealand. Geophys. Res. Lett. 44 (18), 9261-9266. https://doi.org/10.1002/2017gl074641 Henrys, S., Eberhart‐Phillips, D., Bassett, D., Sutherland, R., Okaya, D., Savage, M., Evanzia, D., Stern, T., Sato, H., Mochizuki, K., Iwasaki, T., Kurashimo, E., Seward, A., Wech, A., 2020. Upper plate heterogeneity along the southern Hikurangi margin, New Zealand. Geophys. Res. Lett. 47 (4), e2019GL085511. https://doi.org/10.1029/2019gl085511 Holt, A.F., Condit, C.B., 2021. Slab temperature evolution over the lifetime of a subduction zone. Geochem. Geophys. Geosyst. 22 (6), e2020GC009476. https://doi.org/10.1029/2020gc009476 Holtkamp, S.G., 2017. Slow Slip and Mixed Seismic/Aseismic Moment Release Episodes on the Alaskan-Aleutian Megathrust, EarthScope National Meeting, Anchorage, Alaska. Horowitz, W.L., Steffy, D.A., Hoose, P.J., Turner, R.F., 1989. Geologic report for the Shumagin Planning Area, Western Gulf of Alaska, OCS Report, MMS 89-0097, U.S. Department of the Interior, Minerals Management Service, Alaska OCS Region. Ikari, M.J., Saffer, D.M., Marone, C., 2009. Frictional and hydrologic properties of clay-rich fault gouge. J. Geophys. Res. 114 (B5). https://doi.org/10.1029/2008jb006089 Jarrard, R.D., 1986. Relations among subduction parameters. Rev. Geophys. 24 (2), 217-284. https://doi.org/10.1029/RG024i002p00217 Jiang, Y., González, P.J., Bürgmann, R., 2022. Subduction earthquakes controlled by incoming plate geometry: The 2020 M>7.5 Shumagin, Alaska, earthquake doublet. Earth Planet. Sci. Lett. 584. https://doi.org/10.1016/j.epsl.2022.117447 Katayama, I., Kubo, T., Sakuma, H., Kawai, K., 2015. Can clay minerals account for the behavior of non-asperity on the subducting plate interface? Prog. in Earth and Planet. Sci 2 (30). https://doi.org/10.1186/s40645-015-0063-4 Kodaira, S., Iidaka, T., Kato, A., Park, J.-O., Iwasaki, T., Kaneda, Y., 2004. High pore fluid pressure may cause silent slip in the Nankai Trough. Science 304 (5678), 1295-1298. https://doi.org/10.1126/science.1096535 Kuehn, H., 2019. Along-trench segmentation and down-dip limit of the seismogenic zone at the eastern Alaska-Aleutian subduction zone. Dalhousie University, p. 334. Lay, T., Kanamori, H., Ammon, C.J., Koper, K.D., Hutko, A.R., Ye, L., Yue, H., Rushing, T.M., 64 2012. Depth-varying rupture properties of subduction zone megathrust faults. J. Geophys. Res. 117 (B04311), B04311. https://doi.org/10.1029/2011jb009133 Li, J., Shillington, D.J., Saffer, D.M., Bécel, A., Nedimović, M.R., Kuehn, H., Webb, S.C., Keranen, K.M., Abers, G.A., 2018a. Connections between subducted sediment, pore- fluid pressure, and earthquake behavior along the Alaska megathrust. Geology 46 (4), 299-302. https://doi.org/10.1130/g39557.1 Li, S., Freymueller, J., McCaffrey, R., 2016. Slow slip events and time‐dependent variations in locking beneath Lower Cook Inlet of the Alaska‐Aleutian subduction zone. J. Geophys. Res. Solid Earth 121 (2), 1060-1079. https://doi.org/10.1002/2015jb012491 Li, Z., Wiens, D., Shen, W., Shillington, D., submitted. Along-strike variations of Alaska subduction zone structure and hydration determined from amphibious seismic data. Lindsey, E.O., Mallick, R., Hubbard, J.A., Bradley, K.E., Almeida, R.V., Moore, J.D.P., Bürgmann, R., Hill, E.M., 2021. Slip rate deficit and earthquake potential on shallow megathrusts. Nat. Geosci. 14 (5), 321-326. https://doi.org/10.1038/s41561-021-00736-x Liu, C., Lay, T., Xiong, X., 2022a. The 29 July 2021 Mw 8.2 Chignik, Alaska Peninsula earthquake rupture inferred from seismic and geodetic observations: Re‐rupture of the western 2/3 of the 1938 Rupture Zone. Geophys. Res. Lett. 49 (4). https://doi.org/10.1029/2021gl096004 Liu, C., Zhang, S., Sheehan, A.F., Ritzwoller, M.H., 2022b. Surface wave isotropic and azimuthally anisotropic dispersion across Alaska and the Alaska‐Aleutian subduction zone. J. Geophys. Res. Solid Earth 127 (11), e2022JB024885. https://doi.org/10.1029/2022jb024885 Lynner, C., 2021. Anisotropy-revealed change in hydration along the Alaska subduction zone. Geology 49 (9), 1122-1125. https://doi.org/10.1130/g48860.1 McCaffrey, R., 2002. Crustal block rotations and plate coupling, in: Stein, S., Freymueller, J.T. (Eds.), AGU Geodyn. Series, pp. 101-122. Meyer, B., Chulliat, A., Saltus, R., 2017. Derivation and Error Analysis of the Earth Magnetic Anomaly Grid at 2 arc min Resolution Version 3 (EMAG2v3). Geochem. Geophys. Geosyst. 18 (12), 4522-4537. https://doi.org/10.1002/2017gc007280 Molina-Ormazabal, D., Ampuero, J.-P., Tassara, A., 2023. Diverse slip behavior of velocity- weakening fault barriers. Nat. Geosci. 16 (12), 1200-1207. https://doi.org/10.1038/s41561-023-01312-1 Moreno, M., Haberland, C., Oncken, O., Rietbrock, A., Angiboust, S., Heidbach, O., 2014. Locking of the Chile subduction zone controlled by fluid pressure before the 2010 earthquake. Nat. Geosci. 7 (4), 292-296. https://doi.org/10.1038/ngeo2102 65 Nishikawa, T., Matsuzawa, T., Ohta, K., Uchida, N., Nishimura, T., Ide, S., 2019. The slow earthquake spectrum in the Japan Trench illuminated by the S-net seafloor observatories. Science 365 (6455), 808-813. https://doi.org/10.1126/science.aax5618 Okada, Y., 1992. Internal deformation due to shear and tensile faults in a half-space. Bull. Seismol. Soc. Am. 82 (2), 1018-1040. https://doi.org/10.1785/BSSA0820021018 Onyango, E.A., Worthington, L.L., Schmandt, B., Abers, G.A., 2022. Subduction zone interface structure within the southern MW 9.2 1964 great Alaska earthquake asperity: Constraints from receiver functions across a spatially dense node array. Geophys. Res. Lett. 49 (15), e2022GL098334. https://doi.org/10.1029/2022gl098334 Pesicek, J.D., Zhang, H., Thurber, C.H., 2014. Multiscale Seismic Tomography and Earthquake Relocation Incorporating Differential Time Data: Application to the Maule Subduction Zone, Chile. Bull. Seismol. Soc. Am. 104 (2), 1037-1044. https://doi.org/10.1785/0120130121 Plafker, G., Moore, J.C., Winkler, G.R., Plafker, G., Berg, H.C., 1994. Geology of the southern Alaska margin, The Geology of Alaska. Geological Society of America, p. 0. Reece, R.S., Gulick, S.P.S., Christeson, G.L., Horton, B.K., van Avendonk, H., Barth, G., 2013. The role of farfield tectonic stress in oceanic intraplate deformation, Gulf of Alaska. J. Geophys. Res. Solid Earth 118 (5), 1862-1872. https://doi.org/10.1002/jgrb.50177 Ruppert, N.A., Barcheck, G., Abers, G.A., 2022. Enhanced regional earthquake catalog with Alaska Amphibious Community Seismic Experiment data. Seismol. Res. Lett. 94 (1), 522-530. https://doi.org/10.1785/0220220226 Saffer, D.M., Tobin, H.J., 2011. Hydrogeology and mechanics of subduction zone forearcs: Fluid flow and pore pressure. Annu. Rev. Earth Planet. Sci. 39 (1), 157-186. https://doi.org/10.1146/annurev-earth-040610-133408 Sallares, V., Ranero, C.R., 2019. Upper-plate rigidity determines depth-varying rupture behaviour of megathrust earthquakes. Nature 576 (7785), 96-101. https://doi.org/10.1038/s41586-019-1784-0 Shillington, D.J., Bécel, A., Nedimović, M.R., 2022. Upper plate structure and megathrust properties in the Shumagin Gap near the July 2020 M7.8 Simeonof event. Geophys. Res. Lett. 49 (2), e2021GL096974. https://doi.org/10.1029/2021gl096974 Shillington, D.J., Bécel, A., Nedimović, M.R., Kuehn, H., Webb, S.C., Abers, G.A., Keranen, K.M., Li, J., Delescluse, M., Mattei-Salicrup, G.A., 2015. Link between plate fabric, hydration and subduction zone seismicity in Alaska. Nat. Geosci. 8 (12), 961-964. https://doi.org/10.1038/ngeo2586 66 Suito, H., Freymueller, J.T., 2009. A viscoelastic and afterslip postseismic deformation model for the 1964 Alaska earthquake. J. Geophys. Res. Solid Earth 114 (B11), B11404. https://doi.org/10.1029/2008jb005954 Sun, T., Saffer, D., Ellis, S., 2020. Mechanical and hydrological effects of seamount subduction on megathrust stress and slip. Nat. Geosci. 13 (3), 249-255. https://doi.org/10.1038/s41561-020-0542-0 Takei, Y., 2002. Effect of pore geometry on VP/VS: From equilibrium geometry to crack. J. Geophys. Res. 107 (B2), ECV-6. https://doi.org/10.1029/2001jb000522 Tobin, H.J., Saffer, D.M., 2009. Elevated fluid pressure and extreme mechanical weakness of a plate boundary thrust, Nankai Trough subduction zone. Geology 37 (8), 679-682. https://doi.org/10.1130/g25752a.1 von Huene, R., Miller, J.J., Krabbenhoeft, A., 2021. The Alaska convergent margin backstop splay fault zone, a potential large tsunami generator between the frontal prism and continental framework. Geochem. Geophys. Geosyst. 22 (1), e2019GC008901. https://doi.org/10.1029/2019gc008901 von Huene, R., Miller, J.J., Weinrebe, W., 2012. Subducting plate geology in three great earthquake ruptures of the western Alaska margin, Kodiak to Unimak. Geosphere 8 (3), 628-644. https://doi.org/10.1130/ges00715.1 Vrolijk, P., 1987. Tectonically driven fluid flow in the Kodiak accretionary complex, Alaska. Geology 15 (5), 466-469. https://doi.org/10.1130/0091- 7613(1987)15<466:TDFFIT>2.0.CO;2 Wang, K., Wells, R., Mazzotti, S., Hyndman, R.D., Sagiya, T., 2003. A revised dislocation model of interseismic deformation of the Cascadia subduction zone. J. Geophys. Res. Solid Earth 108 (B1). https://doi.org/10.1029/2001jb001227 Wei, S.S., Ruprecht, P., Gable, S.L., Huggins, E.G., Ruppert, N., Gao, L., Zhang, H., 2021. Along-strike variations in intermediate-depth seismicity and arc magmatism along the Alaska Peninsula. Earth Planet. Sci. Lett. 563, 116878. https://doi.org/10.1016/j.epsl.2021.116878 Witter, R.C., Briggs, R.W., Engelhart, S.E., Gelfenbaum, G., Koehler, R.D., Barnhart, W.D., 2014. Little late Holocene strain accumulation and release on the Aleutian megathrust below the Shumagin Islands, Alaska. Geophys. Res. Lett. 41 (7), 2359-2367. https://doi.org/10.1002/2014gl059393 Xiao, Z., Freymueller, J.T., Grapenthin, R., Elliott, J.L., Drooff, C., Fusso, L., 2021. The deep Shumagin gap filled: Kinematic rupture model and slip budget analysis of the 2020 Mw 7.8 Simeonof earthquake constrained by GNSS, global seismic waveforms, and floating InSAR. Earth Planet. Sci. Lett. 576, 117241. https://doi.org/10.1016/j.epsl.2021.117241 67 Ye, L., Lay, T., Kanamori, H., Yamazaki, Y., Cheung, K.F., 2021. The 22 July 2020 M 7.8 Shumagin seismic gap earthquake: Partial rupture of a weakly coupled megathrust. Earth Planet. Sci. Lett. 562, 116879. https://doi.org/10.1016/j.epsl.2021.116879 Zelt, C.A., 1998. Lateral velocity resolution from three-dimensional seismic refraction data. Geophys. J. Int. 135 (3), 1101-1112. https://doi.org/10.1046/j.1365-246X.1998.00695.x Zhang, H., Thurber, C.H., 2003. Double-difference tomography: The method and its application to the Hayward fault, California. Bull. Seismol. Soc. Am. 93 (5), 1875-1889. https://doi.org/10.1785/0120020190 Zhao, B., Burgmann, R., Wang, D., Zhang, J., Yu, J., Li, Q., 2022. Aseismic slip and recent ruptures of persistent asperities along the Alaska-Aleutian subduction zone. Nat Commun 13 (1), 3098. https://doi.org/10.1038/s41467-022-30883-7 Zhao, D., Katayama, Y., Toyokuni, G., 2022. The Moho, slab and tomography of the East Japan forearc derived from seafloor S-net data. Tectonophysics 837, 229452. https://doi.org/10.1016/j.tecto.2022.229452 68 APPENDIX Text 2.1. Seismic tomography details 2.1.1. Seismic travel-time tomography In this study, we image the high-resolution 3D seismic velocity structure of the Alaska Peninsula and the adjacent regions using body-wave travel-time data recorded by the Alaska Amphibious Community Seismic Experiment (AACSE) and the adjacent stations of the EarthScope Transportable Array and the Alaska Earthquake Center (Ruppert et al., 2022) (Figure 2.2). We conduct a body-wave travel-time tomography imaging on multiple scales. First, we develop a large-scale velocity model using regional and teleseismic data. Second, we use this model as a starting model to image the VP and VP/VS structure down to 200 km depth with only regional seismic data. With this strategy, our final results are largely independent of the starting model. In the first step, we assemble the P- and S-wave arrival times from May 2018 through August 2019 recorded by onshore and offshore seismic stations of the AACSE, EarthScope USArray, and the Alaska regional network (Figure 2.2c). These arrivals are manually picked by the Alaska Earthquake Center (Ruppert et al., 2022). We correct for the reported clock drifts of OBSs (Barcheck et al., 2020) and gather 302,976 P-wave and 129,073 S-wave arrival times from 7,242 local events in the Alaska Peninsula (Figure 2.2c). We also collect teleseismic arrival times consisting of two parts. The first part includes the direct P, S, and depth phases (pP and sS) from Alaska earthquakes recorded by global seismic stations from 2009 to 2019 and reported by the International Seismological Center (ISC, https://doi.org/10.31905/D808B830) (Figure 2.2a). The second part includes the P and S phases from the teleseismic events with magnitudes larger than 6 and epicentral distance from 30° to 90° recorded by the regional stations from May 2018 69 through August 2019 (Figure 2.2b). We filter the waveforms of these events between 0.1~1 Hz for P waves and 0.05~0.15 Hz for S waves, and utilize an adaptive stacking method (Rawlinson and Kennett, 2004) to obtain P and S-wave arrival times. Only the arrival times with small picking errors (< 0.3 s) are kept for further analysis. Thus, the teleseismic dataset includes 272,074 direct P, 9,857 direct S, 2,267 pP, and 30 sS arrival times. We use an improved version of teletomoDD (Pesicek et al., 2014), a teleseismic double- difference tomography package, to invert the regional and teleseismic arrival times for 3-D VP and VS models down to 700 km depth. This tomography algorithm requires two sets of grids of nodes for regional and global scales, respectively. The regional grid extends from sea level to 700 km depth, with horizontal and vertical spacing of 0.3º and 20 km, respectively. The global grid is sparser, with a 5˚ lateral spacing and a 90-km vertical spacing, extending to 2889 km depth. We construct the starting velocity model by combining the top 50 km part of a recent regional VS model constrained by ambient noise and teleseismic surface waves (Li et al., submitted) and a 3D global model TX2019slab below 50 km depth (Lu et al., 2019). We utilize a bi-weight function to down-weight the arrival times of observations with large residuals and apply a first-order smoothing constraint to slowness perturbations to stabilize the tomographic inversion. We perform four iterations of joint inversion for velocity and relocation and two iterations of inversion for only relocation. The damping and smoothing parameters are chosen via a trade-off analysis of model perturbations and travel time residuals. The root-mean-square (RMS) value of differential time residuals decreases from the initial value of 1.476 s to the final value of 0.485 s. Because the S-wave dataset is much smaller than the P-wave dataset, the output VP model has a much higher resolution compared to the VS model. In the second step, we invert for higher-resolution VP and VP/VS models from 0 to 200 km 70 depths using only regional seismic data, using the output models from the first step as the starting model, and tomoDD (Guo et al., 2018; Zhang and Thurber, 2003), a regional double- difference tomography package. We densify the inversion grid to 0.2º spacing in longitude and 0.1º spacing in latitude. In the vertical direction, the grid spacing is 5 km in the top 50 km, 10 km from 50 to 80 km, and 20 km between 80 and 200 km depth. The tomographic inversion includes 9 iterations of joint inversion for velocity and relocation as well as 6 iterations of inversion for only relocation. The VP and VP/VS models are always inverted simultaneously to avoid mapping errors in one model into the other. We choose the regularization parameters according to the trade-off between model perturbations and time residuals. Because the VP and VS models from the first step have different resolutions, using them to construct a starting VP/VS model will produce some unrealistic values and potentially bias inversion in the second step. Therefore, we set the starting VP/VS model to have a uniform value of 1.75. As shown in Figure 2.3, the output VP/VS models given two different starting models are similar. Therefore, we prefer the inverted VP and VP/VS models with a starting VP/VS value of 1.75 (Figures 2.3a, 2.3c, 2.4-2.6). 2.1.2. Seismic tomography resolution tests We perform 4 synthetic tests to evaluate our seismic-velocity model resolution. The first two tests are designed to evaluate the resolution limit of our tomography model using different sizes of input velocity anomalies (Figures 2.10-2.12). In the first checkerboard test, the size of the input velocity anomalies is 78, 77, and 15 km in the east-west, north-south, and vertical directions, respectively (Figure 2.10a). In the second checkerboard test, the anomaly size decreases to 26, 22, and 10 km in the east-west, north-south, and vertical directions, respectively (Figure 2.10b). We construct the synthetic model by adding positive and negative 5% velocity perturbations to the starting VP and VS models. Then we compute synthetic arrival times using 71 the same distribution of stations and events as the real data and add Gaussian-distributed random noise. Finally, we invert the synthetic data using the same inversion parameters as in the real data inversion. We also conduct a third resolution test for the ability to recover a high VP/VS plate interface and a high VP/VS upper crust in the overriding plate by our inversion system (Figure 2.13). The input model consists of a 10-km-thick high VP/VS (~1.89) layer below the slab surface and another 10-km-thick high VP/VS layer in the overriding plate. The ambient model has a constant VP/VS at 1.75. We compute synthetic data and conduct inversions following the same process as the checkerboard test. The fourth resolution test is designed to validate the along-strike variations of VP/VS at the plate interface (Figures 2.14-2.16). We create a 3D input VP/VS model simplified from the real- data inversion results. In the Shumagin segment (cross-section A-A’ in Figures 2.14-2.16), a 10- km-thick high VP/VS layer at the plate interface extends from near the trench to 50 km depth. In the Chignik and Kodiak segments, a 10-km-thick low VP/VS layer at the plate interface at 20–30 km depths is sandwiched by two 10-km thick high VP/VS layers at the plate interface in both up- dip and down-dip directions (cross-section C-C’ and E-E’ in Figures 2.14-2.16, respectively). In the Chirikof segment (cross-section D-D’ in Figures 2.14-2.16), a 10-km-thick low VP/VS layer at the plate interface extends from 10 km to 40 km and in the up-dip direction of a 10-km-thick high layer at the plate interface. The ambient model has a constant VP/VS at 1.75. We compute synthetic data and conduct inversions following the same process as the checkerboard test. Text 2.2. Geological implications of VP/VS anomalies and effects of seismic anisotropy 2.2.1. Geological implications of VP and VP/VS anomalies In this study, we focus on the VP/VS images while also using the VP model to help with 72 interpretation. We estimate the VP and VP/VS of potential compositions on the Alaska slab surface using a MATLAB toolbox developed by Abers and Hacker (2016). The thermal structure of the slab surface is simulated by Syracuse et al. (2010). We investigate most common minerals and rocks in subduction zones, such as zoisite facies, prehnite-pumpellyite facies, lawsonite blueschist facies, greenschist facies, lherzolite, harzburgite, and antigorite. The VP and VP/VS profiles of various minerals and rocks at different depths are shown as the solid or the dashed lines in Figure 2.17. The VP and VP/VS on the slab surface along all cross-sections A-A’ to E-E’ from our tomographic model are plotted for comparison. The extremely high VP/VS (>1.83) along cross-sections A-A’ and B-B’ above 30 km depth cannot be explained by any minerals and rocks. Thus, we suggest that the extremely high VP/VS anomalies may indicate the existence of fluids on the slab surface. The relatively higher VP/VS anomalies along cross-sections A-A’, B- B’, and E-E’ than along cross-sections C-C’ and D-D’ indicate a higher fluid content at the slab surface in the Shumagin and Kodiak segments than in the Chignik and Chirikof segments. In addition, highly damaged rocks and fluid-filled cracks also could result in high VP/VS anomalies (Miller et al., 2021). Therefore, we interpret high VP/VS anomalies as a proxy for fluids in highly deformed rocks and cracks without explicitly estimating the rock porosity. 2.2.2. Effects of seismic anisotropy on our interpretations Rock experiments have reported that seismic anisotropy can increase the measured VP/VS of rocks when the ray path is perpendicular to crack fabric and structural foliation and thus the fast direction of anisotropy (Miller et al., 2021; Wang et al., 2012). Previous seismic studies have observed complex patterns of azimuthal anisotropy in this region (Liu et al., 2022b; Lynner, 2021). In the Shumagin segment, azimuthal distributions of the seismic rays sampling the plate interface and overriding plate are well balanced between the trench-parallel and trench-normal 73 directions (Figures 2.19a, 2.19b). In the Chignik and Chirikof segments, more seismic rays are trench-parallel than trench-normal. Given that the fast direction of azimuthal anisotropy is nearly trench-parallel in the overriding plate and trench-normal in the slab (Liu et al., 2022b), the measured VP/VS may artificially decrease in the overriding plate and increase in the slab. In other words, the true structure may show moderately low VP/VS anomalies in the overriding plate and even lower VP/VS anomalies at the slab interface compared to those shown in Figures 2.4-2.6. In the Kodiak segment, the seismic rays are predominantly trench-normal. With a similar azimuthal anisotropy as in the Chignik and Chirikof segments, the true structure in the Kodiak segment may have even lower VP/VS anomalies in the overriding plate and even higher VP/VS anomalies at the slab interface compared to those shown in Figures 2.4-2.6. Thus, we believe azimuthal anisotropy does not change our interpretations of the along-strike variations in VP/VS. Separately, Li et al. (2018b) suggest that tilted transversely isotropy, possibly caused by structural foliation, widely exists in subducted slabs. However, the seismic rays that sample the plate interface are well balanced between the slab-normal and slab-parallel directions (Figure 2.19c). Therefore, the tilted transversely isotropy in the slab should have minor effects on our measured VP/VS. Previous studies have also reported strong sub-slab anisotropy in the asthenosphere (Song and Kawakatsu, 2012, 2013). However, only 0.047% of the seismic rays used in our regional high-resolution models (Figures 2.4–2.6) travel beyond the bottom of the subducted slab, and most of these ray paths are located beneath the Shumagin segment and to its west (Figure 2.19d). Therefore, the sub-slab anisotropy in the asthenosphere should not affect our VP and VP/VS models along the strike and dip directions. In summary, we conclude that anisotropy has only minor effects on our VP and VP/VS models (Figure 2.19). 74 CHAPTER 3: SLAB MORPHOLOGY, DEHYDRATION, AND SUB-ARC MELTING BENEATH THE ALASKA PENINSULA REVEALED BY BODY-WAVE TOMOGRAPHY An edited version of this chapter is under review at the Journal of Geophysical Research: Solid Earth: Fan Wang, S. Shawn Wei, Natalia A. Ruppert, and Haijiang Zhang. Slab morphology, dehydration, and sub-arc melting beneath the Alaska Peninsula revealed by body-wave tomography 75 Abstract The Alaska Peninsula has a long history of plate subduction and diverse arc volcanism along the strike. The Alaska slab morphology below 200 km depth remains debated due to limited seismic data and thus low tomography resolution in this region. In addition to the diverse arc volcanism on the surface, the sub-arc melting processes that feed these volcanoes are unclear. Here we utilize the newly available regional and teleseismic data to build 3-D high-resolution VP and VS models to 660 km depth. We find that the high-velocity Pacific Plate subducts to the bottom of the mantle transition zone (MTZ) with complex deformation and gaps. In the southwest and northeast, we observe vertical gaps between the high-velocity slab above 200 km and within the MTZ. In the center to northeast, the Pacific Plate is a continuous slab subducted to the MTZ with an abrupt increase in the dip angle at 200 km depth. We also invert for 3-D VP and VP/VS models to 200 km depth with higher resolutions and interpret strong along-strike changes in slab dehydration and sub-arc melting, based on low VP and high VP/VS anomalies. Slab dehydration and sub-arc melting are most extensive below the Pavlof and Shumagin segments in the southwest, weakening below the Chignik and Chirikof segments in the northeast, and strengthening again beneath the Kodiak segment further to the northeast. We propose that the variations of slab hydration at the outer rise significantly influence slab dehydration at greater depths and further control sub-arc melting beneath the Alaska Peninsula. 76 3.1. Introduction Slab morphology in the deep mantle reflects the mantle rheology and dynamics of mantle processes. Subducted slabs can display large variations in its morphology when they interact with mantle flow while sinking into the deep mantle. A slab may stagnate at the base of the MTZ or penetrate into the lower mantle, according to seismic tomography and geodynamic modeling (Billen, 2008; Fukao & Obayashi, 2013; Goes et al., 2017; Lay, 1994; Wei et al., 2020). Furthermore, slabs also display more complex deformation and dynamics, such as buckling, bending, and tearing (Obayashi et al., 2017; Ye et al., 2016; Zhang et al., 2019). Moreover, fossil slabs from ancient subduction events also exist in the mantle at various depths (Cai & Wiens, 2016; Jia et al., 2020; Chen & Michael, 2001). Due to the limited resolution of global and regional tomography models, slab morphology in the deep mantle remains poorly constrained in many subduction zones, leading to large uncertainties in studies of plate reconstruction and mantle dynamics. The subduction of a slab enables efficient material exchange between the Earth’s surface and deep mantle, particularly for volatiles such as water. An oceanic plate gets hydrated in the form of hydrous minerals through hydrothermal reactions (Faccenda, 2014). The most effective hydration of a plate happens at the outer rise where plate-bending faults are created (Faccenda, 2014; Shillington et al., 2015). During subduction, the subducted slab releases water into the mantle wedge and the crust of the overriding plate during dehydration reactions of hydrous minerals when the pressure-temperature conditions are met, possibly triggering intermediate- depth earthquakes (Hacker et al., 2003; Yamasaki & Seno, 2003). The addition of water triggers flux melting of the mantle wedge and slab crust materials beneath arc volcanoes, returning water back to the Earth's surface (Grove et al., 2009; van Keken et al., 2011). Within the overriding 77 plate, hydrous melts migrate upward from crystal-rich mushes in the lower crust to melt- dominated magma chambers in the upper crust (Cashman et al., 2017). Although the reactive flow in the crust can dramatically change the melt chemistry, the mantle source places the upper limit of volatiles entering the transcrustal magmatic system. The Alaska Peninsula section of the Alaska-Aleutian subduction zone shows a long history of plate subduction (Plafker et al., 1994; Rioux et al., 2010; Singer et al., 2007). Yet the deep slab morphology remains unclear due to the complicated subduction history and limited seismic imaging in the past. Plate reconstruction models show multi-stage subductions of the Kula and Pacific Plates since the Cretaceous (Madsen et al., 2006; Müller et al., 2016). The Kula Plate subducted from ~110-80 Myr ago to ~40 Myr ago. One or more mid-ocean ridge systems subducted along the margin between Alaska and the Pacific Northwest during this period. The Pacific Plate subduction started ~48 Myr ago and had a major change of subduction direction relative to North America ~42 Myr ago. Since ~40 Myr ago, the configuration of Pacific-North American plate motion has been stable and similar to the present (Madsen et al., 2006). Therefore, at least a 2,000-km-long Pacific Plate has been subducted beneath the North American Plate given the convergence rate of 50 mm/yr (Pavlis et al., in press). Since intermediate-depth earthquakes only extend to ~200 km depth in this subduction zone (Wei et al., 2021), the locations of the subducted Kula and Pacific Plates below 200 km, and especially the slab morphology in the mantle transition zone (MTZ) beneath the Alaska Peninsula, are unclear (Boyce et al., 2023; Burdick et al., 2017; Gou et al., 2019; Jiang et al., 2018; Martin‐Short et al., 2018; Qi et al., 2007). To better describe the along-strike changes of the Alaska Peninsula section of the Alaska-Aleutian subduction zone, we divide this region into five along-strike segments from southwest to northeast following Wang et al. 2024: the Pavlof, Shumagin, 78 Chignik, Chirikof, and Kodiak segments (Figure 3.1). The morphology of the Alaska slab in the deep mantle remains debated according to previous seismic tomography studies using different datasets and methods (Boyce et al., 2023; Burdick et al., 2017; Gou et al., 2019; Jiang et al., 2018; Madsen et al., 2006; Martin‐Short et al., 2018; Müller et al., 2016; Qi et al., 2007; Shephard et al., 2017). Among P-wave tomography studies, the earlier work by Qi et al. (2007) using travel-time data recorded at inland stations imaged all of Alaska, revealing a high-VP slab subducted to ~400 km depth with a normal dip angle (~45°) beneath the Kodiak segment and a separate flattened high-VP structure in the MTZ to the northwest with a gap between them. They suggested that the flattened high-VP structure in the MTZ and above 400 km represents the extinct Kula Plate and the subducted Pacific Plate down to 300-400 km, respectively, and the gap in between presents the ancient Kula-Pacific spreading center. In contrast, recent travel-time tomography models of North America using data recorded by USArray and global stations (Burdick et al., 2017 and Boyce et al, 2023) suggested a continuous Alaska slab extending to ~600 km depth with a steep dip angle beneath the Alaska Peninsula. Burdick et al. (2017) found the Alaska slab subducting steeply without much folding in the MTZ, whereas Boyce et al. (2023) observed a less continuous high- VP structure in the MTZ. Moreover, Gou et al. (2019) also reported an Alaska slab subducting to ~660 km and flattening in the MTZ using P-wave travel-time data from Alaska inland stations. Most published S-wave tomography studies about slab deep morphology used ambient noise, teleseismic Rayleigh waves, and teleseismic S-wave travel-time data, sometimes combined with teleseismic receiver functions, recorded at Alaska inland stations: Martin‐Short et al. (2018) reported an Alaska slab subducting to ~450-550 km with a steep dip angle (~70˚) beneath the Alaska Peninsula, whereas depth extend of the slab gets shallower to its northeast. In 79 comparison, Jiang et al. (2018) observed that the Alaska slab subducts to ~400-500 km depth with a normal dip angle (~45°) beneath the Kodiak segment. In summary, all the regional seismic tomography studies suggest that the Alaska slab subducts northwestward and does not penetrate into the lower mantle beneath the Alaska Peninsula but the slab shows different morphologies in the MTZ. In contrast, a global tomography model vote map indicates a diffusive south-dipping structure in the lower mantle beneath the Kodiak-Prince William Sound region (Shephard et al., 2017). It remains unknown if the northwest-dipping high-VP structure in the MTZ is connected to the lower mantle velocity structure. At depths shallower than 200 km, the Alaska slab displays large variations in physical and chemical properties along the strike across different segments of the Alaska Peninsula. Previous studies have investigated the changes in sediment thickness, rupture processes of large megathrust earthquakes, slab slip behavior, background seismicity, as well as fluids distribution in the slab and overriding plate (Cordell et al., 2023; Drooff & Freymueller, 2021; Li et al., 2018; Li et al., 2024; Liu et al., 2022; Lynner, 2021; Shillington et al., 2015; Wang et al., 2024; Wei et al., 2021; Xiao et al., 2021). The Alaska slab is presumably more hydrated at the outer rise and thus undergoes more intensive dehydration after subduction in the Shumagin segment to the southwest than in the Chignik and Chirikof segments to the northeast (Li et al., 2024; Lynner, 2021; Shillington et al., 2015; Wang et al., 2024; Wei et al., 2021). Shillington et al. (2015) suggest that the along-strike variation in the slab hydration state, indicated by low-velocity anomalies in the incoming plate crust and uppermost mantle, coincides with the changes in the pre-existing fabric orientation on the incoming Pacific seafloor. In the Pavlof and Shumagin segments, the seafloor fabric has a low angle (10–20˚) with the trench, and thus more plate- bending faults are created to allow for extensive hydration of the incoming plate. In comparison, 80 the seafloor fabric has a high angle (50–70˚) with the trench in the Chignik and Chirikof segments. Consequently, fewer bending faults are developed or reactivated at the outer rise, resulting in low-degree hydration of the incoming plate. Sediment input into the trench also varies along the strike due to the subducted Zodiac fan and Surveyor fan (von Huene et al., 2012; Li et al., 2018). Thicker sediment subducts in the Chignik segment than in the Shumagin segment because of the Zodiac fan, while sediment subduction in the Kodiak segment is controlled by the thick Surveyor fan. Most of these studies focus on slab hydration and dehydration above ~50–80 km depths, but slab dehydration at greater depths and its influences on sub-arc melting are poorly known. Wei et al. (2021) investigated intermediate-depth earthquakes at 50–200 km depths and geochemical signatures of trace elements of arc lavas. With indirect evidence, they hypothesized that along-strike changes in slab dehydration at depths greater than 50 km are affected by both pre-existing fabric in the incoming plate and the subducted sediments. A comprehensive analysis of slab dehydration below 50 km and its relationship to material recycling in the sub-arc is needed given the new seismic data and high- resolution tomography images. Arc volcanoes in the Alaska Peninsula show diverse eruption styles and strong along-strike changes in geochemical signatures, which may be related to the melting processes in the sub-arc mantle and continental crust (Larsen, 2016; Wei et al., 2021). Okmok and Aniakchak Volcanoes in the west form large-volume, melt-rich calderas and erupt effusive to highly explosive caldera- forming events, whereas Augustine Volcano in the east is a smaller-volume, crystal-choked andesite stratovolcano (Global Volcanism Program, 2024). It is under debate whether the along- strike volcanic changes are mainly caused by mantle or crustal processes (Buurman et al., 2014). Previous studies have reported a higher amount of fluids subducted in the southwest segments 81 and a higher volume of sediments subducted in the northeast segments (Wei et al., 2021), implying that the water input into the mantle wedge controls the sub-arc melting processes. However, direct evidence of changes in the sub-arc mantle melting source and its effects on arc volcanism is still missing. In this study, we utilize both regional and teleseismic P- and S-wave travel-time data to image large-scale 3-D VP and VS structures from the Earth’s surface to 660 km depth beneath the Alaska Peninsula, thanks to the newly acquired offshore and onshore data. Then, we develop new VP and VP/VS models to 200 km depth with an even higher resolution. The new 3-D velocity models provide important constraints on slab morphology, slab dehydration, as well as sub-arc melting beneath the Alaska Peninsula. Here, we focus on the slab interior, mantle wedge, and sub-arc crust, as the seismogenic zone above 50 km depth and the fore-arc region has been thoroughly investigated by a separate study (Wang et al., 2024). 3.2. Data and Methods 3.2.1 Data In this study, we used both teleseismic and regional arrival-time data to image the slab and mantle wedge structure at multiple scales. The teleseismic arrival-time dataset consisted of two parts. The first part included 824,382 P-, 10,911 pP-, 55,497 S-, and 827 sS-wave arrival-times from Alaska earthquakes recorded by global seismic stations from 2000 to 2020 and reported by the International Seismological Centre (Figure 3.2b). The second part included P and S waves from teleseismic events with magnitudes larger than 6 and epicentral distances from 30° to 90° recorded by regional stations from May 2018 through August 2019 (Figure 3.2a). We filtered the waveforms of these events at 0.1–1 Hz for P waves and 0.05–0.15 Hz for S waves. Then we utilized an adaptive stacking method (Rawlinson & Kennett, 2004) to automatically pick P- and 82 S-wave arrival-times. Arrival-times with large picking errors (> 0.3 s) were discarded. This second part consisted of 12,220 P- and 11,177 S-wave teleseismic arrival-times that sampled the deep structure beneath the Alaska Peninsula. For the regional arrival-time dataset, we collected P- and S-wave arrival-times of local earthquakes from May 2018 through August 2019 recorded by onshore and offshore seismic stations of the Alaska Amphibious Community Seismic Experiment (AACSE), EarthScope USArray, and the Alaska regional network (Figure 3.2c). These arrivals were manually picked by the Alaska Earthquake Center (Ruppert et al., 2022). We corrected the reported OBS clock drifts (Barcheck et al., 2020) and gathered 302,976 P- and 129,073 S-wave arrival-times from 7,242 local events in the Alaska Peninsula. 3.2.2 Teleseismic double-difference VP and VS tomography In the first stage, we use an improved version of teletomoDD (Pesicek et al., 2010), a teleseismic double-difference tomography package, to invert the regional and teleseismic arrival- times for 3-D VP and VS models down to 660 km depth. We perform 3D ray tracing for both teleseismic and regional data using the pseudo-bending algorithm within a nested regional and global 3D velocity model. The regional velocity model extends from the western end of the Alaska Peninsula to the eastern end of the Kenai Peninsula, with a horizontal spacing of 0.3º. Vertically, the model extends to 660 km depth with a spacing of 20 km for the top 140 km and a larger spacing of 30~50 km at greater depths. We construct the starting regional velocity model by combining the top 50 km of a recent regional model constrained by ambient noise and receiver functions (Li et al., 2024) with the 3D global TX2019slab model (Lu et al., 2019) below 50 km depth. The starting global model is built based on the 3D global TX2019slab (Lu et al., 2019) model with much sparser nodes with a horizontal spacing of 5˚ and a vertical spacing of 83 ~90 km below 100 km depth. The regional model is embedded within the global model to allow for ray tracing of teleseismic waves outside of the regional model. We perform 4 iterations of joint inversions for velocity and relocation and 2 iterations of inversions only for relocation. We down-weight the arrival-times with large residuals using a bi- weight function and apply a first-order smoothing constraint to the slowness of model perturbations to stabilize the inversion. The damping and smoothing parameters are chosen with a trade-off analysis of model perturbations and travel-time residuals. The root means square (RMS) value of differential time residuals decreases from the initial value of 1.081 s to the final value of 0.446 s for P waves, and from 1.462 s to 0.511 s for the S waves. Because the S-wave dataset is much smaller than the P-wave dataset, the output VP model has a much higher resolution compared to the VS model. 3.2.3 Regional-scale double-difference VP and VP/VS tomography for the top 200 km structure In the second stage, we invert the regional dataset for higher-resolution VP and VP/VS models from 0 to 200 km depth using a regional-scale double-difference tomography package tomoDD (Guo et al., 2018; Zhang & Thurber, 2003). The starting 3D VP and VP/VS velocity models are adapted from the VP and VS models from the first stage. We densify the inversion grid to 0.2º spacing in longitude and 0.1º spacing in latitude. In the vertical direction, the grid spacing is 5 km in the top 50 km, 10 km from 50 to 80 km depths, and 20 km between 80 and 200 km depths. Since the VP and VS models from the first stage have different resolutions, using them to create a starting VP/VS model may contain unrealistic pre-existing structures and potentially bias the VP/VS inversion. Therefore, we also try another series of inversions with the starting VP/VS value being uniformly 1.75 to avoid potential biases caused by the starting model. We jointly invert for velocity and relocation in 9 iterations and then invert only for relocation in 6 iterations. We 84 choose the regularization parameters based on the trade-off between model perturbations and travel-time residuals. 3.3. Results 3.3.1 Large-scale VP and Vs structures imaged with teletomoDD Our teletomoDD results from the first stage show VP structure down to 660 km depth and VS structure in a more confined area down to 400 km depth (Figures 3.3 and 3.4). We image a dipping structure with high VP and VS from the trench to at least 200 km depth beneath the Alaska Peninsula, interpreted to be the subducted Alaska slab. The top of the slab above 200 km depth coincides with the Slab2 model (Hayes et al., 2018). Due to the lower ray coverage below 200 km depth, the amplitude of high-VP below 200 km (~1%) is consistently lower than that above 200 km (~3%) and the VS model does not show any high-velocity anomalies below 200 km. The slab morphology revealed by high-velocity anomalies shows obvious along-strike changes. At depths shallower than 200 km, the slab dip angle increases from ~10° in the Kodiak segment (represented by cross-sections E-E’ and F-F’ in Figure 3.4) southwestward to ~30° in the Pavlof segment (represented by cross-section A-A’). In the Chirikof segment (represented by cross-section D-D’), the slab extends to at least 660 km depth and its dip angle increases abruptly from ~45° to nearly ~75° at ~200-300 km depths in the VP model (Figure 3.4a). In the Kodiak segment to the northeast (represented by cross-sections E-E’ and F-F’), the slab seems to be less continuous below 200 km depth, with a few gaps separating high-velocity anomalies down to the bottom of the MTZ. In the Chignik, Shumagin, and Pavlof segments (represented by cross- sections C-C’, B-B’, and A-A’, respectively), the high-velocity slab appears to end at 200 km depth, and there is another high-VP anomaly in the lower MTZ at 500-700 km depth, separated 85 by a ~300-km-wide vertical gap from the slab structure above. From the map view (Figure 3.3), this slab gap is at least ~200 km wide at 300-500 km depth in the Chignik, Shumagin, and Pavlof segments. We observe low-VP and low-VS anomalies in the overriding plate’s crust and mantle wedge above 300 km depth (Figure 3.4). The lowest VP and VS structures are commonly located beneath arc volcanoes in the Alaska Peninsula (Figures 3.3a and 3.4), consistent with previous studies in this region (Eberhart-Phillips et al., 2006; Gou et al., 2019; Martin‐Short et al., 2018; Qi et al., 2007; Tian & Zhao, 2012; Yang & Gao, 2020). We also observe widespread low-VP anomalies in the mantle wedge at 60-300 km depths, extending northwestward to at least ~150 km away from the slab (Figures 3.3, 3.4). Low VP anomalies are also imaged immediately beneath the bottom of the high-VP slab, particularly in the Pavlof, Shumagin, and Kodiak segments. 3.3.2 Synthetic tests for the slab morphology We perform 5 synthetic tests of VP tomography, including 2 checkerboard tests and 3 restoration tests, to evaluate the robustness of the slab morphology from the teletomoDD inversion. The recovered checkerboard models indicate that the slab morphology is well-imaged in our study (Figure 3.5). In the first checkerboard test (Figures 3.5a, 3.5b), the input velocity anomalies have the size of 145, 99, and ~100 km in the east-west, north-south, and vertical directions, respectively. In the second checkerboard test (Figures 3.5c, 3.5d), we decrease the size of input velocity anomalies to ~68, ~67, and 60 km in the east-west, north-south, and vertical directions, respectively. We build the synthetic models by adding positive or negative 5% velocity perturbations to the starting VP and VS models. We then compute synthetic P-wave travel-times by assuming the same distributions of events and stations as in the real data. Finally, we invert the synthetic data with the same inversion strategy and regularization parameters as in 86 the real data inversion. The checkerboard results show that at depths shallower than 200 km sampled by regional P waves, we can achieve a resolution of ~68 km and ~60 km in the horizontal and vertical directions, respectively. At depths below 200 km, we can resolve velocity anomalies with a size of ~150 km horizontally and ~100 km vertically. Overall, we have a better resolution beneath the Shumagin and Kodiak segments due to more local events there (Figure 3.5). By visually comparing the input and recovered models, we choose the number of sampling rays larger than 20 as the threshold for the model resolvability and use this contour to mask out the low-resolution regions in the final VP and VS models. We also conduct 3 restoration tests to validate the slab morphology beneath 200 km (Figure 3.6). In all of these restoration tests, we calculate the synthetic arrival-times based on different configurations of slab morphology and perform the inversions following the same procedure as in the checkerboard tests. The input model in our first resolution test consists of a 100-km-thick high VP slab extending to 200 km depth following the Slab2 model (Hayes et al., 2018), with the ak135-f model (Kennett, 1995) as the background model (Figure 3.6a). The output model shows a well-recovered high-VP slab subducting to 200 km in all cross-sections (Figure 3.6b). The high- VP slab above 200 km would not be mistakenly mapped to the MTZ depths. Therefore, the high- VP anomalies imaged in the MTZ from the real data inversion are reliable structures and are not artifacts caused by smoothing in the tomographic inversions (Figures 3.3, 3.4). The second resolution test is designed to evaluate if we can resolve a continuous slab to the MTZ. Similar to the first restoration test, we extend the hypothetical slab to 660 km depth in the input model (Figure 3.6c). The output model indicates that the continuous high-VP slab structure can be well- recovered to the bottom of the MTZ (Figure 3.6d). If the Alaska slab is a continuous high-VP structure extending to 660 km depth, no gap in the high-velocity slab would be imaged along any 87 cross-sections. Therefore, the gap imaged from the real data inversion indicates that the Alaska slab is not continuous at 200-500 km depths beneath the Pavlof, Shumagin, Chignik, and Kodiak segments. It is worth noting that the recovered high-VP slab has a smaller amplitude and becomes more diffusive below 200 km than the slab structure above, which is also found in the real data inversion. It results from a lower resolution due to lower data coverage below 200 km. In the third resolution test, we impose a slab gap at 200–300 km depths between two high-velocity structures (Figure 3.6e). The output model shows that the gap can be slightly healed with vertical smoothing (Figures 3.6e, 3.6f). Therefore, the observed slab gaps at 200–500 km depths in the real data inversion results could be wider than they appear in the image. In summary, the checkerboard results suggest that we can image the slab morphology and mantle wedge structure at a spatial resolution of ~150 km horizontally at depths below 200 km and achieve a higher resolution of ~68 km above 200 km depth. In addition, we are able to image gaps wider than 100 km in the high-velocity slab in the vertical direction and detect high VP structures in the MTZ from our large-scale teletomoDD VP model. 3.3.3 High-resolution VP and VP/VS structures imaged with regional tomoDD The high-resolution VP and VP/VS models extend from 0 to 200 km depth. We prefer the inverted VP/VS model with a starting VP/VS of 1.75 to avoid any potential bias introduced by the initial model, because the VP/VS model shows similar structures regardless of the choice of the initial model, using either a constant value of 1.75 or the VP/VS values from the first stage of teletomoDD inversions (Figures 3.3 and 3.4). Since our tomography method does not include a discrete slab boundary in the model parameterization, we independently construct a slab surface (i.e., the top of the plate interface) model for further interpretations. This slab surface model starts from the Slab2 model (Hayes et al., 2018), and then we update the slab geometry in the 88 Shumagin, Chignik, and Chirikof segments based on more accurate active-source seismic imaging (Kuehn, 2019). Consequently, the slab surface is ~7 km deeper than that in the Slab2 model in the Chignik and Chirikof segments, and transitions to the Slab2 model within a 50-km- wide buffer zone outside of the active-source imaging region. We observe low VP (< 7.7 km/s) and high VP/VS (> 1.80) anomalies beneath the slab surface with strong along-strike variations in our final models (Figures 3.8, 3.9). In the seismogenic zone at depths shallower than 50 km, VP immediately beneath the slab surface is ~7 km/s, more than 0.8 km/s lower than the averaged VP value in the slab at the same depths. VP/VS shows high values (> 1.84) near the trench in most segments. At depths greater than 50 km, the low VP and high VP/VS anomalies beneath the slab surface extend to different depths and change their amplitude and thickness in different segments. In the Pavlof segment to the southwest (cross- section a-a’ in Figures 3.8 and 3.9), we observe a >10-km-thick low VP (< 7.8 km/s) layer beneath the slab surface extending to ~80 km depth. High VP/VS (~1.84) anomalies are confined within ~20 km beneath the slab surface and extend to ~200 km depth, accompanied by abundant intermediate-depth earthquakes that form a double seismic zone. The upper-plane earthquakes coincide with high VP/VS anomalies, whereas the lower-plane events are located in a region of low VP/VS at ~60-170 km depths. In the Shumagin segment (cross-section b-b’ in Figures 3.8 and 3.9), a 10-km-thick low VP (< 7.9 km/s) layer underlays the slab surface and extends to ~70 km depth. We observe high VP/VS anomalies (~1.82-1.84) within ~10-20 km beneath the slab surface from the trench to ~80 km depth. At depths below ~120 km, VP/VS is neutral-to-high (~1.80- 1.82) beneath the slab surface where seismic activity decreases in this region. In the Chignik segment (cross-sections c-c’ in Figures 3.8 and 3.9), we observe a thick (~20 km) low VP (< 7.7 km/s) layer beneath the slab surface extending to > 90 km depth. VP/VS is generally low in the 89 slab, whereas sporadic high VP/VS (~1.82) anomalies are confined in 10-km-thick narrow zones beneath the slab surface and show lower amplitude than that in the segments to the west. We do not observe a double seismic zone in this region as suggested by Wei et al. (2021) due to the limited events recorded during the AACSE deployment. In the Chirikof segment further to the east (cross-section d-d’ in Figures 3.8 and 3.9), a low VP (< 7.9 km/s) layer beneath the slab surface is <10 km thick and extends to ~80-100 km depths. The slab is overall characterized by low VP/VS, whereas high VP/VS (~1.82-1.84) anomalies are confined in ~10-km-thick regions beneath the slab surface at depths of ~40 km and ~130-150 km. Intraplate earthquakes are abundant at depths below 40 km and extend to at least 180 km depth. In the Kodiak segments to the northern end (cross-sections e-e’ and f-f’ in Figures 3.8 and 3.9), a 10-km-thick low VP (< 7.8 km/s) layer extends to ~60 km depth. We observe high VP/VS anomalies within 20 km beneath the slab surface and extending to ~120 km depth, accompanied by abundant intermediate-depth earthquakes. In summary, we observe a low VP (< 7.9 km/s) layer extending to at least 60 km depth in all the segments. High VP/VS anomalies in the slab are most extensive in the Pavlof and Shumagin segments, and only sporadically found in the Chignik and Chirikof segments. The slab in the Kodiak segment shows high VP/VS anomalies in its top 10-20 km beneath the slab surface. We also image low VP and high VP/VS anomalies in the crust of the overriding plate and the mantle wedge at 50-100 km depths beneath all the segments (Figures 3.8, 3.9). The continental Moho is 30–50 km deep, constrained by an independent study of surface wave and receiver functions (magenta curves in Figures 3.8 and 3.9) (Li et al., 2024). In the overriding crust, we observed lower VP (6.0 km/s, less than 0.3 km/s lower than the mean value of the overriding crust at corresponding depths) and high VP/VS (> 1.82) structures near the trench and beneath the arc volcanoes. Below ~40 km depth in the mantle wedge, we image low VP (< 7.9 km/s) and high 90 VP/VS (> 1.82) anomalies with strong-along strike variations, which show a spatial correlation with the low VP and high VP/VS structures in the slab. In the Pavlof segment to the southwest, VP shows low values (< 7.8 km/s) and VP/VS is high (> 1.84) in the lower crust and at 60-100 km depth in the mantle wedge. The high VP/VS structures extend over broad areas in the mantle wedge (>60 km horizontally and vertically) and into the slab interior. In the Shumagin segment, an extremely low VP (~7.6 km/s) structure extends >100 km horizontally and vertically in the mantle wedge, where VP/VS is slightly higher than average. The overriding crust displays weak low VP and high VP/VS features beneath the Veniaminof volcano. In the Chignik segment, the extremely low VP (~7.6 km/s) structure extends over a large area in the mantle wedge, and high VP/VS (~1.82) anomalies are distributed at 80-120 km depths above the slab. Lower VP and high VP/VS are found in the overriding crust beneath the Aniakchak volcano, coinciding with a cluster of seismic activity. In the Chirikof segment, low VP (~7.8 km/s) structures are found at 40-100 km depth, whereas high VP/VS (1.82-1.84) anomalies are distributed at 80-140 km depths downdip of the low VP anomalies in the mantle wedge. The high VP/VS structure at 150-180 km depths extends from the mantle wedge into the slab interior. We also observe low VP and high VP/VS in the overriding crust beneath the Ukinrek Maars. In the Kodiak segment to the northeast, we find widely distributed low VP (~7.8 km/s) anomalies in the mantle wedge where VP/VS is on average (1.70–1.80). Beneath the Trident volcano, low VP anomalies are confined in the upper crust whereas high VP/VS appears only in the lower crust. In summary, we observe strong and extensive low VP anomalies in the mantle wedge in all segments except for the Pavlof segment, where high VP/VS anomalies are imaged. Low VP anomalies are found in the upper crust beneath active volcanoes, whereas abnormally high VP/VS features are confined in the lower crust beneath the Pavlof and Trident volcanoes. 91 3.3.4 Synthetic tests for the VP and VP/VS inversion We perform 4 synthetic tests of VP and VP/VS tomography, including 2 checkerboard tests and 2 restoration tests, to evaluate our model resolution for the second stage of regional scale tomography inversions (Figures 3.10, 3.11). The tests suggest that we can well image the overriding crust, slab, and mantle wedge from 0 to at least 180 km depth, particularly near earthquake sources (Figure 3.10). In the first checkerboard test (Figures 3.10a-f), the size of the input velocity anomalies is 26, 22 in the east-west and north-south directions. In the vertical direction, the grid spacing is 10 km above 50 km, 20 km from 50 km to 80 km, and 40 km below 80 km. In the second checkerboard test (Figures 3.10g-l), the anomaly size increases to 78 and 77 in the east-west and north-south directions, respectively. In the vertical direction, the grid spacing is 20 km above 40 km, 40 km from 40 km to 80 km, and 80 km below 80 km. We construct the synthetic model by adding alternating positive or negative 5% velocity perturbations to the starting VP and VP/VS models. We then compute synthetic P- and S-wave travel-times assuming the same event and station distributions as in the real data. We invert VP and VP/VS structures using the same inversion parameters as in the real data inversion. Finally, we calculate the semblance value between the recovered and input model at each node, serving as a proxy for model resolution. By visually comparing the input and recovered models, we choose 0.54 as the threshold of semblance value and use this contour to mask out low-resolution regions in the final VP and VP/VS models (Zelt, 1998). As shown in Figure 3.10, we achieve a resolution of 25 and 20 km (10 km above 80 km depth) in the horizontal and vertical directions, respectively, in the resolved region. Therefore, we focus on discussing intermediate-scale (≥30 km horizontally and ≥20 km vertically) structures in our models because some small-scale (<20 km horizontally) anomalies may reflect model resolution changes along the strike and dip 92 directions due to the ununiform distribution of stations and earthquakes. We also conduct 2 restoration tests to evaluate our ability to image a high VP/VS plate interface (Figure 3.11). We calculate synthetic arrival-times for each input model with different thicknesses of the high VP/VS layer, and then perform the tomography inversions following the same procedure as for the real data. In the first restoration test, the input model consists of a 10- km-thick high VP/VS (~1.89) layer below the slab surface subducted to 200 km depth. The ambient model has a constant VP/VS of 1.75. This test shows that the 10-km-thick high VP/VS layer at the plate interface can be well recovered at 15–80 km depths throughout the region (Figures 3.11a–b). In the second restoration test, the input model consists of a 15-km-thick high VP/VS (~1.89) layer below the slab surface and another 10-km-thick high VP/VS layer in the overriding plate. The resolution test confirms that we are able to image the high VP/VS layer subducted beyond 120 km. The 15-km-thick high VP/VS layer is well recovered to 80-120 km depth and not affected by the 10-km-thick high VP/VS layer in the overriding crust (Figures 3.11c, 3.11d). It indicates the high VP/VS anomaly in the overriding plate, representing a possible fore- arc sedimentary basin, will not be mistakenly mapped to greater depths. Therefore, we are confident with the imaging results near the plate interface. In summary, we can image the slab interior and mantle wedge structure with a resolution of ≥30 km horizontally and ≥20 km vertically in all the parts of the resolved regions. Therefore, we are confident with all moderate-scale anomalies described in Section 3.3.3. 3.4. Discussion 3.4.1 Slab morphology beneath the Alaska Peninsula We image a high VP structure extending to as deep as 660 km depth and a high VS structure extending to at least 200 km depth, both interpreted as the subducted Alaska slab. The 93 discrepancy between the VP and VS models results from the different teleseismic ray coverage. Since the VP model has a higher resolution and a larger resolved region compared to the VS model, we rely on the P-wave tomography results to interpret the slab morphology below 200 km depth. It is worth noting that the top of the slab coincides with the Slab2 model (Hayes et al., 2018), even though our tomography inversions are independent of any models of the slab morphology. Our results indicate that the Alaska slab has obvious morphological changes along the strike (Figure 3.12). Above 200 km depth, the slab dip angle decreases northeastward from ~30° to 10°, consistent with previous studies of seismicity (Hayes et al., 2018; Wei et al., 2021). At 200–300 km depths, the Alaska slab changes its dip angle abruptly from ~45° to ~75° in regions where the slab extends beyond 200 km depth. This change in the slab dip angle has been reported in previous P- and S-wave tomography studies (Boyce et al., 2023; Burdick et al., 2017; Martin‐ Short et al., 2018). In all segments, our VP model shows a slab structure in the MTZ, although it is not always connected to the Alaska slab above 200 km depth (Figures 3.4, 3.12). We interpret the lack of high-VP anomalies at 200-500 km depths as a slab gap. In the Pavlof, Shumagin, and Chignik segments, the ~300-km-wide slab gap separates the Alaska slab above 200 km from the high-VP structure in the lower MTZ. In the Chirikof segment, the Alaska slab continuously extends to at least 660 km depth. In the Kodiak segment, the slab gap is imaged at 200–300 km depths (Figure 3.4a). The slab morphology in the MTZ beneath the Alaska Peninsula remains debated in previous studies due to the lack of data and high-resolution tomography in this region. In the Kodiak segment to the northeast, regional tomography studies suggest that the Alaska slab subducts to either the top of the MTZ (Jiang et al., 2018; Qi et al., 2007), or ~450–550 km depths in the 94 MTZ (Martin‐Short et al., 2018), or the bottom of the MTZ (Boyce et al., 2023; Burdick et al., 2017; Gou et al., 2019). Southwest of the Kodiak segment, Boyce et al. (2023) suggest the Alaska slab subducts to the bottom of MTZ while other regional tomography models do not have a sufficient resolution in this region. In addition to seismic tomography, receiver function analyses report contradictory measurements of the MTZ thickness, resulting in different interpretations of the slab morphology in the MTZ. Hao et al. (2023) observe a thinned MTZ and suggest that the Alaska slab has not reached the MTZ in the Alaska Peninsula, whereas Dahm et al. (2017) report a thicker-than-normal MTZ and suggest the existence of slab segments in the MTZ. Van Stiphout et al. (2019) find a complex topography of the 410 discontinuity beneath the Kodiak segment, and interpret it as a complex discontinuity or an artifact caused by 3D velocity corrections. These discrepancies mainly result from the lack of local seismic stations before the AACSE deployment, and most previous studies exclude the Alaska Peninsula from their model parameterizations or put it near the edge of the models. Our results provide new constraints on the variations of deep slab morphology in the MTZ beneath the Alaska Peninsula. The restoration tests confirm that we are able to resolve a slab subducted to the bottom of the MTZ or a slab gap at 300–500 km depths (Figure 3.6). Our VP model indicates that the Alaska slab reaches the MTZ, and we do not observe a flattened slab feature or a southward extension dipping structure in the MTZ. Although the Alaska slab appears to be continuously extending to the MTZ in the Chirikof segment, we observe slab gaps wider than 100 km vertically in the Pavlof, Shumagin, Chignik, and Kodiak segments (Figure 3.12). These slab gaps have not been resolved by previous studies, as they have limited resolution at the MTZ depths due to the lack of offshore data. The strong along-strike variations in the deep slab morphology reflect the complicated 95 subduction history along the Alaska-Aleutian Trench, including the subduction of the Kula Plate and Pacific Plate beneath the Alaska Peninsula. The location of the subducted Kula Plate remains debated (Fuston & Wu, 2020; Qi et al., 2007). Since the Pacific Plate has been subducting at a constant speed of ~50 mm/yr beneath the North American Plate since 40 Ma age (Pavlis et al., in press), at least a 2,000-km-long Pacific Plate should have subducted into the Alaska-Aleutian Trench. Our tomography results show a maximum slab length of ~1,000 km along the subduction direction assuming no slab gaps, folding, or buckling. Besides, we observe the Alaska slab increases its thickness by less than twofold and is much less diffuse than previous models in the MTZ (Boyce et al., 2023; Burdick et al., 2017; Gou et al., 2019). Therefore, the high-VP structure in the MTZ should represent the subducted Pacific Plate rather than the Kula Plate, and the prominent slab gap above the MTZ may imply significant slab tearing in the upper mantle. Although the slab is continuous in the Chirikof segment, it is significantly thinner than above and below, possibly implying an intermediate stage where the Alaska slab is being torn at 200–500 km depths. Qi et al. (2007) reported a high VP structure in the MTZ northwest of the subducted Pacific Plate with a gap in between in the Kodiak segment. They interpret this high VP structure as the subducted Kula Plate. However, we do not observe such high-velocity structures north of the Alaska slab or other high-velocity anomalies in the MTZ beneath the Alaska Peninsula (Figure 3.3). The subducted Kula Plate may have sunk into the lower mantle or moved outside of our study area. 3.4.2 Along-strike changes in slab dehydration Our VP and VP/VS models, along with the relocated earthquake distribution, reveal details of the Alaska slab dehydration. VP is affected by pressure, temperature, and lithology, whereas VP/VS is more sensitive to the existence of fluids or hydrous minerals (Takei, 2002). In a 96 subduction zone setting, low VP and high VP/VS anomalies are usually associated with the presence of fluids or hydrous minerals such as serpentine (Nakajima et al., 2001; Peacock et al., 2011; Zhang et al., 2004). In contrast, high VP and low VP/VS anomalies are usually interpreted as a cold, dry subducted plate. The combination of low VP and low VP/VS is caused by a greater reduction of VP than VS. Zhang et al. (2004) suggest that the low VP and low VP/VS anomalies in the lower plane of the Honshu double seismic zone are related to the forsterite-enstatite-H2O formation from serpentine dehydration where fluids are present in spheroids. The combination of high VP and high VP/VS can be caused by fluid-filled cracks in cold rocks (Nugraha et al., 2019). Seismic anisotropy potentially has strong effects on VP/VS (Miller et al., 2021; Wang et al., 2012). Based on seismic ray path distribution, we conclude that anisotropy has minor effects on VP/V at the plate interface and overriding plate above 50 km depth (see detailed discussions in Wang et al. (2024)). However, anisotropy at greater depths may artificially increase or decrease VP/VS where ray paths are dominated by one direction. Besides, dehydration reactions of hydrous minerals in a slab release fluids at 50–300 km depths and help facilitate earthquakes by reducing effective normal stress where brittle failure is impossible in dry conditions (Hacker et al., 2003; Yamasaki & Seno, 2003). Therefore, VP, VP/VS, and the abundance of intermediate-depth earthquakes can provide important constraints on the dehydration processes in a slab, although cautions need to be taken regarding the effects of anisotropy on VP/VS. We observe low VP and high VP/VS anomalies in the slab interior with strong variations along the strike, indicating different degrees of slab dehydration (Figures 3.8, 3.9). Within the slab crust, the ultra-low-VP layer extends to 60–70 km depths, indicating the onset of crustal eclogitization (Yamasaki & Seno, 2003). In this layer, the extremely low VP (more than 0.8 km/s slower than the mean value) and high VP/VS (>1.82) anomalies indicate not only hydrous crustal 97 minerals but also aqueous fluids released from sediment compaction (Figures 3.8, 3.9). The slab mantle dehydration shows variations in depth and intensity, indicated by high VP/VS anomalies at various depths and with different amplitudes and thicknesses. The high VP/VS anomalies within the slab are unlikely caused by seismic anisotropy because the seismic rays sampling this region are dominantly parallel to the slab, in which case mineral foliation should artificially decrease rather than increase VP/VS (Miller et al., 2021). In the Pavlof segment to the west (cross-section a-a’ in Figures 3.8 and 3.9), slab dehydration extends to at least ~200 km depth where high VP/VS anomalies are imaged ~20 km beneath the slab surface in the slab-normal direction. They indicate that intensive dehydration happens to ~13 km beneath the slab Moho, as a result of extensive hydration of the slab crust and mantle through abundant normal faults in the outer rise (Shillington et al., 2015). The abundant intermediate-depth earthquakes and the well-developed double seismic zone are also consistent with deep slab dehydration towards ~200 km depth (Wei et al., 2021). It is worth noting that the lower plane of the Pavlof double seismic zone is within a region of low VP and low VP/VS, similar to the observation in northern Honshu (Zhang et al., 2004), Hikurangi (Reyners et al., 2006), and Hokkaido subduction zones (Nakajima et al., 2009), suggesting a low degree of slab mantle dehydration with fluids in isolated pockets. In the Shumagin segment to the northeast (cross-section b-b’ in Figures 3.8 and 3.9), slab dehydration extends to ~80 km depth and ~10-20 km beneath the slab surface, consistent with the distribution of intermediate-depth earthquakes. In the Chignik and Chirikof segments further to the northeast (cross-sections c-c’ and d-d’ in Figures 3.8 and 3.9), slab dehydration is not obvious from ~15 km to 40 km depths, implying a relatively dry plate interface at the seismogenic depths (Wang et al., 2024). The significantly lower number of earthquakes at 15~40 km also indicates a low degree of slab dehydration in a 98 relatively dry slab. A thin (5–10 km thick) moderately high VP/VS layer is found at 40-60 km depth where the slab intersects with the overriding plate Moho, probably related to the dehydration of sediments. Previous active-source imaging reports a thick sediment layer subducted in the Chignik segment (Li et al., 2018). Slab dehydration observed at 40-60 km in this region is consistent with numerical simulations of sediment dehydration in a relatively cold subduction setting (Abers et al., 2017). We also observe low-degree dehydration in the slab crust at ~80–110 km and ~120–150 km depths in the Chignik and Chirikof segments, respectively. Besides, the overall low level of intermediate-depth seismic activity also indicates that slab dehydration is less intensive in the Chignik and Chirikof segments compared to the other segments (Wei et al., 2021). In the Kodiak segment to the northeast end (cross-sections e-e’ and f-f’ in Figures 3.8 and 3.9), slab dehydration extends to at least 150 km depth and 10-40 km beneath the slab surface, leading to abundant intermediate-depth earthquakes in the slab crust and mantle (Shillington et al., 2015; Wei et al., 2021). The strong along-strike variations of slab dehydration at depths are related to the varying degree of hydration of the incoming Pacific Plate at the outer rise. In the Pavlof and Shumagin segments to the southwest, active-source seismic imaging and surface-wave tomography have reported a highly hydrated incoming plate near the trench where the preexisting fabric orientation has a low angle with the trench (Feng, 2021; Li et al., 2024; Shillington et al., 2015). Abundant normal faults at the outer rise allow deep hydration of the slab crust and mantle. Consequently, the highly hydrated slab undergoes extensive dehydration to 200 km depth in the Pavlof and Shumagin segments. In the Chignik and Chirikof segments to the northeast, the incoming plate is less hydrated because the seafloor fabric orientation has a high angle with the trench. Much fewer normal faults are developed and thick sediments overlay the slab crust 99 (Shillington et al., 2015). Thus, slab dehydration at depths below 50 km is limited due to a relatively dry slab. In the Kodiak segment to the northeastern end, the incoming plate is extensively hydrated (Wang et al., 2024), probably related to the subduction of the Kodiak- Bowie Ridge seamounts and fracture zones (von Huene et al., 2021) and the relatively low angle (30–40˚) between the trench and seafloor fabrics. The highly hydrated slab leads to extensive dehydration beyond 100 km depth. In summary, we observe strong along-strike variations in slab dehydration beneath the Alaska Peninsula, which correlate with the changes in the hydration state of the incoming Pacific Plate. This remarkable correlation suggests that the along-strike variations in the slab hydration state can persist from the trench to great depths. The slab gets hydrated in different degrees at the outer rise due to the changes in the abundance of normal faults, sediments, seamounts, and fracture zones in the incoming plate. After subduction, slab dehydration at depths can critically affect the fluid distribution in the plate interface, subducting slab, and overriding plate, facilitating large variations in the slab slip behavior at the seismogenic zone and intermediate- depth seismicity (Wang et al., 2024; Wei et al., 2021). 3.4.3 Sub-arc melting processes Our VP and VP/VS models also provide constraints on sub-arc melting processes in the mantle wedge and the overriding plate. Low VP and high VP/VS anomalies in the sub-arc imply partial melt or serpentinization whereas high VP and high VP/VS anomalies in the sub-arc crust may be caused by volcanic intrusive rocks along with potential magmatic conduits (Nakajima & Hasegawa, 2003; Nugraha et al., 2019). In other cases, low VP and low VP/VS anomalies in the mantle wedge may indicate high temperature and strong anisotropy of peridotite with oriented minerals (Hacker & Abers, 2012), whereas high VP and low VP/VS anomalies in the lower crust 100 may indicate ultramafic rocks (Shillington et al., 2013). We observe large along-strike variations of VP and VP/VS structures in the mantle wedge and the sub-arc crust (Figures 3.9, 3.10, 3.14). In the sub-arc crust, low VP (~6.0 km/s, more than 0.2 km/s lower than the mean VP at the same depth) and high VP/VS (> 1.82) anomalies are located right beneath the arc volcanoes, implying melting processes in the continental crust. Most of the low VP and high VP/VS anomalies extend into the mantle wedge, indicating the deep origin of arc magmatism. The VP/VS value in the mantle wedge is likely to be underestimated due to seismic anisotropy (Hacker & Abers, 2012). In the Pavlof segment to the southwest, thick (>40 km) and widely (~150 km) distributed low VP and high VP/VS anomalies are found in the mantle wedge and lower crust of the overriding plate (Figure 3.14), suggesting extensive melting processes in the sub-arc. It is consistent with the Pavlof volcano being among the most active volcanoes in the Alaska Peninsula (Global Volcanism Program, 2024). The extensive sub-arc melting provides abundant magma that ascends and erupts through Pavlof Volcano. The corner of the mantle wedge seaward of the arc is also characterized by low VP and high VP/VS, unlike observations in several other subduction zones (Abers et al., 2017), indicating that the fore-arc mantle wedge corner is serpentinized by the fluids released from the slab above 80 km depth. Besides, our VP/VS model indicates extensive sub-arc flux melting extends to ~100 km beneath the Pavlof volcano, consistent with a strong deep fluid signature reported by previous geochemical analysis (Wei et al., 2021). In the Shumagin segment, low VP and moderately high VP/VS anomalies are observed immediately above the slab surface, probably indicating a moderate degree of flux-melting in the mantle wedge. In the Chignik and Chirikof segments, despite the widespread low VP anomalies in the mantle wedge, high VP/VS anomalies are confined at 80-140 km depths immediately above 101 the slab. These anomalies suggest a hot mantle wedge with limited flux melting due to a limited fluid supply from slab dehydration, consistent with lower Ph/Th signatures at Aniakchak Volcano and Ukinrek Maars (Wei et al., 2021). Although the low VP/VS value (~1.72) at the fore- arc mantle wedge corner may be underestimated due to anisotropy (Hacker & Abers, 2012), this corner appears to be less serpentinized in the Chignik and Chirikof segments than in the Pavlof segment. Low VP and moderately high VP/VS anomalies are also observed in the overriding plate crust immediately beneath Aniakchak and Ukinrek, indicating partial melts feeding these volcanoes. In the Kodiak segment to the northeastern end, the mantle wedge is characterized by low VP and moderately high VP/VS anomalies, indicating a high temperature and moderate flux melting. Strong high VP/VS anomalies (up to 1.9) are confined in the lower crust (Figure 3.14), indicating a widespread partial molten mushy region with a volume of 45000–80000 km3 in the lower crust beneath volcanoes from Mount Douglas to the northeast to Mount Martin to the southwest. This broad mush zone beneath Trident is consistent with geophysical observations beneath Mount St. Helens (Bedrosian et al., 2018) and Long Valley (Biondi et al., 2023) and the conceptual model of deep crustal hot zones based on petrological studies (Annen et al, 2005). The upper crust of the overriding plate has a confined low VP zone beneath volcanoes, suggesting more silicic magma reservoirs at shallow depths. It is worth noting that the upper crust of this region has the lowest VP/VS (~1.65), similar to the observation in northeastern Japan (Nakajima et al., 2001). This low VP/VS value can only be explained with abundant quartz, which is surprising as all volcanoes in this region consist of andesites or basaltic andesites (Global Volcanism Program, 2024). We suspect that the extremely low VP/VS value reflects strong anisotropy in the upper crust, which may be related to the compressional stress state restricting the magma rising towards the upper crust (Larsen, 2016). 102 3.5. Conclusions In this study, we systematically investigate the slab morphology, slab dehydration, and its relationship with the sub-arc melting in the Alaska Peninsula section of the Alaska-Aleutian subduction zone. We find that the Pacific Plate subducts to ~660 km depth with strong along- strike variations in the slab morphology below 200 km. The Alaska slab has several >100-km- wide vertical gaps in the Pavlof, Shumagin, Chignik, and Kodiak segments, while it subducts continuously to the bottom of the MTZ in the Chirikof segment with an abrupt increase in the slab dip angle at 200 km depth. According to plate reconstruction models (e.g., Fuston & Wu, 2020), the high-VP slab structures in the MTZ are parts of the subducted Pacific Plate. We suspect that the Kula Plate may have subducted to the lower mantle or outside of our study area. We also observe strong along-strike variations in slab dehydration down to 200 km depth and sub-arc melting in the mantle wedge and overriding plate. We propose that the along-strike changes in sub-arc melting are mainly controlled by the changes in slab dehydration, which is further linked to the variations of slab hydration at the outer rise in different segments (Figure 3.15). In the Pavlof and Shumagin segments where the slab is most hydrated through abundant outer rise normal faults (Shillington et al., 2015), slab dehydration after subduction is most intensive among all segments and persists to at least 10 km beneath the slab Moho and to ~200 km in depth. The prominent slab dehydration triggers extensive flux melting, generating large volumes of partial melts in the mantle wedge and overriding crust to feed Pavlof and Veniaminof volcanos. In the Chignik and Chirikof segments, the Alaska slab is relatively dry due to a low degree of hydration at the outer rise and a large volume of sediment subduction (Li et al., 2015; Li & Freymueller, 2018; Shillington et al., 2015). Consequently, slab dehydration is less intensive than in the Pavlof and Shumagin segments, inducing a lower degree of sub-arc melting 103 beneath Aniakchak volcano and Ukinrek Maars. In the Kodiak segment to the northeast end, the Alaska slab is moderately hydrated due to the subduction of the seamounts and fracture zones (von Huene et al., 2021). Extensive slab dehydration extends ~10 km in the slab mantle and to ~120-140 km in depth. Intensive flux melting triggers the widely spread sub-arc melting in the mantle wedge and the overriding crust beneath Trident and Augustine volcanoes. 104 FIGURES Figure 3.1. Bathymetry and topography of the Alaska Peninsula and adjacent region. The orange curve with ticks indicates the Alaska-Aleutian Trench where the Pacific Plate is subducting into the deep Earth at a convergence rate of 63 mm/yr (Bird, 2003). The magenta curves illustrate the slab surface (top of the plate interface) at 20, 40, 60, 80, and 100 km depths. The slab surface geometry is based on a model constrained by active-source seismic experiments (Kuehn, 2019) combined with the Slab2 model (Hayes et al., 2018). Black dashed contours outline the rupture zones of historical megathrust earthquakes determined from their aftershock distributions (Davies et al., 1981), whereas black solid contours show 1-m slip areas of previous megathrust earthquakes (Elliott et al., 2022; Freymueller et al., 2021; Xiao et al., 2021). Red 105 Figure 3.1 (cont’d) triangles indicate Holocene volcanoes (Global Volcanism Program, 2024), with the names of 6 of them labeled. The gray dashed lines and blue bars illustrate the 5 segments discussed in this paper, following Wei et al. (2021) and Wang et al. (2024). The red and black box indicate the area where large-scale teleTomoDD and regional-scale tomoDD are performed, respectively. 106 Figure 3.2. Event and station distributions used in this study. (a and b) Distribution of teleseismic events (red dots) and global stations (blue crosses) used in the first step of tomography for a large-scale velocity model. (c) Distribution of regional stations (blue triangles) and local events (red dots) used in the second step of tomography for high-resolution VP and VP/VS models. Other features are the same as in Figure 3.1. 107 Figure 3.3. Maps of the large-scale VP model (a) at 50, 150, 200, 300, 400, 500, and 600 km depths, and VS model (b) at 50, 150, and 200 km depths imaged with teletomoDD. Anomalies are referenced to the 1-D Ak135-f model (Kennett, 1995). Bold dark-green lines indicate the cross-sections in Figures 3.4 and 3.6. The black curves show the coastlines. The red contours in (a) and (b) outline 1% of dVP and 0.8% of dVS, respectively. The thick dashed lines 108 Figure 3.3 (cont’d) show the boundary between adjacent segments in this study. The gray dashed lines and blue bars illustrate the 5 segments discussed in this paper (Pav: Pavlof segment; Shu: Shumagin segment; Chi: Chignik segment; Chr: Chirikof segment; Kod: Kodiak segment). 109 Figure 3.4. VP and Vs models imaged with teletomoDD along 6 cross-sections in Figure 3.3. The black curves show the slab surface geometry based on a model constrained by active-source seismic experiments (Kuehn, 2019) combined with the Slab2 model (Hayes et al., 2018). Black 110 Figure 3.4 (cont’d) dots indicate the relocated earthquakes. The low-resolution regions are masked according to the resolution tests. Topography/bathymetry is plotted on the top of each panel with vertical exaggeration. The segment and volcano names are shown atop the topography/bathymetry. 111 Figure 3.5. Checkerboard tests for the VP model imaged with teletomoDD. The input models are checkerboard-shaped velocity anomalies in 3D with dVP varying from -5% to 5%. (a) and (b) show the first resolution test result. The size of the input velocity anomalies is 145, 99, and ~100 in the longitude, latitude, and depth directions, respectively. (c) and (d) show the second resolution test result. The anomaly size decreases to ~68, ~67, and 60 in the longitude, latitude, and depth directions, respectively. Model resolvability contours of 20 (green) in the first checkerboard test are chosen to mask out low-resolution regions. The left, middle, and right column represent the checkerboard results in horizontal, latitude, and longitude directions, respectively. The dashed lines in the left column show the location of cross-sections in the middle and right columns. The black curves and dots in cross-sections show the slab surface and relocated earthquakes, respectively. 112 Figure 3.6. Synthetic tests with different input models for the VP model imaged with teletomoDD. (a) and (b) show the input and recovered models if the slab subducts only to 200 km depth. (c) and (d) show the input and recovered models if the slab subducts to 660 km depth 113 Figure 3.6 (cont’d) as a continuous feature. (e) and (f) show the input and recovered models if the slab subducts to 600 km with a gap at 200-300 km depths. The black curves and white dots in cross-sections show the slab surface and relocated earthquakes, respectively. 114 Figure 3.7. Maps of the high-resolution VP (a) and VP/VS (b) models imaged with regional tomoDD at 30, 50, and 100 km depths. The black dashed lines and blue bars in (a) illustrate the 5 segments discussed in this paper (Pav: Pavlof segment; Shu: Shumagin segment; Chi: Chignik 115 Figure 3.7 (cont’d) segment; Chr: Chirikof segment; Kod: Kodiak segment). Bold darkgreen and black lines in (b) indicate the cross-sections in Figures 3.8, 3.9, 3.11 and 3.14. Dark yellow contours on the Pacific Plate outline seafloor magnetic anomalies (Meyer et al., 2017). Other features are the same in Figure 3.1. 116 Figure 3.8. VP model imaged with tomoDD along 6 cross-sections shown in Figure 3.7. In each panel, the black curve shows the slab surface based on the Slab2 (Hayes et al., 2018) model combined with the geometry changes constrained by previous active-source seismic studies (Kuehn, 2019), where the gray curve shows the slab Moho. The dark-green and red dashed curves show the bottom of the sediment layer and the overriding plate Moho, respectively, based on an independent surface-wave tomography study (Li et al., 2024). Regions with low resolutions are masked. Topography/bathymetry is plotted on the top of each panel with vertical exaggeration. The relocated earthquake distribution is shown as black dots in each cross-section. The VP structure above 50 km is plotted relative to the average VP of the structure above the slab surface at the same depth. 117 Figure 3.9. VP/VS model imaged with tomoDD along 6 cross-sections shown in Figure 3.7. Other features are the same as in Figure 3.8. 118 Figure 3.10. Checkerboard tests for the VP (a, c, d, g, i, j) and VP/VS (b, e, f, h, k, l) models imaged with tomoDD. The input models are checkerboard-shaped velocity anomalies in 3D 119 Figure 3.10 (cont’d) with dVP varying from -5% to 5% and d(VP/VS) varying from -8% to 8%. (a–f) The first resolution test result. The size of the input velocity anomalies is 26, 22, and 10 km in the longitude, latitude, and depth directions, respectively. (g–l) The second resolution test result. The anomaly size increases to 78, 77, and 15 km in the longitude, latitude, and depth directions, respectively. Model resolvability contours of 0.54 (purple) in the second checkerboard test are chosen to mask out low-resolution regions. The gray dashed lines in the 30 km depth map view (a and g) indicate the locations of cross-sections in c-f and i-l. 120 Figure 3.11. Two synthetic tests for a high VP/VS layer with different thicknesses immediately beneath the slab surface. In each panel, the black curve shows the slab surface based on the Slab2 (Hayes et al., 2018) model combined with the geometry changes constrained by previous active-source seismic studies (Kuehn, 2019). (a–b) Input and recovered models assuming the high VP/VS layer is 10 km thick. (c–d) Input and recovered models for a 15-km- thick high VP/VS layer atop the slab and a 10-km thick high VP/VS layer on top of the overriding plate mimicking a hypothetical thick sediment layer. 121 Figure 3.12. Slab morphology from Earth’s surface to the bottom of MTZ. The Alaska slab subducts to 660 km beneath the Kodiak and Chirikof segments, while the slab has a ~ 300 km long vertical gap between the high-velocity slab above and the high-velocity structure at the lower MTZ beneath the Chignik, Shumagin, and Pavlof segments. The segment boundaries are plotted at the topography/bathymetry map on the top. 122 Figure 3.13. The high-resolution VP and VP/VS structure at the slab surface and 10 km beneath the slab surface. Other features are the same as in Figure 3.7. 123 Figure 3.14. Cross-sections of the VP and VP/VS models beneath the volcanic arc. The cross- section (g-g’) location is shown in Figures 3.7 and 3.13. The black dashed, black solid, and purple dashed curves represent the location of the slab Moho, slab surface, and the overriding plate Moho, respectively. The locations of volcanoes are shown as red triangles on the top with the names of selected volcanoes. The VP structure above 50 km is plotted relative the average VP of the structure above the slab surface at the same depth. 124 Figure 3.15. Cartoon for the slab dehydration and sub-arc melting. Red triangles on the top indicate the representative volcanoes in each segment. The blue-shaded color indicates the slab hydration state at various degrees at depths. Red regions indicate partial melts in the mantle wedge and overriding plate revealed by the VP and VP/VS model. 125 Figure 3.16. Maps of VP and VP/VS models using different starting VP/VS values at constant depths of 30, 60, 90, and 120 km. (a) Final VP model. (b) Final VP/VS model inverted using a uniform value of 1.75 as the starting model. (c) VP/VS model inverted using the VP/VS values from the first step as the starting model. Other features are the same as in Figure 3.7. 126 REFERENCES Abers, G. A., van Keken, P. E., & Hacker, B. R. (2017), The cold and relatively dry nature of mantle forearcs in subduction zones, Nat. Geosci., 10(5), 333-337, doi:10.1038/ngeo2922. Annen, C., Blundy, J. D., & Sparks, R. S. J. (2005), The genesis of intermediate and silicic magmas in deep crustal hot zones, J. Petrol, 47(3), 505-539, doi: 10.1093/petrology/egi084. Barcheck, G., Abers, G. A., Adams, A. N., Bécel, A., Collins, J., Gaherty, J. B., et al. (2020), The Alaska amphibious community seismic experiment, Seismol. Res. Lett., 91(6), 3054- 3063, doi:10.1785/0220200189. Bedrosian, P. A., Peacock, J. R., Bowles-Martinez, E., Schultz, A., & Hill, G. J. (2018), Crustal inheritance and a top-down control on arc magmatism at Mount St Helens, Nat Geosci, 11(11), 865-870, doi: 10.1038/s41561-018-0217-2. Billen, M. I. (2008), Modeling the dynamics of subducting slabs, Annu. Rev. Earth Planet. Sci., 36(1), 325-356, doi:10.1146/annurev.earth.36.031207.124129. Biondi, E., Zhu, W., Li, J., Williams, E. F., & Zhan, Z. (2023), An upper-crust lid over the Long Valley magma chamber, Sci. Adv., 9(42), eadi9878, doi:10.1126/sciadv.adi9878. Bird, P. (2003), An updated digital model of plate boundaries, Geochem. Geophy. Geosy., 4(3), doi:10.1029/2001gc000252. Boyce, A., Liddell, M. V., Pugh, S., Brown, J., McMurchie, E., Parsons, A., et al. (2023), A New P‐Wave Tomographic Model (CAP22) for North America: Implications for the Subduction and Cratonic Metasomatic Modification History of Western Canada and Alaska, J. Geophys. Res. Solid Earth, 128(3), doi:10.1029/2022jb025745. Burdick, S., Vernon, F. L., Martynov, V., Eakins, J., Cox, T., Tytell, J. et al. (2017), Model Update May 2016: Upper‐Mantle Heterogeneity beneath North America from Travel‐ Time Tomography with Global and USArray Data, Seismol. Res. Lett., 88(2A), 319-325, doi:10.1785/0220160186. Buurman, H., Nye, C. J., West, M. E., & Cameron C. (2014), Regional controls on volcano seismicity along the Aleutian arc, Geochem. Geophy. Geosy., 15(4), 1147-1163, doi: 10.1002/2013gc005101. Cai, C., & Wiens, D. A. (2016), Dynamic triggering of deep earthquakes within a fossil slab, Geophys. Res. Lett., 43(18), 9492-9499, doi:10.1002/2016gl070347. Cashman, K. V., Sparks, R. S. J., & Blundy, J. D. (2017), Vertically extensive and unstable magmatic systems: A unified view of igneous processes, Science, 355(6331), eaag3055, 127 doi:10.1126/science.aag3055. Chen, Wang-Ping., & Michael, B. R. (2001), Evidence for a large-scale remnant of subducted lithosphere beneath Fiji, Science, 292(5526), 2475-2479. Cordell, D., Naif, S., Evans, R., Key, K., Constable, S., Shillington, D., & Bécel, A. (2023), Forearc seismogenesis in a weakly coupled subduction zone influenced by slab mantle fluids, Nat. Geosci., doi:10.1038/s41561-023-01260-w. Dahm, H. H., Gao, S. S., Kong, F., & Liu, K. H. (2017), Topography of the Mantle Transition Zone Discontinuities Beneath Alaska and Its Geodynamic Implications: Constraints From Receiver Function Stacking, J. Geophys. Res. Solid Earth, 122(12), doi:10.1002/2017jb014604. Davies, J., Sykes, L., House, L., & Jacob, K. (1981), Shumagin seismic gap, Alaska Peninsula: History of great earthquakes, tectonic setting, and evidence for high seismic potential, J. Geophys. Res., 86(B5), 3821-3855, doi:10.1029/JB086iB05p03821. Drooff, C., & Freymueller, J. T. (2021), New Constraints on Slip Deficit on the Aleutian Megathrust and Inflation at Mt. Veniaminof, Alaska From Repeat GPS Measurements, Geophys. Res. Lett., 48(4), e2020GL091787, doi:10.1029/2020gl091787. Eberhart-Phillips, D., Christensen, D. H., Brocher, T. M., Hansen, R., Ruppert, N. A., Haeussler, P. J., et al. (2006), Imaging the transition from Aleutian subduction to Yakutat collision in central Alaska, with local earthquakes and active source data, J. Geophys. Res. Solid Earth, 111(B11). doi:10.1029/2005jb004240. Elliott, J. L., Grapenthin, R., Parameswaran, R. M., Xiao, Z., Freymueller, J. T., & Fusso, L. (2022), Cascading rupture of a megathrust, Sci. Adv., 8(18), eabm4131, doi:10.1126/sciadv.abm4131. Faccenda, M. (2014), Water in the slab: A trilogy, Tectonophysics, 614, 1-30, doi:10.1016/j.tecto.2013.12.020. Feng, L. (2021), Amphibious shear wave structure beneath the Alaska‐Aleutian subduction zone from ambient noise tomography, Geochem. Geophys. Geosyst., 22(5), doi:10.1029/2020gc009438. Freymueller, J. T., Suleimani, E. N., & Nicolsky, D. J. (2021), Constraints on the slip distribution of the 1938 Mw 8.3 Alaska Peninsula earthquake from tsunami modeling, Geophys. Res. Lett., 48(9), e2021GL092812, doi:10.1029/2021gl092812. Fukao, Y., & Obayashi, M. (2013), Subducted slabs stagnant above, penetrating through, and trapped below the 660 km discontinuity, J. Geophys. Res. Solid Earth, 118(11), 5920- 5938, doi:10.1002/2013jb010466. 128 Fuston, S., & Wu, J. (2020), Raising the Resurrection plate from an unfolded-slab plate tectonic reconstruction of northwestern North America since early Cenozoic time, GSA Bulletin, 133(5-6), 1128-1140, doi:10.1130/b35677.1. Global Volcanism Program, 2024. [Database] Volcanoes of the World (v. 5.1.7; 26 Apr 2024). Distributed by Smithsonian Institution, compiled by Venzke, E. doi:10.5479/si.GVP.VOTW5-2023.5.1 Goes, S., Agrusta, R., van Hunen, J., & Garel, F. (2017), Subduction-transition zone interaction: A review, Geosphere, 13(3), 644-664, doi:10.1130/ges01476.1. Gou, T., Zhao, D., Huang, Z., & Wang, L. (2019), Aseismic Deep Slab and Mantle Flow Beneath Alaska: Insight From Anisotropic Tomography, J. Geophys. Res. Solid Earth, 124(2), 1700-1724, doi:10.1029/2018jb016639. Grove, T. L., Till, C. B., Lev, E., Chatterjee, N., & Medard, E. (2009), Kinematic variables and water transport control the formation and location of arc volcanoes, Nature, 459(7247), 694-697, doi:10.1038/nature08044. Guo, H., Zhang, H., & Froment, B. (2018), Structural control on earthquake behaviors revealed by high-resolution Vp/Vs imaging along the Gofar transform fault, East Pacific Rise, Earth Planet. Sci. Lett., 499, 243-255, doi:10.1016/j.epsl.2018.07.037. Hacker, B. R., & Abers, G. A. (2012), Subduction Factory 5: Unusually low Poisson's ratios in subduction zones from elastic anisotropy of peridotite, J. Geophys. Res. Solid Earth, 117, B06308, doi:10.1029/2012jb009187. Hacker, B. R., Peacock, S. M., Abers, G. A., & Holloway, S. D. (2003), Subduction factory 2. Are intermediate-depth earthquakes in subducting slabs linked to metamorphic dehydration reactions?, J. Geophys. Res. Solid Earth, 108(B1), doi:10.1029/2001jb001129. Hao, S., Shearer, P., & Liu, T. (2023), The Upper‐Mantle Structure Beneath Alaska Imaged by Teleseismic S‐Wave Reverberations, J. Geophys. Res. Solid Earth, 128(6), doi:10.1029/2023jb026667. Hayes, G. P., Moore, G. L. P., Portner D. E., Hearne, M. F., Flamme, H., Furtney, M., et al. (2018), Slab2, a comprehensive subduction zone geometry model, Science, 362(6410), 58-61, doi:10.1126/science.aat4723. Jia, Z., Shen, Z., Zhan, Z., Li, C., Peng, Z., & Gurnis, M. (2020), The 2018 Fiji M 8.2 and 7.9 deep earthquakes: One doublet in two slabs, Earth Planet. Sci. Lett., 531, doi:10.1016/j.epsl.2019.115997. Jiang, C., Schmandt, B., Ward, K. M., Lin, F. C., & Worthington, L. L. (2018), Upper Mantle Seismic Structure of Alaska From Rayleigh and S Wave Tomography, Geophys. Res. 129 Lett., 45(19), doi:10.1029/2018gl079406. Kennett, B. L., Engdahl. R., & Buland, R. (1995), Constraints on seismic velocities in the earth from travel times, Geophys. J. Int, 122(1), 108-124. Kuehn, H. (2019), Along-trench segmentation and down-dip limit of the seismogenic zone at the eastern Alaska-Aleutian subduction zone, 334 pp, Dalhousie University. Larsen, J. F. (2016), Unraveling the diversity in arc volcanic eruption styles: Examples from the Aleutian volcanic arc, Alaska, J. Volcanol. Geotherm. Res, 327, 643-668, doi:10.1016/j.jvolgeores.2016.09.008. Lay, T. (1994), The fate of desending slabs, Annu. Rev. Earth Planet. Sci, 22, 33-61. Li, J., Shillington, D. J., Bécel, A., Nedimović, M. R., Webb, S. C., Saffer, D. M., et al. (2015), Downdip variations in seismic reflection character: Implications for fault structure and seismogenic behavior in the Alaska subduction zone, J. Geophys. Res. Solid Earth, 120(11), 7883-7904, doi:10.1002/2015jb012338. Li, J., Shillington, D. J., Saffer, D. M., Bécel, A., Nedimović, M. R., Kuehn, H., et al. (2018), Connections between subducted sediment, pore-fluid pressure, and earthquake behavior along the Alaska megathrust, Geology, 46(4), 299-302, doi:10.1130/g39557.1. Li, S., & Freymueller, J. T. (2018), Spatial variation of slip behavior beneath the Alaska Peninsula along Alaska-Aleutian subduction zone, Geophys. Res. Lett., 45(8), 3453-3460, doi:10.1002/2017gl076761. Li, Z., Wiens, D. A., Shen, W., & Shillington, D. J. (2024), Along-strike variations of Alaska subduction zone structure and hydration determined from amphibious seismic data, J. Geophys. Res. Solid Earth, e2023JB027800, doi:10.1029/2023JB027800. Liu, C., Zhang, S., Sheehan, A. F., & Ritzwoller, M. H. (2022), Surface wave isotropic and azimuthally anisotropic dispersion across Alaska and the Alaska‐Aleutian subduction zone, J. Geophys. Res. Solid Earth, 127(11), e2022JB024885, doi:10.1029/2022jb024885. Lu, C., Grand, S. P., Lai, H., & Garnero, E. J. (2019), TX2019slab: A new P and S tomography model incorporating subducting slabs, J. Geophys. Res. Solid Earth, 124(11), 11549- 11567, doi:10.1029/2019jb017448. Lynner, C. (2021), Anisotropy-revealed change in hydration along the Alaska subduction zone, Geology, 49(9), 1122-1125, doi:10.1130/g48860.1. Madsen, J. K., Thorkelson, D. J., Friedman, R. M., & Marshall, D. D. (2006), Cenozoic to Recent plate configurations in the Pacific Basin: Ridge subduction and slab window magmatism in western North America, Geosphere, 2(1), doi:10.1130/ges00020.1. 130 Martin‐Short, R., Allen, R., Bastow, I. D., Porritt, R. W., & Miller, M. S. (2018), Seismic Imaging of the Alaska Subduction Zone: Implications for Slab Geometry and Volcanism, Geochem. Geophys. Geosyst., 19(11), 4541-4560, doi:10.1029/2018gc007962. Meyer, B., Chulliat, A., & Saltus, R. (2017). Derivation and error analysis of the earth magnetic anomaly grid at 2 arc min resolution version 3 (EMAG2v3). Geochem. Geophys. Geosyst, 18, 4522–4537, doi: 10.1002/2017GC007280 Miller, P. K., Saffer, D. M., Abers, G. A., Shillington, D. J., Bécel, A., Li, J., et al. (2021), P‐ and S‐wave velocities of exhumed metasediments from the Alaskan subduction zone: Implications for the in situ conditions along the megathrust, Geophys. Res. Lett., 48(20), e2021GL094511, doi:10.1029/2021gl094511. Müller, R. D., Seton, M., Zahirovic, S., Williams, S. E., Matthews, K. J., Wright, N. M., et al. (2016), Ocean basin evolution and global-scale plate reorganization events since Pangea breakup, Annu. Rev. Earth Planet. Sci., 44(1), 107-138, doi:10.1146/annurev-earth- 060115-012211. Nakajima, J., & Hasegawa, A. (2003), Tomographic imaging of seismic velocity structure in and around the Onikobe volcanic area, northeastern Japan: implications for fluid distribution, J. Volcanol. Geotherm. Res, 127(1-2), 1-18, doi:10.1016/s0377-0273(03)00155-0. Nakajima, J., Matsuzawa, T., Hasegawa, A., & Zhao, D. (2001), Three‐dimensional structure of VP, VS, and VP/VS beneath northeastern Japan: Implications for arc magmatism and fluids, J. Geophys. Res. Solid Earth, 106(B10), 21843-21857, doi:10.1029/2000jb000008. Nakajima, J., Tsuji, Y., Hasegawa, A., Kita, S., Okada, T., & Matsuzawa, T. (2009), Tomographic imaging of hydrated crust and mantle in the subducting Pacific slab beneath Hokkaido, Japan: Evidence for dehydration embrittlement as a cause of intraslab earthquakes, Gondwana Research, 16(3-4), 470-481, doi:10.1016/j.gr.2008.12.010. Nugraha, A. D., Indrastuti, N., Kusnandar, R., Gunawan, H., McCausland, W., Aulia, A. N., et al. (2019), Joint 3-D tomographic imaging of Vp, Vs and Vp/Vs and hypocenter relocation at Sinabung volcano, Indonesia from November to December 2013, J. Volcanol. Geotherm. Res, 382, 210-223, doi:10.1016/j.jvolgeores.2017.09.018. Obayashi, M., Fukao, Y., & Yoshimitsu, J. (2017), Unusually deep Bonin earthquake of 30 May 2015: A precursory signal to slab penetration?, Earth Planet. Sci. Lett., 459, 221-226, doi:10.1016/j.epsl.2016.11.019. Pavlis, G., Jadamec, M., Mann, M. E., Yang, X., Schaeffer, A. W., Wei, S. S, et al. (in press), Synthesis of the seismic structure of the greater Alaska region: Subducting slab geometry, in Tectonics and seismicity of Alaska and Western Canada: EarthScope and beyond, edited by N. Ruppert, M. Jadamec, J. Freymueller, American Geophysical Union. 131 Peacock, S. M., Christensen, N. I., Bostock, M. G., & Audet, P. (2011), High pore pressures and porosity at 35 km depth in the Cascadia subduction zone, Geology, 39(5), 471-474, doi:10.1130/g31649.1. Pesicek, J. D., Thurber, C. H., Zhang, H., DeShon, H. R., Engdahl, E. R., & Widiyantoro, S. (2010), Teleseismic double-difference relocation of earthquakes along the Sumatra- Andaman subduction zone using a 3-D model, J. Geophys. Res., 115(B10), doi:10.1029/2010jb007443. Plafker, G., Moore, J. C., Winkler, G. R., Plafker, G., & Berg, H. C. (1994), Geology of the southern Alaska margin, in The Geology of Alaska, edited, p. 0, Geological Society of America, doi:10.1130/dnag-gna-g1.389. Qi, C., Zhao, D., & Chen, Y. (2007), Search for deep slab segments under Alaska, Phys. Earth Planet. Inter., 165(1-2), 68-82, doi:10.1016/j.pepi.2007.08.004. Rawlinson, N., & Kennett, B. L. N. (2004), Rapid estimation of relative and absolute delay times across a network by adaptive stacking, Geophys. J. Int., 157(1), 332-340, doi:10.1111/j.1365-246X.2004.02188.x. Reyners, M., Eberhart-Phillips, D., Stuart, G., & Nishimura, Y. (2006), Imaging subduction from the trench to 300 km depth beneath the central North Island, New Zealand, with Vp and Vp/Vs, Geophys. J. Int., 165(2), 565-583, doi:10.1111/j.1365-246X.2006.02897.x. Rioux, M., Mattinson, J., Hacker, B., Kelemen, P., Blusztajn, J., Hanghøj, K., et al. (2010), Intermediate to felsic middle crust in the accreted Talkeetna arc, the Alaska Peninsula and Kodiak Island, Alaska: An analogue for low‐velocity middle crust in modern arcs, Tectonics, 29(3), doi:10.1029/2009tc002541. Ruppert, N. A., Barcheck, G., & Abers, G. A. (2022), Enhanced regional earthquake catalog with Alaska Amphibious Community Seismic Experiment data, Seismol. Res. Lett., 94(1), 522-530, doi:10.1785/0220220226. Shephard, G. E., Matthews, K. J., Hosseini, K., & Domeier, M. (2017), On the consistency of seismically imaged lower mantle slabs, Sci Rep, 7(1), 10976, doi:10.1038/s41598-017- 11039-w. Shillington, D. J., Bécel, A., Nedimović, M. R., Kuehn, H., Webb, S. C., Abers, G. A., et al. (2015), Link between plate fabric, hydration and subduction zone seismicity in Alaska, Nat. Geosci., 8(12), 961-964, doi:10.1038/ngeo2586. Shillington, D. J., Van Avendonk, H. J. A., Behn, M. D., Kelemen, P. B., & Jagoutz, O. (2013), Constraints on the composition of the Aleutian arc lower crust from VP/VS, Geophys. Res. Lett., 40(11), 2579-2584, doi:10.1002/grl.50375. Singer, B. S., Jicha, B. R., Leeman, W. P., Rogers, N. W., Thirlwall, M. F., Ryan, J., et al. (2007), 132 Along-strike trace element and isotopic variation in Aleutian Island arc basalt: Subduction melts sediments and dehydrates serpentine, J. Geophys. Res., 112(B6), doi:10.1029/2006jb004897. Takei, Y. (2002), Effect of pore geometry on VP/VS: From equilibrium geometry to crack, J. Geophys. Res., 107(B2), ECV-6, doi:10.1029/2001jb000522. Tian, Y., & Zhao, D. (2012), Seismic anisotropy and heterogeneity in the Alaska subduction zone, Geophys. J. Int., 190(1), 629-649, doi:10.1111/j.1365-246X.2012.05512.x. van Keken, P. E., Hacker, B. R., Syracuse, E. M., & Abers, G. A. (2011), Subduction factory: 4. Depth-dependent flux of H2O from subducting slabs worldwide, J. Geophys. Res., 116(B1), doi:10.1029/2010jb007922. van Stiphout, A. M., Cottaar, S., & Deuss, A. (2019), Receiver function mapping of mantle transition zone discontinuities beneath Alaska using scaled 3-D velocity corrections, Geophys. J. Int., 219(2), 1432-1446, doi:10.1093/gji/ggz360. von Huene, R., Miller, J. J., & Weinrebe, W. (2012), Subducting plate geology in three great earthquake ruptures of the western Alaska margin, Kodiak to Unimak, Geosphere, 8(3), 628-644, doi: 10.1130/ges00715.1. von Huene, R., Miller, J. J., & Krabbenhoeft, A. (2021), The Alaska convergent margin backstop splay fault zone, a potential large tsunami generator between the frontal prism and continental framework, Geochem. Geophys. Geosyst., 22(1), e2019GC008901, doi:10.1029/2019gc008901. Wang, F., Wei, S. S., Drooff, C., Elliott, J. L., Freymueller, J. T., Ruppert, N. A., et al. (2024), Fluids control along-strike variations in the Alaska megathrust slip, Earth Planet. Sci. Lett., 633, 118655, doi: 10.1016/j.epsl.2024.118655. Wang, X. Q., Schubnel, A., Fortin, J., David, E. C., Guéguen, Y., & Ge, H. K. (2012), High Vp/Vs ratio: Saturated cracks or anisotropy effects?, Geophys. Res. Lett., 39(L11307), L11307, doi:10.1029/2012gl051742. Wei, S. S., Shearer, P. M., Lithgow-Bertelloni, C., Stixrude, L., & Tian D. (2020), Oceanic plateau of the Hawaiian mantle plume head subducted to the uppermost lower mantle, Science, 370(6519), 983-987, doi: 10.1126/science.abd0312. Wei, S. S., Ruprecht, P., Gable, S. L., Huggins, E. G., Ruppert, N., Gao, L., et al. (2021), Along- strike variations in intermediate-depth seismicity and arc magmatism along the Alaska Peninsula, Earth Planet. Sci. Lett., 563, 116878, doi:10.1016/j.epsl.2021.116878. Xiao, Z., Freymueller, J. T., Grapenthin, R., Elliott, J. L., Drooff, C., & Fusso, L. (2021), The deep Shumagin gap filled: Kinematic rupture model and slip budget analysis of the 2020 Mw 7.8 Simeonof earthquake constrained by GNSS, global seismic waveforms, and 133 floating InSAR, Earth Planet. Sci. Lett., 576, 117241, doi:10.1016/j.epsl.2021.117241. Yamasaki, T., & Seno, T. (2003), Double seismic zone and dehydration embrittlement of the subducting slab, J. Geophys. Res. Solid Earth, 108(B4), doi:10.1029/2002jb001918. Yang, X., & Gao, H. (2020), Segmentation of the Aleutian‐Alaska Subduction Zone Revealed by Full‐Wave Ambient Noise Tomography: Implications for the Along‐Strike Variation of Volcanism, J. Geophys. Res. Solid Earth, 125(11), doi:10.1029/2020jb019677. Ye, L., Lay, T., Zhan, Z., Kanamori, H., & Hao, J.-L. (2016), The isolated 680 km deep 30 May 2015 MW 7.9 Ogasawara (Bonin) Islands earthquake, Earth Planet. Sci. Lett., 433, 169- 179, doi:10.1016/j.epsl.2015.10.049. ∼ Zelt, C. A. (1998), Lateral velocity resolution from three-dimensional seismic refraction data, Geophys. J. Int., 135(3), 1101-1112, doi:10.1046/j.1365-246X.1998.00695.x. Zhang, H., & Thurber, C. H. (2003), Double-difference tomography: The method and its application to the Hayward fault, California, Bull. Seismol. Soc. Am., 93(5), 1875-1889, doi:10.1785/0120020190. Zhang, H., Thurber, C. H., Shelly, D., Ide, S., Beroza, G. C., & Hasegawa, A. (2004), High- resolution subducting-slab structure beneath northern Honshu, Japan, revealed by double- difference tomography, Geology, 32(4), 361-364, doi:10.1130/G20261.1. Zhang, H., Wang, F., Myhill, R., & Guo, H. (2019), Slab morphology and deformation beneath Izu-Bonin, Nat Commun, 10(1), 1310, doi:10.1038/s41467-019-09279-7. 134 CHAPTER 4: DEEP SLAB MORPHOLOGY CONSTRAINED BY BODY-WAVE TOMOGRAPHY IN THE TONGA SUBDUCTION ZONE An edited version of this chapter is under preparation for the Earth and Planetary Science Letters: Fan Wang, S. Shawn Wei, Douglas A. Wiens, Aubreya Adams, and Haijiang Zhang. Deep slab morphology constrained by body-wave tomography in the Tonga subduction zone 135 Abstract The Tonga subduction zone with the fastest plate convergence speed globally hosts about two- thirds of the world’s deep-focus earthquakes (300-700 km). An active back-arc basin with the fastest back-arc spreading rate in the world has developed as a consequence of rapid trench retreat. The 3D slab morphology in the mantle transition zone remains debated and its influences on deep earthquakes and mantle wedge processes require further investigations. Here we use seismic body-wave arrival-times from both manually picked and machine-learning-powered datasets recorded by regional temporary ocean-bottom seismographs and global land-based stations to image the Tonga subduction zone to 700 km depth. We observe the high-velocity Tonga slab extending into the lower mantle south of 19° S but flattening in the mantle transition zone (MTZ) north of 19° S, with complex slab deformation in the mid-mantle. A remnant slab segment from the fossil Vitiaz Trench is present in the MTZ beneath the Fiji Plateau north of 20° S, with a gap of over 50 km separating it from the Tonga slab to the east. Low-velocity anomalies are imaged in the slab at 350–550 km depths, potentially indicating a metastable olivine wedge or hydrous minerals. In addition, we image sub-slab low-velocity anomalies in the upper mantle, possibly indicating the sub-slab asthenosphere entrained by subduction or the influence of the Samoan mantle plume. Our high-resolution model offers valuable insights into the evolution of the Tonga subduction zone. 136 4.1. Introduction An oceanic plate subducts into the deep mantle at convergence margins due to negative gravitational buoyancy (Fukao and Obayashi, 2013; Goes et al., 2017). A subducted slab goes through complicated geophysical and geochemical processes such as slab rollback and mineral phase transformations, resulting in slab deformation at great depths (Billen, 2020; Guest et al., 2003). Therefore, a slab can display large morphological variations, from flattening at certain depths to penetrating into the lower mantle, according to seismic tomography and geodynamic modeling (Billen, 2020; Fukao and Obayashi, 2013; Goes et al., 2017; Lay, 1994). Small-scale complexities of slab morphology such as bulking, bending, tearing, and the formation of slab holes have also been reported in all subduction zones (Myhill, 2013; Obayashi et al., 2017; Tang et al., 2014; Xi et al., 2024a; Ye et al., 2016; Zhang et al., 2019). Moreover, the subduction of an ancient slab may lead to an isolated fossil slab in the mantle, indicating a long and complex subduction history (Cai and Wiens, 2016; Chen and Brudzinski, 2001; Jia et al., 2020). The Tonga subduction zone marks the fastest convergence margin globally, where the Pacific Plate subducts westward beneath the Indo-Australian Plate along the Tonga Trench at a rate of 240 mm/year in the north (16˚S) and 164 mm/year in the south (21˚S) (Figure 4.1) (Bevis et al., 1995; DeMets et al., 2010). The fast subduction of the Pacific Plate is accompanied by rapid retreat of the Tonga trench at a rate of ~16 cm/yr near the northern end (16˚S) and > 10 cm/yr at 20 ˚S (Bevis et al., 1995; DeMets et al., 2010). Intense seismicity activities have been reported towards the depth at the bottom of MTZ with more than two-thirds of worldwide deep earthquakes (300-700 km) happening in Tonga (Frohlich, 2006; Wei et al., 2017; Xi et al., 2024b). The abundant deep earthquakes in Tonga offer an ideal laboratory to investigate detailed 3D morphological features of the slab at deep depths and their implications for deep earthquake 137 mechanisms. The distribution of deep earthquakes and slab morphology show strong complexities and along-strike variations at the depths of the MTZ in the Tonga subduction zone. Previous studies on the stress state of deep earthquakes show a homogeneous down-dip compressional stress regime in the Tonga slab in the south of 21° S and large variations of the stress orientations between the intermediate and deep portions of the slab in the north of 21° S (Bonnardot et al., 2009). Regional and global seismic tomography studies have reported deflected slab in the MTZ before sinking into the lower mantle in the north, while further south, the slab penetrates the 660 km discontinuity into the lower mantle (Conder and Wiens, 2006; Fukao and Obayashi, 2013; Goes et al., 2017; van der Hilst, 1995). A detached remanent slab from subduction along the fossil Vitiaz trench is reported in the MTZ beneath Fiji by analysis of deep earthquakes and waveform modeling (Cai and Wiens, 2016; Chen and Brudzinski, 2001; Chen and Brudzinski, 2003; Fan et al., 2019; Jia et al., 2020). A recent geodynamic modeling study suggests that the collision of a relic slab from subduction along the Vanuatu trench with the Tonga slab has a major impact on the slab morphology in the MTZ and deep earthquakes in the north Tonga (Liu et al., 2021). In comparison, Peng and Stegman (2024) indicate that the slab kink at the MTZ can be simulated without introducing a relic slab, instead using a hybrid velocity condition and reduced upper mantle viscosity. A joint study of 3D anisotropic tomography and 3D petrological-thermo-mechanical modeling further implies the Samoan plume significantly influences the stagnancy of the Tonga slab in the north and contributes to a slab tear in the south, where the Pacific Plate penetrates the 660 km discontinuity into the lower mantle (Chang et al., 2016). In summary, a flattened slab is observed in northern Tonga while the slab penetrates the bottom of MTZ in the south, but it remains unclear on the location of the relic slabs. The detailed 138 structure of the Tonga slab in the MTZ is not yet fully understood, particularly regarding how a stagnant slab transitions into one that penetrates the lower mantle. In this study, we utilize both regional and teleseismic P- and S-wave arrival-time data to image 3-D VP and VS structures from the Earth’s surface to 700 km depth in the Tonga subduction zone. A large portion of the data comes from multiple temporary deployments of offshore ocean bottom seismographs (OBSs) and land stations in Tonga and Fiji. The new 3-D velocity models provide important constraints on the slab morphology in the upper mantle and melting process beneath the back-arc spreading centers. They also provide valuable insights into the mechanisms of deep earthquakes in the Tonga subduction zone. 4.2. Data and Methods 4.2.1 Data In this study, we assembled regional and teleseismic arrival-times recorded by local stations and global stations, to image 3D P- and S-wave velocity structures of the Tonga subduction zone. The arrival-time data were collected during three key periods—1993-1995, 2001-2002, and 2009-2010—when OBSs and/or land stations were deployed across the Tonga-Lau-Fiji region. From 1993 to 1995, 12 broadband stations were deployed in the islands of Fiji, Tonga, and Niue during the 2-year Southwest Pacific Seismic Experiment (SPASE) (Wiens et al., 1995). A related deployment of 25 OBSs took place in the Lau Basin from September to December 1994 (LABATTS) (Koper et al., 1999; Zhao et al., 1997). In 2001-2002, two arrays of land stations, each consisting of 12 stations, were deployed across the Fiji and Tonga islands as part of the Seismic Arrays in Fiji and Tonga (SAFT) experiment (Tibi and Wiens, 2005). During 2009- 2010, 17 land stations and 49 OBSs were deployed across the Fiji-Lau-Tonga region in the Lau Spreading Center Imaging project (Wei et al., 2016; Zha et al., 2014). We assembled two 139 datasets from these deployments, leveraging all available local stations within the Tonga subduction zone as well as global stations. The first dataset came from the manually picked catalog while the second was generated from a newly developed machine-learning (ML) workflow (Xi et al., 2024b). The primary difference between these two datasets is that the arrival time catalog during 2009-2010 deployment in the ML-based dataset is derived from machine learning picks. For the regional datasets in both the manually picked and ML-based catalogs, we collected P- and S-wave arrival-times of Tonga earthquakes recorded by onshore and offshore seismic stations in the 1993-1995 SPASE, 1994 LABATTS, 2001-2002 SAFT, and 2009-2010 Lau Spreading Center Imaging experiments (Figure. 4.2c). The manually picked dataset included a total of 41,520 P-wave and 17,394 S-wave arrival times. For the ML-based dataset, it consists of manual picks during 1993-1995, 1994 and 2001-2002 deployments and the newly available machine-learning picks during 2009-2010. We selected the P- and S- arrivals from the local events that were shallower than 450 km and recorded by land stations and OBSs from the 2009- 2010 Lau Spreading Center Imaging project (Xi et al., 2024b). These ML picks greatly improved the data coverage in the mantle wedge. To avoid an overrepresentation of deep earthquakes in the Tonga subduction zone, we excluded events deeper than 450 km from the machine learning catalog, thereby reducing dominant sampling paths from the numerous deep events. The number of P- and S-wave arrival-times increased by 9 fold in the ML-based dataset than in the manually picked dataset for the events shallower than 450 km. In total, we assembled 219,090 P- and 61,291 S-wave arrival-times from local events for the regional dataset in our ML-based dataset. We assembled two types of teleseismic arrival-times for both manually picked and ML- based catalog. These two catalogs consist of the same teleseismic arrival times from same local 140 events. The first part included 243,392 P-, 8,522 pP, 9,644 S- and 512 sS-wave arrival times from the Tonga earthquakes recorded by global seismic stations from the International Seismological Centre (ISC) during 1993-1995, 2001-2002, and 2009-2010 (Figure 4.2). We selected the P- and S- arrival-times which followed similar empirical time-distance curves, and excluded the outliners incorrectly reported as P- or S- waves by ISC. The second part included P and S waves from the teleseismic events with magnitudes larger than 6 and epicentral distances from 30° to 90° recorded by regional stations in the 1993-1995 SPASE, 2001-2002 SAFT, and 2009-2010 Lau Spreading Center Imaging experiments. We first downloaded continuous waveforms from the EarthScope Consortium and calculated the spectrum of the seismograms within 20 s before and 120 s after the theoretical P- and S- arrival times. Then, we manually picked the P and S arrival-times from the filtered (0.1-2 Hz for P waves and 0.05-0.3 Hz for S waves) seismograms, and the spectrograms were used for accurate picking when the phase onsets were not clear in the seismograms for some OBSs. The second part consisted of 3,693 P- and 1,507 S wave teleseismic arrival-times that sampled the upper mantle structures in the Tonga subduction zone. In summary, we compiled a total of 288,605 P, 28,545 S, 8,522 pP, and 512 sS arrival-times for the manually picked dataset. For the ML-based dataset, we assembled 446,175 P, 72,442 S, 8,522 pP, and 512 sS arrival-times. 4.2.2 Double-difference tomography for VP and VS We utilize a teleseismic double-difference tomography package teletomoDD (Pesicek et al., 2010) to invert the regional and teleseismic arrival-times for 3D VP and VS models from the Earth surface down to ~700 km depth. Both the manually picked and the ML-based dataset are utilized to invert for 3D velocity models, therefore two velocity models are produced (Figures 4.3, 4.6). 141 The ML-based dataset is used to determine whether incorporating additional arrival times from machine learning methods can improve the resolution of slab and mantle structures. First, we build a regional 3D velocity model nested in a global 3D velocity model to trace the ray paths for both teleseismic and regional arrival-times. The regional model extends from the western end of Fiji (-183.7º E) to the Tonga trench (-172.1º E) and spans from -22.3º S to -15º S with a horizontal spacing of 0.3º and 0.4º in the longitude and latitude directions, respectively. We construct the starting regional velocity model by integrating the top 200 km of the regional Rayleigh-wave model (Wei et al., 2016) with an 8% high-velocity slab model aligning with the seismicity distribution below 200 km, with the ak135-f model (Kennett et al., 1995) as the background model. Incorporating the synthetic high-velocity slab into the starting regional model enhances the robustness of ray tracing for numerous deep earthquakes in Tonga. Additionally, the low-velocity structures in the mantle wedge significantly reduce the average data misfit by approximately 0.9 seconds before inverting for the velocity models. We construct the global model with a much sparser horizontal grid spacing of 5˚ and a vertical spacing of ~90 km below 100 km depth based on the 3D global TX2019slab model (Lu et al., 2019). The 3D global model enables ray tracing of teleseismic waves that travel outside the regional model domain. We perform 6 iterations to refine the 3D velocity models, comprising 4 joint inversions for velocity and relocation, and 2 inversions only for relocation. Due to much smaller number of arrival-times from teleseismic events compared to the regional dataset and teleseismic arrival- times from local events, we first down-weight the regional and teleseismic arrival-times of local events by half and up-weight the arrival-times of teleseismic events by three fold. These adjustments are implemented prior to any inversions to reduce the impact of the imbalanced data volume across the different datasets. During inversions, we down-weight arrival-times with large 142 residuals using a bi-weight function and apply a first-order smoothing constraint to the slowness of model perturbations to stabilize the inversion. The damping and smoothing parameters are chosen with a trade-off analysis of model perturbations and travel-time residuals. The root means square (RMS) value of time residuals decreases from 1.590 s to 0.716 s with the manually picked dataset, whereas the RMS value decreases from 1.620 s to 0.648 s with the ML-based dataset. 4.3. Results 4.3.1 VP and VS models imaged with the manually picked dataset Our VP and VS models imaged with the manually picked dataset show a high-velocity anomaly extending to the MTZ, interpreted as the subducted Tonga slab with obvious along- strike morphological changes (Figures 4.3, 4.4). The top of the high-velocity slab generally coincides with the relocated earthquakes above 400 km depth. This high-velocity anomaly dips into the MTZ with a high angle (~45˚) and flattens in the MTZ north of ~19 ˚S (H2 in Figures 4.3, 4.6). South of ~19 ˚S, the steeply dipping slab penetrates the 660-km discontinuity (H3 in Figures 4.3, 4.6). We also observe high-velocity (~2%-3% higher than the reference ak135-f model) anomalies at ~500-700 km depth beneath the Fiji Plateau (H1 in Figures 4.3, 4.6), separated by a >50 km gap in high-velocity anomalies from the subducted Tonga slab to the east. The high-velocity anomalies to the west have a slightly higher amplitude and seem to be clearer in our VP model than the VS model due to the much more P-wave arrival-times available for seismic imaging (H1 in Figures 4.3, 4.6). Besides, much fewer earthquakes are reported in the vicinity of the high-velocity anomalies to the west than the ones in the flattened Tonga slab where most deep earthquakes are located. Below 250 km depth, we observe a thin layer (<50 km) of low-velocity anomalies (~-1%) within the subducted high-velocity slab at ~450-500 km depth (L1 in Figures 4.3, 4.6) and low- 143 velocity anomalies below the slab at ~200–700 km depths (L2 in Figures 4.3, 4.6), where the high-angle dipping slab starts to flatten to the west in the MTZ (Figure 4.3). The distribution of deep earthquakes also shows a major change from ~45˚ dip to be vertical dip in the vicinity of the low-velocity anomalies in the slab. L1 is most obvious in the cross-section B-B’ from our VP model beneath the dense 2009–2010 OBS array. The amplitude of L1 in our VS model is weaker due to the much fewer S-wave arrival-times available. L2 in the bottom of the slab extends to the bottom of MTZ and spreads more than 300 km along the strike. Sub-slab L2 is primarily located below 200 km depth in the VP model, where the VS model shows lower resolution for sub-slab L2 below ~300-400 km due to fewer S-wave arrival-times available. It is interesting to note that sub- slab L2 is more clearly observed above 400 km depth in the cross-section A-A’ to the north, where the Samoan mantle plume is nearby. Above 250 km depth, we also observe strong low-velocity anomalies L3 at <200 km depths above the westward dipping subducted slab in our VP and VS models (Figure 4.3). These low- velocity anomalies are >10% slower than the reference ak135-f model and extend from ~10 km depth beneath back-arc spreading centers towards deeper depths to the west. Clusters of earthquakes shallower than 50 km depth are observed in the vicinity of the strongest low-velocity anomalies immediately beneath the ELSC and VRF (cross-sections B-B’ and C-C’). Moreover, we find high-velocity anomalies extend from the Earth's surface to ~50 km depth beneath the Fiji Plateau-Lau Ridge, overlaying thick low-velocity anomalies. These thick low-velocity anomalies beneath the Fiji Plateau-Lau Ridge are connected to the strong and shallow low-velocity anomalies beneath the spreading centers to the east. It should be noted that the prominent low- velocity anomalies above 200 km in our VP and VS models mainly come from the initial models built on a Rayleigh-wave tomography model. These low-velocity regions are primarily sampled 144 by parallelly distributed ray paths from deep earthquakes, with fewer contributions from shallower ones, leading to only minor changes in velocities after inversions. However, our 3D seismic imaging with the body-wave arrival-times reveals an increased amplitude of these low- velocity anomalies beneath the ELSC and VRF compared to the Rayleigh-wave velocity model. 4.3.2 VP and VS anomalies imaged with the ML-based dataset In the VP and VS models imaged with the ML-based dataset, we observe similar slab morphology, mantle wedge, and overriding plate structures as those in the models based on the manually picked dataset (Figures 4.6, 4.7). However, major differences are observed in the VP models in the structures near the top of the MTZ between these two datasets. In our VP model from the ML-based dataset, we observe larger low-velocity anomalies (~100 km horizontally and vertically) near the top of the high-velocity slab at 400-500 km depth along the cross-section B- B’ (L1 in Figure 4.6). These low-velocity anomalies exhibit a higher amplitude (~2%) and are more pronounced near where the slab becomes flattened at 400-500 km along the cross-sections A-A’ and C-C’ compared to the manually picked dataset. Besides, large-scale low-velocity anomalies (> 150 km) with amplitudes of approximately -2% are noted above the flattened high- velocity slab beneath the Fiji Plateau at 300-500 km depths in the VP model. Subtle differences are also observed in the mantle wedge structures between the VP models from these two datasets. For example, the thickness of low-velocity anomalies beneath the Fiji Platea-Lau Ridge is reduced by more than 50 km in the VP model from the ML-based dataset. Furthermore, the double seismic zone is more pronounced, primarily due to the greater number of shallow earthquakes in the ML-based dataset. 145 4.3.3 Synthetic tests for model resolution using manually picked dataset and machine learning dataset We perform 4 synthetic tests of VP and VS tomography, including 2 checkerboard tests with the manually picked dataset and two with the ML-based dataset to evaluate the model resolution (Figures 4.5, 4.8). We design two sets of input models with varying sizes of velocity anomalies to evaluate the resolution of models using different datasets (Figure 4.5). The first input velocity anomalies have the size of ~200 km and 200 km in the horizontal and vertical directions, respectively. In our second checkerboard test, we decrease the size of input velocity anomalies to ~66 km, ~88 km, and 80 km in the east-west, north-south, and vertical directions, respectively. We build the synthetic models by adding positive or negative 5% velocity perturbations to the starting VP and VS models. We then compute synthetic P-wave and S-wave travel-times by assuming the same distributions of events and stations as in the real data. Finally, we invert the synthetic data with the same inversion strategy and regularization parameters as in the real data inversion. The checkerboard results for the manually picked dataset show that we can achieve a resolution of ~66 km, ~88 km, and 80 km in the east-west, north-south, and vertical directions, respectively, and a higher resolution in the vicinity of seismicity. By visually comparing the input and recovered models, we choose the number of sampling rays larger than 500 and 100 as the threshold for the model resolvability of the VP and VS models, respectively. Then we use these contours to mask out the low-resolution regions in the final VP and VS models for both the manually picked and ML-based datasets. In comparison, the checkerboard tests using the ML- based dataset show a more accurate recovery of alternating high- and low-velocity anomaly patterns in the slab and beneath the spreading centers than the ones using the manually picked 146 dataset, indicating improved resolution. In summary, the checkerboard results indicate that we can image the slab morphology and velocity anomalies beneath the spreading centers and the Fiji Plateau-Lau at a spatial resolution of ~66 km, ~88 km, and 80 km in the east-west, north-south, and vertical directions, respectively, for both the manually picked and ML-based datasets. The machine learning dataset provides higher resolution due to a more accurate recovery of checkerboard patterns with the same input model. Besides, the VP model shows overall better recovery compared to the VS model, reflecting its higher resolution, likely due to the greater number of P-wave arrival-times compared to S- wave arrival-times in our study. 4.4. Discussions 4.4.1 Deep slab morphology in the MTZ The westward dipping high-velocity (>3%) anomalies are interpreted as the Tonga slab subducting into the MTZ with complicated slab morphology in both our VP and VS models from manually picked and ML-based datasets (Figures 4.3 and 4.6). The slab morphology imaged from these two datasets shares key similarities in the main features discussed below, with differences primarily in the amplitude of the velocity anomalies. Therefore, we visualize the slab morphology in 3D by plotting the 2.5% high-velocity isosurface for structures above 300 km depth and the 4.0% high-velocity isosurface for those below 300 km depth using the VP model from the manually picked dataset (Figure 4.9). The subducted Tonga slab exhibits a distinct change in morphology near ~19º S within the MTZ. North of ~19º S, the high-velocity slab descends with a steep dip angle (~45º) above ~450 km depth, then flattens at around ~550 km depth (Figures 4.3, 4.6, 4.9). In comparison, south of 19º S, the Tonga slab maintains a steep dip angle, descending deeper and reaching the base of the 147 MTZ (Figures 4.3, 4.6, 4.9). The unusually thick (>200 km) high-velocity anomalies in the MTZ along the cross-section B-B’ in Figures 4.3 and 4.6 reveal a transition from the flattened slab (H2) north of 19˚S to the dipping slab (H3) south of 19˚S. We observe a slab tear around 19º S at the bottom of the MTZ, accompanied by an abrupt change in the distribution of deep earthquakes in its vicinity at 500-600 km depth (Figure 4.9). This slab tear appears to develop due to the sharp transition from a flat slab (H2) to a dipping slab (H3) at ~19º S within the MTZ. North of ~18-19º S, the deep earthquakes trend in the NNW direction, while south of 19º S, their trend shifts to the NNE direction (Figures 4.4, 4.7). Notably, these deep earthquakes extend further west beneath the Lau Ridge at ~17.5º S compared to those at ~18-19º S, despite the Tonga Trench being slightly farther west at ~18-19º S (Figures 4.4, 4.7, 4.9). This pattern of deep seismicity is consistent with the presence of a flattened slab to the north and a significant change in slab morphology around 19º S, supporting a slab tear in the MTZ at this latitude (Cai & Wiens, 2016). The complex slab morphology in the MTZ observed in our models aligns with previous regional and global tomography results, while providing higher resolution of the slab structure at 19º S. Regional and global tomography studies have reported a flattened slab within the MTZ near the 410-discontinuity in northern Tonga, transitioning to a steep slab that penetrates into the lower mantle to the south (Amaru, 2007; Fukao and Obayashi, 2013; Lu et al., 2019; van der Hilst, 1995), which is consistent with our findings. Fukao and Obayashi (2013) proposed a slab tear in the MTZ where the Tonga slab intersects with the Kermadec slab. However, our observations place this transition and associated slab tear at approximately 19º S, further north than previously reported (Fukao and Obayashi, 2013; Lu et al., 2019; van der Hilst, 1995). A regional 2D body-wave tomography study shows that the slab morphology transitions from 148 shallow-dipping at around 500 km to nearly vertical at the bottom of MTZ (Conder and Wiens, 2006). We interpret this along-dip transition in slab morphology as resulting from the overlap between the flattened slab to the north of ~19º S and the steeply dipping slab to the south since their 2D profile cuts across the transition region. These overlapped slab segments lead to thickening and morphological changes observed in the MTZ. In their imaging results (Conder and Wiens, 2006), the western extent of the flattened slab is shown as reduced high-velocity anomalies in the MTZ with part of area masked due to lower resolution. Our 3D VP and VS models reveal a clear flattened slab in the north of 19º S and a steeply dipping slab south of 19º S, with an associated slab tear near the transition region. We interpret the complex deep slab morphology as indicative of strong slab deformation in the MTZ, likely caused by non-uniform trench retreat rates and pre-existing structures within the subducted slab. Regions of intense deformation with high strain rates are suggested to experience abundant deep earthquakes, while areas with slow deformation in colder slabs may host fewer earthquakes (Billen, 2020). The observed slab tear and the high frequency of deep earthquakes in its vicinity support this hypothesis, indicating significant deformation of the Tonga slab in the MTZ. North of 19º S, the subducted Tonga slab flattens in the MTZ, likely due to a rapid trench retreat rate. This observation is consistent with previous geodynamic studies, which suggest that a fast trench retreat rate plays a primary role in slab flattening within the MTZ (Christensen, 1996; Li et al., 2019; Peng and Stegman, 2024). South of 19º S, the slab subducts at a steep angle into the lower mantle, likely caused by a reduced trench retreat rate, leading to the formation of a slab tear near potential pre-existing weak zones in the subducted slab. Two significant bathymetric features exist in the subducted oceanic plate between 17º S and 19º S: the Capricorn Seamount (~18.5º S) and the Niue Trough (17-18º S). The Capricorn Seamount is suggested to 149 have formed during the Miocene as part of a partially subducted hotspot chain (Brodie, 1965; Crawford et al., 2003), while the Niue Trough is a 400-km-long bathymetric depression potentially created by rifting-related processes, lithospheric tearing, or plume-trench interactions (Pockalny et al., 2009). Nevertheless, these pre-existing structures are suggested to reduce the rigidity of the ocean plate in their vicinity by four to five orders of magnitude within 100 km of the trench (Arredondo and Billen, 2012). Consequently, the subducted Tonga slab may tear in the vicinity of these weak zones under strong deformation influenced by the non-uniform trench retreat rate. Recent geodynamic modeling also highlights the significant role of seamount subduction in slab tearing (Cheng et al., 2019). Although the precise location of these pre- existing structures at MTZ depth remains uncertain, the clockwise rotation of the Tonga Trench associated slab subduction may place these weak zones near the slab tear where the slab morphology changes (Schellart and Spakman, 2012). We observe high-velocity anomalies extend westward beneath the Fiji Plateau in the MTZ north of 20º S, with a 50-km-long region of neutral to <1% high-velocity anomalies separating them from the eastward high-velocity Tonga slab (H1 in Figures 4.4, 4.7). In contrast to the abundant deep seismicity east of this region where the flattened slab exists, only a few deep earthquakes, including the 2009 Mw 7.3 Fiji earthquake, are reported in H1 (Figures 4.3, 4.6, 4.9). We interpret H1 as a segment of a fossil slab that appears to be isolated and is not attached to the subducted Tonga slab in the MTZ. The origin of the fossil slab beneath the Fiji Plateau remains unclear, and its interpretation varies widely in seismology and geodynamic modeling. Previous seismic studies on the seismicity and focal mechanisms of deep earthquakes (Chen and Brudzinski, 2001; Richards et al., 2011) and geodynamic modeling on the slab morphology (Liu et al., 2021) suggest a 150 detached slab segment from the fossil Vitiaz Trench or Vanuatu Trench beneath the Fiji Plateau. In contrast, recent geodynamic modeling indicates that significant changes in slab morphology at 410 km depth can be simulated without a detached slab segment when considering a rapid trench retreat and a weak upper mantle (Peng and Stegman, 2024). We interpret the isolated remnant slab in the MTZ beneath the Fiji Plateau as a fossil slab from the Vitiaz trench, based on our high-resolution imaging and supporting evidence from plate reconstruction (Figures 4.1, 4.9). According to our 3D VP and VS models, the modeling of Peng and Stegman (2024) may have been influenced by blurred seismic imaging results in previous studies, which likely sample the flattened slab to the north and the steeply dipping slab to the south. Besides, we have a good resolution in the region of the remnant slab, suggesting a distinct fossil slab structure from the subducted Tonga slab (Figures 4.5, 4.8). According to previous plate reconstruction models (Falvey, 1975; Hamburger and Isacks, 1987; Packham, 1973), the subduction of the Pacific Plate beneath the Tonga and Vitiaz Trenches occurred during the Tertiary through the late Miocene when the New Hebrides, Fiji, and Tonga-Lau Ridge were connected. In the late Miocene, the Indo-Australian Plate began to subduct beneath the New Hebrides, accompanied by a clockwise rotation of the New Hebrides and retreat of the Vanuatu and Tonga Trenches (Hamburger and Isacks, 1987). This tectonic change in the late Miocene established the modern-day subduction along the Vanuatu and Tonga Trenches. However, plate reconstructions in the Tonga-Vitiaz-Vanuatu subduction zones carry significant uncertainty due to most of the region being below sea level and the relative sparsity of data (Schellart and Spakman, 2012). Nevertheless, the subduction of the Indo-Australian Plate is much younger than the slab subducted at the Vitiaz Trench. Therefore, the Vanuatu slab is expected to subduct to a shallower depth compared to the Vitiaz slab. Near the Vanuatu Trench, the Wadati-Benioff zone 151 extends to ~300-400 km depth, suggesting the Vanuatu slab subducts to <400 km depth (https://www.isc.ac.uk/isc-ehb/regions/index.php?country=VANU). Additionally, deep outboard earthquakes (>500 km) are observed near the current Vitiaz Trench to the northeast and >300 km away from the seismicity in the Vanuatu slab (Chen and Brudzinski, 2001). The shallow depth limit of seismicity and their significant distance from the deep outboard earthquakes in the Vanuatu subduction zone exclude the possibility of the Vanuatu slab being present in the MTZ beneath the Fiji Plateau. Therefore, we prefer a remnant slab from the fossil Vitiaz Trench beneath Fiji in our VP and VS models. Our 3D velocity models and interpretation of the fossil Vitiaz slab provide substantial justification for the occurrence of large deep earthquakes with distinct characteristics in the Tonga-Fiji region. For example, the Mw 7.9 Tonga-Fiji earthquake on September 6, 2018, occurred near the western edge of the Tonga slab at the bottom of the MTZ, where deep earthquakes are relatively sparse (Figures 4.3, 4.6, 4.9). This Mw 7.9 earthquake has a significantly lower rupture speed (2.5 km/s) than the doublet 19 August 2018 Mw 8.2 earthquake (4.1 km/s) and the nearby 9 November 2009 Mw 7.3 earthquake (4.6 km/s) (Figures 4.1, 4.3, 4.6, 4.9) (Cai and Wiens, 2016; Fan et al., 2019; Jia et al., 2020). It also generated much fewer aftershocks compared to the 9 March 1994 Mw 7.6 and 2018 Mw 8.2 deep earthquakes to its east (Fan et al., 2019; Jia et al., 2020; Wiens and McGuire, 2000). This earthquake is suggested to nucleate in a warm environment, with most of its energy dissipated by rupture melting (Fan et al., 2019; Jia et al., 2020). These distinct characteristics of the 2018 Mw 7.9 earthquake can be satisfactorily explained by our tomography imaging. The Mw 7.9 earthquake is located at the western edge of the Tonga slab, separated by a velocity gap from the Vitiaz slab to the west, where warm thermal condition is expected. Additionally, the aftershocks of the 2009 Mw 7.3 152 deep earthquakes appear to follow a linear feature corresponding to the orientation of the Vitiaz Trench (Cai and Wiens, 2016), supporting the presence of deep earthquakes in the fossil Vitiaz slab in our tomography imaging. Regarding the low-intensity deep earthquakes in the fossil Vitiaz slab, the large distance between the fossil slab and the actively subducting Tonga slab may suggest limited deformation of the Vitiaz slab. Therefore, a lower strain rate near the slowly deformed fossil slab could contribute to a reduced frequency of deep earthquakes (Billen, 2020). 4.4.2 Sub-slab low-velocity anomalies in the MTZ The low-velocity anomalies observed beneath the bottom of the slab at 300-700 km depth in our VP and VS models suggest entrainment of the asthenosphere or the Samoan mantle plume from the north. It is worth noting that our VS model has a lower resolution below 600 km, therefore, we rely on the VP model for the discussion of deep sub-slab low-velocity anomalies (Figures 4.3, 4.4, 4.6, 4.7). The sub-slab low-velocity anomalies may indicate entrained asthenosphere during the fast subduction of the Tonga slab. Previous geodynamic modeling shows that significant volumes (> 100 km thick) of the sub-slab asthenosphere can be entrained into the deep mantle when the plate convergence exceeds ~4 cm/yr (Liu and Zhou, 2015). Given the rapid subduction rate of the Tonga slab (~22 cm/yr) (Bevis et al., 1995), this process likely facilitates the effective entrainment of asthenosphere material into the MTZ. Similar deep entrainment of sub-slab asthenosphere has been observed in other subduction zones, such as the Japan Trench, where the entrained material escapes through slab holes in the Pacific Plate and contributes to magma generation beneath the Changbaishan volcanoes (Huang and Zhao, 2006; Tang et al., 2014). Besides, the large trench-parallel sub-slab anisotropy in the Tonga-Kermadec subduction zone supports the presence of entrained asthenosphere in the deep mantle (Long and Silver, 2009; 153 Song and Kawakatsu, 2012). The sub-slab low-velocity anomalies could also reflect the influence of the Samoan mantle plume from the north. Our VP and VS models show widely spread low-velocity anomalies along the strike, particularly at deep depths of 400-600 km depths (Figures 4.4, 4.7). Previous seismic and geochemistry studies in northern Tonga indicate that materials from the Samoan mantle plume can escape through a slab window at the northern edge of the Tonga slab, significantly affecting the geochemical signatures of volcanic rocks and mantle flow patterns (Falloon et al., 2008; Smith et al., 2001). The deep-originated mantle plume may contribute to the sub-slab low- velocity anomalies we observed in our VP and VS models. However, the path of the Samoan plume remains largely uncertain due to the limited number of seismic stations in the northern region. 4.4.3 Low-velocity anomalies in the slab at 400-500 km depth We observe low-velocity anomalies at 400-500 km depth within the slab and the deep earthquakes in their vicinity (Figures 4.3 and 4.6). These anomalies are thick (~10 km) and more pronounced (~100 km thick and ~2-3%) in the ML-based dataset than in the manually picked dataset (Figures 4.3 and 4.6). Given the fast and cold nature of the subducted Tonga slab, these low-velocity zones may correspond to metastable olivine near the cold core of the slab. However, the narrow width of the metastable olivine may not be resolvable by traditional tomography techniques with earthquake travel times (Koper et al., 1998). Therefore, the ~10 km wide low-velocity anomalies may be the presence of hydrous minerals within the Tonga slab. Thermobarometric models suggest that superhydrous aluminous silica can remain stable at MTZ depths at the coldest subduction zone, such as Tonga (Cerpa et al., 2022). The cold Tonga slab can provide feasible thermal conditions for the hydrous minerals to subduct to significant depths. 154 Previous seismic studies indicate that hydrous minerals can reach ~300 km depth in the Tonga slab based on the intermediate-deep earthquakes in the double-seismic zone and converted waves at the slab surface (Wang and Wei, 2020; Wang et al., 2017). The low-velocity anomalies inside the Tonga slab at 400-500 km may indicate that hydrous minerals can be present at even greater depths in this cold subduction environment, possibly along with metastable olivine. 4.5. Conclusions In this study, we build 3D high-resolution VP and VS models based on two different datasets, one from the manually picked dataset and another from the ML-based catalog. We observe similar morphological features of the Tonga slab in both the MTZ and mantle wedge structures, but the ML-based models reveal more pronounced low-velocity structures at 350-500 km depth above the slab. These low-velocity structures inside the slab potentially indicate the metastable olivine wedge or hydrous minerals at 350-500 km depth. Additionally, the slab morphology exhibits a distinct change in the MTZ at ~19º S. North of ~19º S, the Tonga slab subducts at a steep angle (~45º) into the deep mantle and becomes flattened within the MTZ. South of ~19º S, the steeply dipping slab penetrates the 660 km discontinuity into the lower mantle. This complex slab morphology and associated slab tear suggest significant deformation of the Tonga slab in the MTZ, likely due to ununiform fast trench retreat and preexisting structures in the ocean plate. Beneath the Fiji Plateau north of 20° S, we observe a remnant slab segment in the MTZ with a gap of over 50 km separating it from the Tonga slab to the east. We interpret this remnant segment as a slab from the fossil Vitiaz Trench, where subduction ceased 5-10 million years ago. The high-resolution 3D VP and VS models developed in this study reveal intricate local-scale deep slab morphology, providing critical insights into mantle dynamics within a complex tectonic setting. The coexistence of slab tears with sub-arc low-velocity anomalies—possibly 155 resulting from entrained asthenosphere and mantle plumes—indicates complex interactions between the Tonga slab and mantle flow, which are key to understanding the tectonic evolution of the Southwestern Pacific. 156 FIGURES Figure 4.1. Bathymetry and topography of the Tonga-Lau-Fiji region. The major bathymetric features on the oceanic plate, Capricorn seamount and Niue Tough, are noted with white text on the map. The red inverted triangles show the location of arc volcanoes since Holocene. The black curve with ticks indicates the Tonga Trench where the Pacific Plate is subducting beneath the Indo-Australian Plate. The position of the Vitiaz Trench in the north is also noted by the black curve with ticks. Dark-green solid lines show the locations of cross- sections in Figures 4.3 and 4.6. Magenta contours illustrate the slab surface (top of the plate interface) at 50, 100, 200, 300, 400, 500, and 600 km depths from the Slab2 model (Hayes et al., 2018). Red beach balls show the focal mechanisms of four deep earthquakes in 1994, 2009, and 2018 (Fan et al., 2019; Jia et al., 2020). Black curves in the Lau Basin indicate back-arc 157 Figure 4.1 (cont’d) spreading centers: CLSC, Central Lau Spreading Center; ELSC, East Lau Spreading Center; LETZ, Lau Extensional Transform Zone; FRSC, Fonualei Rift and Spreading Center; FSC, Futuna Spreading Center; MTJ, Mangatolu Triple Junction; VFR, Valu Fa Ridge; and PR, Peggy Ridge. Black contours show bathymetry of 1 km. The inset shows the Tonga-Lau-Fiji region on a global map. 158 Figure 4.2. Event and station distributions used in this study. (a and b) Distribution of teleseismic events (red dots) and global stations (blue crosses) used in seismic tomography. (c 159 Figure 4.2 (cont’d) and d) Distribution of regional stations and local earthquakes in the manually picked and ML- based datasets, respectively. Regional stations from different seismic experiments are labeled with distinct symbols, and earthquakes are color-coded by depth. The other features are the same as in Figure 4.1. 160 Figure 4.3. VP and VS models along 4 cross-sections imaged with the manually picked dataset. The velocity perturbations are plotted relative to the ak135-f model (Kennett et al., 1995). Note that the color scales are different above and below 250 km. In each panel, the black 161 Figure 4.3 (cont’d) curve shows the slab surface based on the Slab2 model (Hayes et al., 2018). Black open circles indicate the relocated earthquakes within a 50-km distance from each cross-section. The low- resolution regions are masked according to the resolution tests. Topography/bathymetry is plotted on top of each panel with vertical exaggeration. The red beach balls show the focal mechanisms of 4 large deep earthquakes in 1994, 2009, and 2018. The H1-H3 and L1-L3 mark the target regions discussed in the main text. 162 Figure 4.4. Maps of the VP model (a) and VS model (b) at 200, 300, 400, 500, 600, and 700 km depths imaged with the manually picked dataset. The velocity perturbations are plotted relative to the ak135-f model (Kennett et al., 1995). The H1-H3 and L1-L3 mark the target regions discussed in the main text. Other features are the same as in Figure 4.1. 163 Figure 4.5. Checkerboard tests of the VP and VS models using the manually picked dataset. The input models are checkerboard-shaped velocity anomalies in 3D with dVP varying from -5% to 5%. (a) and (b) show the first resolution test result. The size of the input velocity anomalies is ~200 km in both horizontal and vertical directions. (c) and (d) show the second resolution test result. The anomaly size decreases to ~66 km, ~88 km, and 80 km in the east-west, north-south, 164 Figure 4.5 (cont’d) and vertical directions, respectively. Model resolvability contours of 500 (magenta) in the first checkerboard test are chosen to mask out low-resolution regions. The left, middle, and right columns represent the checkerboard results in latitude, longitude, and horizontal directions, respectively. The black dots show the relocated earthquakes while the dark-green lines indicate the locations of cross-sections in Figures 4.3 and 4.6. The red beach balls show the focal mechanisms of the large deep earthquakes in Tonga in the corresponding depth. 165 Figure 4.6. VP and VS models along 4 cross-sections imaged with the ML-based dataset. Other features are the same as in Figure 4.3. 166 Figure 4.7. Horizontal maps of slices at different depths imaged with the ML-based dataset. Other features are the same as in Figure 4.4. 167 Figure 4.8. Checkerboard tests of the VP and VS models using the ML-based dataset. Other features are the same as in Figure 4.5. 168 Figure 4.9. 3-D images of the Tonga slab morphology. The isosurfaces for high-velocity (2.5% and 4%) structures are color-coded with blue-to-white for the anomalies above and below 300 kms, respectively, based on depth. Surface volcanoes are marked with red triangles, and other surface features are consistent with those shown in Figure. 4.1. 169 REFERENCES Amaru, M.L., 2007. Global travel time tomography with 3-D reference models. Ph.D. (274), 174. Arredondo, K.M., Billen, M.I., 2012. Rapid weakening of subducting plates from trench-parallel estimates of flexural rigidity. Phys. Earth Planet. Inter. 196-197, 1-13. DOI: https://doi.org/10.1016/j.pepi.2012.02.007 Bevis, M.G., Taylor, F.W., Schutz, B.E., Recy, J., Isacks, B.L., Helu, S., Singh, R., Kendrick, E., Stowell, J., Taylor, B., Calmant, S., 1995. Geodetic observations of very rapid convergence and back-arc extension at the Tonga arc. Nature 374 (6519), 249-251. DOI: https://doi.org/10.1038/374249a0 Billen, M.I., 2020. Deep slab seismicity limited by rate of deformation in the transition zone. Sci. Adv 6 (22), eaaz7692. DOI: https://doi.org/10.1126/sciadv.aaz7692 Bonnardot, M.-A., Régnier, M., Christova, C., Ruellan, E., Tric, E., 2009. Seismological evidence for a slab detachment in the Tonga subduction zone. Tectonophysics 464 (1-4), 84-99. DOI: https://doi.org/10.1016/j.tecto.2008.10.011 Brodie, J.W., 1965. Capricorn seamount, south-west Pacific Ocean. Trans. R. Soc. N. Z. Geol. 3 (10), 151-158. Cai, C., Wiens, D.A., 2016. Dynamic triggering of deep earthquakes within a fossil slab. Geophys. Res. Lett 43 (18), 9492-9499. DOI: https://doi.org/10.1002/2016gl070347 Cerpa, N.G., Arcay, D., Padrón-Navarta, J.A., 2022. Sea-level stability over geological time owing to limited deep subduction of hydrated mantle. Nature Geoscience 15 (5), 423- 428. DOI: https://doi.org/10.1038/s41561-022-00924-3 Chang, S.-J., Ferreira, A.M., Faccenda, M., 2016. Upper- and mid-mantle interaction between the Samoan plume and the Tonga-Kermadec slabs. Nat Commun 7, 10799. DOI: https://doi.org/10.1038/ncomms10799 Chen, W.P., Brudzinski, M.R., 2001. Evidence for a large-scale remnant of subducted lithosphere beneath Fiji. Science 292 (5526), 2475-2479. DOI: https://doi.org/10.1126/science.292.5526.2475 Chen, W.P., Brudzinski, M.R., 2003. Seismic anisotropy in the mantle transition zone beneath Fiji‐Tonga. Geophys. Res. Lett 30 (13). DOI: https://doi.org/10.1029/2002gl016330 Cheng, Z., Ding, W., Faccenda, M., Li, J., Lin, X., Ma, L., Fang, P., Ding, H., 2019. Geodynamic effects of subducted seamount at the Manila Trench: Insights from numerical modeling. Tectonophysics 764, 46-61. DOI: https://doi.org/10.1016/j.tecto.2019.05.011 170 Christensen, U.R., 1996. The influence of trench migration on slab penetration into the lower mantle. Earth Planet. Sci. Lett 140 (1-4), 27-39. Conder, J.A., Wiens, D.A., 2006. Seismic structure beneath the Tonga arc and Lau back‐arc basin determined from joint Vp, Vp/Vs tomography. Geochemistry, Geophys. Geosystems 7 (3). DOI: https://doi.org/10.1029/2005gc001113 Crawford, W.C., Hildebrand, J.A., Dorman, L.M., Webb, S.C., Wiens, D.A., 2003. Tonga Ridge and Lau Basin crustal structure from seismic refraction data. J. Geophys. Res. Solid Earth 108 (B4). DOI: https://doi.org/10.1029/2001jb001435 DeMets, C., Gordon, R.G., Argus, D.F., 2010. Geologically current plate motions. Geophys. J. Int 181 (1), 1-80. DOI: https://doi.org/10.1111/j.1365-246X.2009.04491.x Falloon, T.J., Danyushevsky, L.V., Crawford, A.J., Meffre, S., Woodhead, J.D., Bloomer, S.H., 2008. Boninites and adakites from the northern termination of the Tonga trench: Implications for adakite petrogenesis. Journal of Petrology 49 (4), 697-715. DOI: https://doi.org/10.1093/petrology/egm080 Falvey, D.A., 1975. Arc reversals and a tectonic model for the North Fiji Basin. Bull. Austr. Soc. Explor. Geophys. 6, 47-49. Fan, W., Wei, S.S., Tian, D., McGuire, J.J., Wiens, D.A., 2019. Complex and diverse rupture processes of the 2018 Mw 8.2 and Mw 7.9 Tonga‐Fiji deep earthquakes. Geophys. Res. Lett 46 (5), 2434-2448. DOI: https://doi.org/10.1029/2018gl080997 Frohlich, C., 2006. Deep earthquakes. Cambridge University Press, Cambridge, UK. Fukao, Y., Obayashi, M., 2013. Subducted slabs stagnant above, penetrating through, and trapped below the 660 km discontinuity. J. Geophys. Res. Solid Earth 118 (11), 5920-5938. DOI: https://doi.org/10.1002/2013jb010466 Goes, S., Agrusta, R., van Hunen, J., Garel, F., 2017. Subduction-transition zone interaction: A review. Geosphere 13 (3), 644-664. DOI: https://doi.org/10.1130/ges01476.1 Guest, A., Schubert, G., Gable, C.W., 2003. Stress field in the subducting lithosphere and comparison with deep earthquakes in Tonga. J. Geophys. Res. Solid Earth 108 (B6). DOI: https://doi.org/10.1029/2002jb002161 Hamburger, M.W., Isacks, B.L., 1987. Deep earthquakes in the southwest Pacific: A tectonic interpretation. J. Geophys 92 (B13), 13841-13854. DOI: https://doi.org/10.1029/JB092iB13p13841 Huang, J., Zhao, D., 2006. High-resolution mantle tomography of China and surrounding regions. J. Geophys 111 (B9), B09305. DOI: https://doi.org/10.1029/2005jb004066 171 Jia, Z., Shen, Z., Zhan, Z., Li, C., Peng, Z., Gurnis, M., 2020. The 2018 Fiji M 8.2 and 7.9 deep earthquakes: One doublet in two slabs. Earth Planet. Sci. Lett 531. DOI: https://doi.org/10.1016/j.epsl.2019.115997 Kennett, B.L.N., Engdahl, E.R., Buland, R., 1995. Constraints on seismic velocities in the earth from travel times. Geophys. J. Int 122 (1), 108-124. Koper, K.D., Wiens, D.A., Dorman, L., Hildebrand, J., Webb, S., 1999. Constraints on the origin of slab and mantle wedge anomalies in Tonga from the ratio of S to P velocities. J. Geophys. Res. Solid Earth 104 (B7), 15089-15104. DOI: https://doi.org/10.1029/1999jb900130 Koper, K.D., Wiens, D.A., Dorman, L.M., Hildebrand, J.A., Webb, S.C., 1998. Modeling the Tonga slab: Can travel time data resolve a metastable olivine wedge? J. Geophys. Res. Solid Earth 103 (B12), 30079-30100. DOI: https://doi.org/10.1029/98jb01517 Lay, T., 1994. The fate of desending slabs. Annu. Rev. Earth Planet. Sci 22, 33-61. Li, Z.-H., Gerya, T., Connolly, J.A.D., 2019. Variability of subducting slab morphologies in the mantle transition zone: Insight from petrological-thermomechanical modeling. Earth Sci Rev 196. DOI: https://doi.org/10.1016/j.earscirev.2019.05.018 Liu, H., Gurnis, M., Leng, W., Jia, Z., Zhan, Z., 2021. Tonga slab morphology and stress variations controlled by a relic slab: Implications for deep earthquakes in the Tonga‐Fiji region. Geophys. Res. Lett 48 (7), e2020GL091331. DOI: https://doi.org/10.1029/2020gl091331 Liu, L., Zhou, Q., 2015. Deep recycling of oceanic asthenosphere material during subduction. Geophys. Res. Lett 42 (7), 2204-2211. DOI: https://doi.org/10.1002/2015gl063633 Long, M.D., Silver, P.G., 2009. Mantle flow in subduction systems: The subslab flow field and implications for mantle dynamics. J. Geophys. Res. Solid Earth 114 (B10), B10312. DOI: https://doi.org/10.1029/2008jb006200 Lu, C., Grand, S.P., Lai, H., Garnero, E.J., 2019. TX2019slab: A new P and S tomography model incorporating subducting slabs. J. Geophys. Res. Solid Earth 124 (11), 11549-11567. DOI: https://doi.org/10.1029/2019jb017448 Martinez, F., Taylor, B., 2002. Mantle wedge control on back-arc crustal accretion. Nature 416, 417-420. DOI: https://doi.org/10.1038/416417a Myhill, R., 2013. Slab buckling and its effect on the distributions and focal mechanisms of deep- focus earthquakes. Geophysical Journal International 192 (2), 837-853. DOI: https://doi.org/10.1093/gji/ggs054 Obayashi, M., Fukao, Y., Yoshimitsu, J., 2017. Unusually deep Bonin earthquake of 30 May 172 2015: A precursory signal to slab penetration? Earth Planet. Sci. Lett 459, 221-226. DOI: https://doi.org/10.1016/j.epsl.2016.11.019 Packham, G.H., 1973. A speculative Phanerozoic history of the Southwest Pacific, in: Coleman, P.J. (Ed.), The Western Pacific: Island Arcs, Marginal Seas, Geochemistry. University of Western Australia Press, Nedlands. Peng, D., Stegman, D., 2024. Modeling subduction with extremely fast trench retreat. ESS Open Archive. DOI: https://doi.org/10.22541/essoar.171288715.53824453/v1 Pesicek, J.D., Thurber, C.H., Zhang, H., DeShon, H.R., Engdahl, E.R., Widiyantoro, S., 2010. Teleseismic double-difference relocation of earthquakes along the Sumatra-Andaman subduction zone using a 3-D model. J. Geophys 115 (B10). DOI: https://doi.org/10.1029/2010jb007443 Pockalny, R.A., Druken, K.A., Kincaid, C.R., Kelley, K.A., Graham, D.W., 2009. Implications of plume-trench interaction models for the formation of Niue Trough. AGU EOS Transactions, Abstract No. V31D–2001. Richards, S., Holm, R., Barber, G., 2011. When slabs collide: A tectonic assessment of deep earthquakes in the Tonga-Vanuatu region. Geology 39 (8), 787-790. DOI: https://doi.org/10.1130/g31937.1 Schellart, W.P., Spakman, W., 2012. Mantle constraints on the plate tectonic evolution of the Tonga–Kermadec–Hikurangi subduction zone and the South Fiji Basin region. Australian Journal of Earth Sciences 59 (6), 933-952. DOI: https://doi.org/10.1080/08120099.2012.679692 Smith, G.P., Wiens, D.A., Fischer, K.M., Dorman, L.M., Webb, S.C., Hildebrand, J.A., 2001. A complex pattern of mantle flow in the Lau backarc. Science 292 (5517), 713-716. DOI: https://doi.org/10.1126/science.1058763 Song, T.-R.A., Kawakatsu, H., 2012. Subduction of oceanic asthenosphere: Evidence from sub- slab seismic anisotropy. Geophys. Res. Lett 39 (17), L17301. DOI: https://doi.org/10.1029/2012gl052639 Tang, Y., Obayashi, M., Niu, F., Grand, S.P., Chen, Y.J., Kawakatsu, H., Tanaka, S., Ning, J., Ni, J.F., 2014. Changbaishan volcanism in northeast China linked to subduction-induced mantle upwelling. Nature Geoscience 7 (6), 470-475. DOI: https://doi.org/10.1038/ngeo2166 Tibi, R., Wiens, D.A., 2005. Detailed structure and sharpness of upper mantle discontinuities in the Tonga subduction zone from regional broadband arrays. J. Geophys. Res. Solid Earth 110 (B6). DOI: https://doi.org/10.1029/2004jb003433 van der Hilst, R.D., 1995. Complex morphology of subducted lithosphere in the mantle beneath 173 the Tonga trench. Nature 374, 154-157. DOI: https://doi.org/10.1038/374154a0 Wang, F., Wei, S.S., 2020. Deep low-velocity layer in the Tonga slab, American Geophysical Union. Wang, Z., Zhao, D., Liu, X., Chen, C., Li, X., 2017. P and S wave attenuation tomography of the Japan subduction zone. Geochemistry, Geophys. Geosystems 18 (4), 1688-1710. DOI: https://doi.org/10.1002/2017gc006800 Wei, S.S., Wiens, D.A., 2018. P-wave attenuation structure of the Lau back-arc basin and implications for mantle wedge processes. Earth Planet. Sci. Lett 502, 187-199. DOI: https://doi.org/10.1016/j.epsl.2018.09.005 Wei, S.S., Wiens, D.A., 2020. High bulk and shear attenuation due to partial melt in the Tonga‐ Lau back‐arc mantle. J. Geophys. Res. Solid Earth 125 (1), e2019JB017527. DOI: https://doi.org/10.1029/2019jb017527 Wei, S.S., Wiens, D.A., van Keken, P.E., Cai, C., 2017. Slab temperature controls on the Tonga double seismic zone and slab mantle dehydration. Sci. Adv 3 (1), e1601755. DOI: https://doi.org/10.1126/sciadv.1601755 Wei, S.S., Zha, Y., Shen, W., Wiens, D.A., Conder, J.A., Webb, S.C., 2016. Upper mantle structure of the Tonga‐Lau‐Fiji region from Rayleigh wave tomography. Geochemistry, Geophys. Geosystems 17 (11), 4705-4724. DOI: https://doi.org/10.1002/2016gc006656 Wiens, D.A., McGuire, J.J., 2000. Aftershocks of the March 9, 1994, Tonga earthquake: The strongest known deep aftershock sequence. J. Geophys 105 (B8), 19067-19083. DOI: https://doi.org/10.1029/2000jb900097 Wiens, D.A., Shore, P.J., McGuire, J.J., Roth, E., Bevis, M.G., Draunidalo, K., 1995. The southwest Pacific seismic experiment. IRIS Newsl 14 (1), 1-4. Xi, Z., Chen, M., Wei, S.S., Li, J., Zhou, T., Wang, B., Kim, Y., 2024a. EARA2024: A new radially anisotropic seismic velocity model for the crust and upper mantle beneath East Asia and Northwestern Pacific subduction zones. Geophys. J. Int, ggae302. DOI: https://doi.org/10.1093/gji/ggae302/7740793 Xi, Z., Wei, S.S., Zhu, W., Beroza, G.C., Jie, Y., Saloor, N., 2024b. Deep learning for deep earthquakes: insights from OBS observations of the Tonga subduction zone. Geophys. J. Int 238 (2), 1073-1088. DOI: https://doi.org/10.1093/gji/ggae200 Ye, L., Lay, T., Zhan, Z., Kanamori, H., Hao, J.-L., 2016. The isolated 680 km deep 30 May 2015 MW 7.9 Ogasawara (Bonin) Islands earthquake. Earth Planet. Sci. Lett 433, 169- 179. DOI: https://doi.org/10.1016/j.epsl.2015.10.049 ∼ Zha, Y., Webb, S.C., Wei, S.S., Wiens, D.A., Blackman, D.K., Menke, W., Dunn, R.A., Conder, 174 J.A., 2014. Seismological imaging of ridge–arc interaction beneath the Eastern Lau Spreading Center from OBS ambient noise tomography. Earth Planet. Sci. Lett 408, 194- 206. DOI: https://doi.org/10.1016/j.epsl.2014.10.019 Zhan, Z., 2020. Mechanisms and implications of deep earthquakes. Annu Rev Earth Planet Sci 48 (1), 147-174. DOI: https://doi.org/10.1146/annurev-earth-053018-060314 Zhang, H., Wang, F., Myhill, R., Guo, H., 2019. Slab morphology and deformation beneath Izu- Bonin. Nat Commun 10 (1), 1310. DOI: https://doi.org/10.1038/s41467-019-09279-7 Zhao, D., Xu, Y., Wiens, D.A., Dorman, L.M., Hildebrand, J.A., Webb, S., 1997. Depth extent of the Lau back-arc spreading center and its relation to subduction processes. Science 278 (5336), 254-257. DOI: https://doi.org/10.1126/science.278.5336.254 175