IMPLICATIONS OF FLOW AND THERMOCHEMICAL HETEROGENEITIES IN THE EARTH’S LOWERMOST MANTLE FOR PROPERTIES OF EARTH’S LOWER MANTLE By Jiaxin Zhang A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Geological Sciences – Doctor of Philosophy 2023 ABSTRACT Mantle convection is an important heat and material transport mechanism between the Earth’s deep interior and the surface. The Earth’s deep mantle is seismically heterogeneous and the lower boundary layer of the large-scale mantle dynamics but remains poorly resolved by seismic observations and poorly understood theoretically. Seismic heterogeneities in the Earth’s deep mantle include two Large Low Shear Velocity Provinces (LLSVPs) surrounded by subducted slabs with faster-than-average seismic velocities and numerous patches of ultra-low-velocity- zones (ULVZs) (Chapter 1). The overarching goal of this dissertation is to understand the origin, nature, and evolution of these structures using numerical geodynamic models (Chapter 2). Chapter 3 studies how to link lowermost mantle flow patterns with lower mantle slab deformation styles. Slab deformational styles in the Earth’s lower mantle are often under-resolved in seismic tomography models. In lab experiments, it was found that slab deformation was controlled by slab strength. The strong seismic anisotropy observed in the lowermost mantle is hypothesized to be induced by slab deformation and may be informative in the lowermost mantle flow patterns. To correlate seismic anisotropy with deformation, the first step is to link deformation with flow patterns. I performed 3D spherical calculations and found that slab strength affected lower mantle slab deformation styles and lowermost mantle flow patterns. However, slab trench length has a trade-off effect on lowermost mantle flow patterns, making it difficult to infer lower mantle slab deformation styles solely from lowermost mantle flow patterns. Chapter 4 investigates if the symmetry of ULVZ cross-sectional shapes implies the viscosity of LLSVPs with respect to the background mantle. The shape of ULVZs is controlled by their properties (e.g., density, viscosity) and their interaction with the LLSVPs and the surrounding mantle. Therefore, ULVZ shapes provide information about the rheological properties of the LLSVPs. In geodynamic models, ULVZs and LLSVPs material usually have the same composition-dependent viscosity as the background mantle. However, this might be unreasonable if they originate from compositionally distinctive materials. I incorporated the composition- dependence of viscosity into calculations and found that the composition-dependence of viscosity affects the viscous coupling among LLSVP materials, ULVZ materials, and the background mantle. The viscous coupling, therefore, controls the symmetry of ULVZ shapes. In Chapter 5, I explored how dense thermochemical piles with increased pile intrinsic viscosity respond to changing convective flow patterns in terms of lateral mobility and the stability of their cross-sectional morphologies. Whether LLSVPs are fixed at their current position for up to a few hundred million years is under debate. In geodynamic models, LLSVPs, if caused by thermochemical piles, appeared as passive structures swept around by subducted slabs. However, other parameters, such as compositional viscosity contrast between piles and the background mantle (as explored in this study), have not been sufficiently explored in these models. Using the method developed in Chapter 4, I perturb upwelling and downwelling convective patterns in models with increased pile intrinsic viscosity and examine the long-term stability of thermochemical piles. I found that piles are mobile even if their compositional viscosity increase by 5,000x, but their cross-sectional morphology are more stable as pile viscosity increases. Chapter 6 concludes these research projects and discusses perspectives for future work. This thesis is dedicated to my parents and fiancé. Thank you for always being supportive and faithful to me. iv ACKNOWLEDGEMENTS Firstly and mostly, I want to thank my Ph.D. advisor, Allen McNamara, for supporting, inspiring, and believing in me. My first interest in Geodynamics was inspired by one of your papers. Later on, during my Ph.D., A lot of our small talks on earth science have inspired me. Thank you for challenging and guiding me to develop scientific ideas, become a better scientific writer and presenter, independently learn and develop numerical tools, and, most importantly, obtain a big-picture view of the research questions to ask. You have always firmly believed in my abilities and that I can be a geodynamicist, even when I was confused about my career path. I want to thank the rest of my Ph.D. committee: Min Chen, Susannah Dorfman, Kevin Mackey, and Songqiao (Shawn) Wei, for all the support and important advice you have given to my research and future career. I obtained much knowledge and perspectives from you outside of my field. Min, I wish you were on my defense, as your passion for science and helping other people have always inspired me. Your faith in and caring for me helped me through some challenging moments during my Ph.D. I want to thank my friends Ying Wang, Yike Shen, Luwen Wan, Yunfei Shao, Ze Zhang, and Zheng Li. The joy and support you have brought to my life are meaningful to me. I feel fortunate always to have you around when needed, even if I have not met with some of you in person for years. I want to thank my geophysics colleagues at Michigan State University for our science conversations and your help in school life. The lab is a welcoming and inclusive environment! I want to give special thanks to Heidi Krauss for the energy you brought to the lab that ignited me. I want to thank my dogs, Ollie and Molly, for your unconditional love for me. Raising you is like a raising experience for my internal self. v I want to thank my parents for supporting me, having my back, and putting me in a math and science context as I grew up. I still remember the endless times you were solving a math problem or physics problem with me at midnight, even though some of them were naive. You gave me the foundation to grow up trying to be a scientist. Lastly, I want to thank my husband, Guanyuan Shuai, for supporting me. You are always there when I need you. You did not hesitate to support me when I had an opportunity to do the research I am interested in within a country different from where you work, even if it is at least a 20-hour flight away and we just got married. Life is much sweeter with you. vi TABLE OF CONTENTS LIST OF ABBREVIATIONS ........................................................................................................ ix CHAPTER 1: INTRODUCTION ............................................................................................. 1 1.1 Background and Rationale ............................................................................................... 1 1.2 Objectives of the Dissertation .......................................................................................... 5 REFERENCES ......................................................................................................................... 11 CHAPTER 2: GOVERNING EQUATIONS AND NUMERICAL METHODS ................... 16 2.1 Governing Equations ...................................................................................................... 17 2.2 Numerical Methods ........................................................................................................ 24 REFERENCES ......................................................................................................................... 26 LOWERMOST MANTLE FLOW PATTERNS CAUSED BY SLABS OF CHAPTER 3: DIFFERENT STRENGTH AND TRENCH LENGTH IMPINGING ON THE CORE-MANTLE BOUNDARY…………………………………………………………………………………….27 3.1 Abstract .......................................................................................................................... 27 Introduction .................................................................................................................... 28 3.2 3.3 Methods .......................................................................................................................... 36 3.4 Results ............................................................................................................................ 43 Lowermost mantle flow patterns induced by slabs as a function of slab strength .. 43 Lowermost mantle flow patterns induced by slabs as a function of trench length . 48 3.5 Discussions ..................................................................................................................... 51 3.6 Conclusions .................................................................................................................... 54 REFERENCES ......................................................................................................................... 55 APPENDIX ............................................................................................................................... 62 3.4.1 3.4.2 CHAPTER 4: ULVZ CROSS-SECTIONAL SHAPE GIVES CLUE TO THE ORIGIN AND INTRINSIC RHEOLOGICAL NATURE OF THEMSELVES AND LLSVPS .......................... 67 4.1 Abstract .......................................................................................................................... 67 4.2 Introduction .................................................................................................................... 68 4.3 Method ........................................................................................................................... 71 4.4 Results ............................................................................................................................ 72 4.5 Discussion and Conclusion ............................................................................................ 78 REFERENCES ......................................................................................................................... 81 APPENDIX ............................................................................................................................... 85 ARE INTRINSICALLY MORE VISCOUS THERMOCHEMICAL PILES CHAPTER 5: STATIONARY AT THE CORE-MANTLE BOUNDARY? ....................................................... 99 5.1 Abstract .......................................................................................................................... 99 5.2 Introduction .................................................................................................................. 100 5.3 The Experiment Setup .................................................................................................. 105 5.4 Methods ........................................................................................................................ 108 5.5 Results .......................................................................................................................... 111 5.6 Discussion and Conclusion .......................................................................................... 126 vii REFERENCES ....................................................................................................................... 133 APPENDIX ............................................................................................................................. 143 CHAPTER 6: CONCLUSIONS AND PERSPECTIVES .......................................................... 148 REFERENCES ....................................................................................................................... 153 viii LIST OF ABBREVIATIONS 2D Two-Dimensional 3D Three-Dimensional Br Bridgmanite CMB Core-Mantle Boundary CPO Crystal Preferred Orientation Fp Ferropericlase LIPs Large Igneous Provinces LLSVPs Large Low Shear Velocity Provinces MORBs Mid-Ocean Ridge Basalts OIBs Oceanic Island Basalts PGZ Plume Generation Zones pPV Post-Perovskite SPO Shape Preferred Orientation ULVZ Ultra-Low Velocity Zones ix CHAPTER 1: INTRODUCTION 1.1 Background and Rationale The dynamical processes within the Earth’s interior determine Earth’s surface evolution and impact surface hazardous geological activities. Near the Earth’s surface, plate tectonics is the process that generates many of the Earth’s surface geological phenomena, such as mountain building, motion of tectonic plates, nearly all volcanism, and most earthquakes. Furthermore, it is critical to maintain the Earth’s geomagnetic field (e.g., Olson, 2016), which shelters the Earth from erosion and solar radiation (e.g., Tarduno et al., 2010). This protectant ensures a secure environment for the Earth to evolve as a habitable planet and is generated by another dynamic process, geodynamo, in the Earth’s core. The Earth’s mantle convection drives plate tectonics and controls the geodynamo through core heat loss (e.g., Biggin et al., 2012). Mantle convection is one of the essential mechanisms of heat and material transport between the Earth’s deep interior and the surface. Understanding the Earth’s mantle convection gives an insight into what processes are necessary to develop and maintain the habitability of the Earth and how the entire planet functions (e.g., Jellinek & Jackson, 2015). The lowermost few hundred kilometers of the Earth’s mantle, as the bottom boundary layer of mantle convection, is not well-constrained due to the lack of direct observations but is one of the most complicated regions on the Earth. Thus, understanding properties, structures, and dynamic processes within the Earth’s deep mantle often requires a combination of advances in geochemistry, geodesy, paleomagnetics, mineral physics, seismology, and geodynamics. To the first order, global seismic tomography models reveal two Large Low Shear Velocity Provinces (LLSVPs) with a height of a few hundred kilometers to up to a thousand kilometers, lying beneath the Pacific and Africa at the base of the mantle (e.g., Dziewonski et al., 2010; French 1 & Romanowicz, 2014; Grand, 2002; Hosseini et al., 2020; Houser et al., 2008; Lekic et al., 2012; Megnin & Romanowicz, 2000; Ritsema et al., 2011; Schuberth et al., 2009; Simmons et al., 2010). Figure 1.1 Seismic tomography model S20RTS. Isosurfaces representing -0.6% shear-wave velocity anomaly that highlights the presence of large low-shear-velocity provinces (LLSVPs) beneath the Pacific and Africa. The background is colored by shear-wave velocity variation, with a color scale at the bottom. Taken from Bull et al., 2009. Figure 1.1 illustrates a 3D seismic tomography model of the LLSVPs on a map view of the model at the lowermost mantle depth. All models agree the seismic velocity of LLSVPs is a few percent slower than the average seismic velocity at the depth. In tomography models, the LLSVPs are surrounded by higher-than-average seismic velocity regions, which are hypothesized to be subducted slabs or their remnants (e.g., Dziewonski et al., 2010; Grand, 2002; Ritsema et al., 2011; Romanowicz, 2003). At a scale of at least one order magnitude smaller than that of the LLSVPs, numerous small thin patches of ultra-low velocity zones (ULVZs), with a thicknesses ranging from a few to 100 km and lateral sizes ranging from tens to around a thousand km, are detected by regional waveform studies near the core-mantle boundary (CMB) (e.g., Yu & Garnero, 2018). These ULVZs have -10% to -30% seismic velocity perturbations 2 and are located near, within, and sometimes outside the LLSVPs (Figure 1.2). In the next section (Chapter 1.2), I introduce the hypothesis for the origin of the LLSVPs. Figure 1.2 Map view of results for ULVZ detecting investigations near the CMB, showing ULVZ detection, none ULVZ detection, and ULVZ detection with uncertainties. Pink regions are LLSVP regions. Taken from Yu and Garnero, 2018. and ULVZs, as well as other related observations in more detail. In some of the hypotheses, the LLSVPs and ULVZs are hypothesized to be caused by compositional distinctive materials. 3 These heterogeneities add complexity to the large-scale mantle convection patterns and have an impact on the CMB heat flux (e.g., Li et al., 2018; Nakagawa et al., 2010), which is crucial to our understanding of the heat budget and thermochemical evolution of the Earth, and the nature of mantle convection. With the advance of lab technology and rapid growth of seismic data, there is still fair uncertainty to first-order features (LLSVPs) in the deep mantle and key properties that govern the dynamics in the lower mantle. For instance, the size and volume of the LLSVPs are not in agreement among papers. The percentage of LLSVP area over CMB area ranges from around 20% (e.g., Burke et al., 2008) to 50% (e.g., Garnero & McNamara, 2008). The volume of LLSVPs is estimated to be 1.6% (e.g., Burke et al., 2008) and up to 8% (e.g., Cottaar & Lekic, 2016). In addition, the density excess of LLSVPs is not in agreement among literatures. Some argue that LLSVPs have a small amount of density excess (e.g., Ishii & Tromp, 1999; Lau et al., 2017) while others argue that they are slightly less dense than the background mantle (e.g., Koelemeijer et al., 2017). Both size and density excess, however, are crucial parameters that input to geodynamic models and determine whether and how long LLSVPs can be stable along the CMB. Lastly, the viscosity of the lower mantle is important to understand mantle convection but as poorly constrained as to have 3-4 orders of magnitude uncertainty. The overarching goal of this dissertation is to better understand the fate and properties, especially rheological properties, of the heterogeneous seismic structures in the Earth’s deep mantle (LLSVPs, ULVZs, and slabs). I perform numerical geodynamical experiments under assumptions (details described in Chapter 2) for mantle convection. The governing equations for mantle convection and the numerical methods used to solve these equations are illustrated in Chapter 2. I examine how the deformation, evolution, and geometries of the seismically 4 heterogeneous structures in the Earth’s deep mantle are controlled by their rheological and physical properties, such as viscosity and density. The hypotheses tested in this dissertation are to not only provide a better fundamental understanding of the dynamical processes associated with the heterogeneities but also bridge geodynamic models and observations, or future geophysical observations if they are not yet available, so that we may eventually infer the properties and dynamical processes in the Earth’s deep mantle from observations. In the next section (Chapter 1.2), I introduce the more specific objectives of Chapters 3-5. 1.2 Objectives of the Dissertation Chapters 3-5 are inspired by the following important questions for understanding the origin, nature, and evolution of the lowermost mantle heterogeneities: What is the fate of the slabs sinking to the deep mantle? What causes the LLSVPs? What is the viscosity of the LLSVPs? Can we use the ULVZs to infer the properties of the LLSVPs? How will the LLSVPs evolve? Are the LLSVPs stable? Is there a mechanism to explain the possible stability of the LLSVPs? The equations and numerical methods applied to the geodynamic studies in Chapters 3-5 are described in Chapter 2. In Chapter 3, I investigated if lowermost mantle flow patterns can reflect slab trench length and the viscosity contrast between slabs and the background mantle (slab strength). Subducted slabs are one of the driving forces for mantle convection and plate tectonics. Lower mantle slab dynamics are critical to studying how materials are transported and recycled between the surface and the CMB. However, in seismic tomography models, slabs in the Earth’s lower mantle are often under-resolved, hindering the understanding of how they deform and what they look like in the lower mantle. Seismic anisotropy is informative in the mantle flow patterns and the deformation in the region. The lowermost few hundred kilometers of the Earth’s mantle exhibit strong seismic 5 anisotropy, hypothesized to be caused by lattice-preferred orientation induced by the deformation of slabs impinging on the CMB or shape-preferred orientation. In this Chapter, I aimed to ultimately bridge lower mantle slab deformation and lowermost mantle seismic anisotropy observations, so I focused on the former hypothesis. I conducted 3D spherical convection calculations using CitcomCU (e.g., Zhong, 2006) and modified surface boundary conditions to induce subducting slabs. I found that weaker slabs and slabs with shorter trench lengths were deformed under pure shear, inducing azimuthal-divergent lowermost mantle flow patterns. On the contrary, stronger slabs and slabs with longer trench lengths underwent buckling or bending, inducing bilateral lowermost mantle flow patterns. Combined with the lowermost mantle seismic anisotropy observations, this work is valuable in examining the dynamical behavior and morphology of slabs in the lower mantle. Moreover, it provides clues about the properties of slabs that reach the CMB. The origin and evolution of LLSVPs are under debate. Several conceptual models of large- scale mantle structure and convection have been generated based on these seismic observations (Figure 1.3). There are hypothesized models that argue LLSVPs are caused by thermal structures, which are either large megaplumes or clusters of smaller-scale plumes (Figure 1.3a, e.g., Thompson & Tackley, 1999; Ritsema et al., 2007; Davies et al., 2012; Schubert et al., 2004). However, how these isochemical models can satisfy seismic evidence of compositionally heterogeneous LLSVPs is unclear (e.g., McNamara, 2019). Other conceptual models hypothesize that the mantle is undergoing large-scale thermochemical convection and LLSVPs are caused by intrinsically more dense compositional reservoirs within the lowermost mantle. If the reservoir is intrinsically denser than the surrounding mantle within a few percent, the reservoirs will form thermochemical domes or superplumes. Superplumes rise and sink periodically as they cool and 6 heat (Figure 1.3b, c, e.g., Davaille & Vatteville, 2005; Tackley, 2012). If, on the other hand, the intrinsic density of the reservoir is even greater, hot and dense material at the core-mantle boundary (CMB) is passively swept or advected to the upwelling center by cold and more viscous subducted slabs and forms vertically stable thermochemical piles with ridge-like shape and sharp boundaries (Figure 1.3d, e.g., Tackley, 1998; McNamara and Zhong, 2005). In this model, the thermochemical piles are caused by the long-term accumulations of subducted oceanic crust (e.g., Brandenburg et al., 2008; Brandenburg & van Keken, 2007; Christensen & Hofmann, 1994; G. F. Davies, 2008; Nakagawa & Tackley, 2006, 2011) or a combination of early differentiation processes and early accumulations of subducted oceanic crust (e.g., Bull et al., 2009; Jellinek & Manga, 2004; McNamara & Zhong, 2004; Tackley, 1998, 2000; Zhang et al., 2010). The hypothesis that the LLSVPs are “thermochemical piles” is more satisfactory to seismic observations towards the LLSVPs (i.e., sharp boundaries, tomography models, and compositional heterogeneity, details discussed in Chapters 4 and 5) and are widely accepted. ULVZs are hypothesized to be caused by partial melting, ultra-dense materials distinct from the LLSVP compositions, or a combination of both. If LLSVPs and ULVZs are hypothesized to be compositionally distinct from the background mantle, they not only have a different composition but should also have a different mineralogical grain size. Both consequences affect the intrinsic viscosity contrast between LLSVPs and the background mantle, leading to different rheological formulations for piles and the ambient mantle. In particular, the LLSVPs may have an intrinsically higher viscosity than the ambient mantle. Because it is currently impossible to directly constrain the grain size/rheologies of the lower mantle from real Earth samples, it is important to study how the increased intrinsic viscosity of thermochemical piles would affect the dynamical behavior of piles (LLSVPs) and 7 ultra-dense patches (ULVZs). However, composition-dependent viscosity is underestimated and usually not included in conventional thermochemical convection models. Using the ratio tracer method (e.g., Tackley & King, 2003), I incorporate a multi-composition thermochemical extension into CitcomCU. I applied composition-dependent rheology formulation in this extension and performed calculations with this extension code for the studies in Chapters 4 and 5. The increased pile intrinsic viscosity might affect the shape of ULVZs. The shape of ULVZs is controlled by ULVZ properties (e.g., density, viscosity) and their interaction with the LLSVPs and the surrounding mantle (e.g., M. Li et al., 2017). Therefore, ULVZ shapes provide information about the lowermost mantle flow and rheological properties of the LLSVPs. In Chapter 4, I explored if the symmetry of ULVZ cross-sectional shapes implies the viscosity of LLSVPs with respect to the ambient mantle. I found that if piles had the same or less intrinsic viscosity as the ambient mantle, ULVZ accumulations at pile edges were thicker on the outboard side of piles, leading to an asymmetric ULVZ shape. However, if compositional reservoirs had higher intrinsic viscosity than the surrounding mantle, ULVZs at pile edges typically became more symmetric or exhibited the opposite asymmetry, thinner on the outboard side of piles. In addition, I found that ULVZs within the pile interior are consistently symmetrical. The finding is important, as it argues that even if ULVZs originate from ultra-dense materials, they can accumulate into symmetrical patches at LLSVP edges, just as if they are caused by partial melting (e.g., M. Li et al., 2017). Furthermore, high-resolution seismic images of ULVZs may be available in the future to advance the understanding of the viscosity, nature, and origin of LLSVPs. 8 Figure 1.3 Conceptual models of the large-scale convection that can result in an LLSVP. (a) In an isochemical mantle, upwellings and plumes can be concentrated away beneath downwellings, and upwellings can be an LLSVP. (b) Former oceanic crust can introduce small-scale heterogeneities into the lowermost mantle. (c) A nearly neutrally buoyant thermochemical structure can have vertical sides. If this structure is unstable (meta-stable pile), it may break away from the CMB as a superplume. (d) with larger density, the thermochemical pile becomes a stable structure with sloped sides and plumes coming off bridges. From Zhao et al., 2015. The increased pile intrinsic viscosity might also influence the stability of the LLSVPs. Such stability usually contains two aspects: lateral mobility and stability of the cross-sectional morphology of the LLSVPs. In Chapter 5, I investigated how dense thermochemical piles with increased pile intrinsic viscosity respond to changing convective flow patterns in terms of lateral mobility and morphological stability. Some paleomagnetic studies (e.g., Steinberger & Torsvik, 2008, 2012) hypothesized that LLSVPs could be fixed at their current position for up to a few hundred million years. On the contrary, a recent study (e.g., Flament et al., 2022) indicated that mobile LLSVPs could non-uniquely explain the paleomagnetic data. In geodynamic models, LLSVPs appeared as passive structures swept around by subducted slabs. However, other parameters, such as compositional viscosity contrast between piles and the ambient mantle, have not been sufficiently explored in these models. Therefore, it’s necessary and significant to explore if composition-dependent rheology could cause the stability of piles against changing mantle flow 9 patterns. I found that the increased intrinsic viscosity of thermochemical piles does not cause them to be more resistant to changing downwelling patterns. However, piles with higher intrinsic viscosity are more resistant to changes in their cross-sectional morphology in response to changing upwelling flows in the background mantle. This study shows that if LLSVPs are indeed fixed, another mechanism must be found to explain the lateral immobility of piles. Chapter 6 summarizes the findings in Chapters 3-5 and recommends some future studies related to this dissertation. 10 REFERENCES Biggin, A. J., Steinberger, B., Aubert, J., Suttie, N., Holme, R., Torsvik, T. H., van der Meer, D. G., & van Hinsbergen, D. J. J. (2012). Possible links between long-term geomagnetic variations and whole-mantle convection processes. Nature Geoscience, 5(8), 526–533. https://doi.org/10.1038/ngeo1521 Brandenburg, J. P., Hauri, E. H., van Keken, P. E., & Ballentine, C. J. (2008). A multiple-system study of the geochemical evolution of the mantle with force-balanced plates and thermochemical effects. Earth and Planetary Science Letters, 276(1–2), 1–13. https://doi.org/10.1016/J.EPSL.2008.08.027 Brandenburg, J. P., & van Keken, P. E. (2007). Deep storage of oceanic crust in a vigorously convecting mantle. Journal of Geophysical Research: Solid Earth, 112(6). https://doi.org/10.1029/2006JB004813 Bull, A. L., McNamara, A. K., & Ritsema, J. (2009). Synthetic tomography of plume clusters and thermochemical piles. Earth and Planetary Science Letters, 278(3–4), 152–162. https://doi.org/10.1016/J.EPSL.2008.11.018 Burke, K., Steinberger, B., Torsvik, T. H., & Smethurst, M. A. (2008). Plume Generation Zones at the margins of Large Low Shear Velocity Provinces on the core–mantle boundary. Earth and Planetary Science Letters, 265(1–2), 49–60. https://doi.org/10.1016/J.EPSL.2007.09.042 Christensen, U. R., & Hofmann, A. W. (1994). Segregation of subducted oceanic crust in the convecting mantle. In JOURNAL OF GEOPHYSICAL RESEARCH (Vol. 99, Issue B10). https://doi.org/10.1029/93JB03403 Cottaar, S., & Lekic, V. (2016). Morphology of seismically slow lower-mantle structures. Geophysical Journal International, 207(2), 1122–1136. https://doi.org/10.1093/gji/ggw324 Davaille, A., & Vatteville, J. (2005). On the transient nature of mantle plumes. Geophysical Research Letters, 32(14), n/a-n/a. https://doi.org/10.1029/2005GL023029 Davies, D. R., Goes, S., Davies, J. H., Schuberth, B. S. A., Bunge, H. P., & Ritsema, J. (2012). Reconciling dynamic and seismic models of Earth’s lower mantle: The dominant role of thermal heterogeneity. Earth and Planetary Science Letters, 353–354, 253–269. https://doi.org/10.1016/J.EPSL.2012.08.016 Davies, G. F. (2008). Episodic layering of the early mantle by the ‘basalt barrier’ mechanism. Earth and Planetary Science Letters, 275(3–4), 382–392. https://doi.org/10.1016/j.epsl.2008.08.036 Dziewonski, A. M., Lekic, V., & Romanowicz, B. A. (2010). Mantle Anchor Structure: An argument for bottom up tectonics. Earth and Planetary Science Letters, 299(1–2), 69–79. https://doi.org/10.1016/j.epsl.2010.08.013 11 Flament, N., Bodur, Ö. F., Williams, S. E., & Merdith, A. S. (2022). Assembly of the basal mantle structure beneath Africa. 603(March). https://doi.org/10.1038/s41586-022-04538-y French, S. W., & Romanowicz, B. A. (2014). Whole-mantle radially anisotropic shear velocity structure from spectral-element waveform tomography. Geophysical Journal International, 199(3), 1303–1327. https://doi.org/10.1093/gji/ggu334 Garnero, E. J., & McNamara, A. K. (2008). Structure and Dynamics of Earth’s Lower Mantle. Science, 320(5876), 626–628. https://doi.org/10.1126/science.1148028 Grand, S. P. (2002). Mantle shear–wave tomography and the fate of subducted slabs. Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 360(1800), 2475–2491. https://doi.org/10.1098/rsta.2002.1077 Hosseini, K., Sigloch, K., Tsekhmistrenko, M., Zaheri, A., Nissen-Meyer, T., & Igel, H. (2020). Global mantle structure from multifrequency tomography using P, PP and P-diffracted waves. Geophysical Journal International, 220(1), 96–141. https://doi.org/10.1093/gji/ggz394 Houser, C., Masters, G., Shearer, P., & Laske, G. (2008). Shear and compressional velocity models of the mantle from cluster analysis of long-period waveforms. Geophysical Journal International, 174(1), 195–212. https://doi.org/10.1111/J.1365-246X.2008.03763.X/3/174- 1-195-FIG016.JPEG Ishii, M., & Tromp, J. (1999). Normal-mode and free-air gravity constraints on lateral variations in velocity and density of earth’s mantle. Science, 285(5431), 1231–1236. https://doi.org/10.1126/SCIENCE.285.5431.1231/SUPPL_FILE/1040466S9_THUMB.GIF Jellinek, A. M., & Jackson, M. G. (2015). Connections between the bulk composition, geodynamics and habitability of Earth. Nature Geoscience, 8(8), 587–593. https://doi.org/10.1038/ngeo2488 Jellinek, A. M., & Manga, M. (2004). Links between long-lived hot spots, mantle plumes, D″, and plate tectonics. Reviews of Geophysics, 42(3). https://doi.org/10.1029/2003RG000144 Koelemeijer, P., Deuss, A., & Ritsema, J. (2017). Density structure of Earth’s lowermost mantle from Stoneley mode splitting observations. Nature Communications, 8. Lau, H. C. P., Mitrovica, J. X., Davis, J. L., Tromp, J., Yang, H. Y., & Al-Attar, D. (2017). Tidal tomography constrains Earth’s deep-mantle buoyancy. Nature, 551(7680), 321–326. https://doi.org/10.1038/nature24452 Lekic, V., Cottaar, S., Dziewonski, A., & Romanowicz, B. (2012). Cluster analysis of global lower mantle tomography: A new class of structure and implications for chemical heterogeneity. Earth and Planetary Science Letters, 357–358, 68–77. https://doi.org/10.1016/J.EPSL.2012.09.014 12 Li, M., McNamara, A. K., Garnero, E. J., & Yu, S. (2017). Compositionally-distinct ultra-low velocity zones on Earth’s core-mantle boundary. Nature Communications, 8(1), 1–9. https://doi.org/10.1038/s41467-017-00219-x Li, M., Zhong, S., & Olson, P. (2018). Linking lowermost mantle structure, core-mantle boundary heat flux and mantle plume formation. Physics of the Earth and Planetary Interiors, 277, 10–29. https://doi.org/10.1016/j.pepi.2018.01.010 McNamara, A. K. (2019). A review of large low shear velocity provinces and ultra low velocity zones. Tectonophysics, 760, 199–220. https://doi.org/10.1016/j.tecto.2018.04.015 McNamara, A. K., & Zhong, S. (2004). Thermochemical structures within a spherical mantle: Superplumes or piles? Journal of Geophysical Research: Solid Earth, 109. https://doi.org/10.1029/2003JB002847 McNamara, A. K., & Zhong, S. (2005). Thermochemical structures beneath Africa and the Pacific Ocean. Nature, 437(7062), 1136–1139. https://doi.org/10.1038/NATURE04066 Megnin, C., & Romanowicz, B. (2000). The three-dimensional shear velocity structure of the mantle from the inversion of body, surface and higher-mode waveforms. Geophysical Journal International, 143(3), 709–728. https://doi.org/10.1046/J.1365- 246X.2000.00298.X/2/M_143-3-709-IEQ052.JPEG Nakagawa, T., & Tackley, P. J. (2006). Three-dimensional structures and dynamics in the deep mantle: Effects of post-perovskite phase change and deep mantle layering. Geophysical Research Letters, 33(12). https://doi.org/10.1029/2006GL025719 Nakagawa, T., & Tackley, P. J. (2011). Effects of low-viscosity post-perovskite on thermo- chemical mantle convection in a 3-D spherical shell. Geophysical Research Letters, 38(4). https://doi.org/10.1029/2010GL046494 Nakagawa, T., Tackley, P. J., Deschamps, F., & Connolly, J. A. D. (2010). The influence of MORB and harzburgite composition on thermo-chemical mantle convection in a 3-D spherical shell with self-consistently calculated mineral physics. Earth and Planetary Science Letters, 296(3–4), 403–412. https://doi.org/10.1016/j.epsl.2010.05.026 Olson, P. (2016). Mantle control of the geodynamo: Consequences of top‐down regulation. Geochemistry, Geophysics, Geosystems, 17(5), 1935–1956. https://doi.org/10.1002/2016GC006334 Ritsema, J., Deuss, A., Van Heijst, H. J., & Woodhouse, J. H. (2011). S40RTS: A degree-40 shear-velocity model for the mantle from new Rayleigh wave dispersion, teleseismic traveltime and normal-mode splitting function measurements. Geophysical Journal International, 184(3), 1223–1236. https://doi.org/10.1111/j.1365-246X.2010.04884.x Ritsema, J., McNamara, A. K., & Bull, A. L. (2007). Tomographic filtering of geodynamic models: Implications for model interpretation and large-scale mantle structure. Journal of Geophysical Research, 112(B1), B01303. https://doi.org/10.1029/2006JB004566 13 Romanowicz, B. (2003). Global mantle tomography: Progress status in the past 10 years. Annual Review of Earth and Planetary Sciences, 31, 303–328. https://doi.org/10.1146/annurev.earth.31.091602.113555 Schubert, G., Masters, G., Olson, P., & Tackley, P. (2004). Superplumes or plume clusters? Physics of the Earth and Planetary Interiors, 146(1–2), 147–162. https://doi.org/10.1016/j.pepi.2003.09.025 Schuberth, B. S. A., Bunge, H.-P., & Ritsema, J. (2009). Tomographic filtering of high- resolution mantle circulation models: Can seismic heterogeneity be explained by temperature alone? Geochemistry, Geophysics, Geosystems, 10(5), n/a-n/a. https://doi.org/10.1029/2009GC002401 Simmons, N. A., Forte, A. M., Boschi, L., & Grand, S. P. (2010). GyPSuM: A joint tomographic model of mantle density and seismic wave speeds. Journal of Geophysical Research: Solid Earth, 115, 12310. https://doi.org/10.1029/2010JB007631 Steinberger, B., & Torsvik, T. H. (2008). Absolute plate motions and true polar wander in the absence of hotspot tracks. Nature, 452(7187), 620–623. https://doi.org/10.1038/NATURE06824 Steinberger, B., & Torsvik, T. H. (2012). A geodynamic model of plumes from the margins of Large Low Shear Velocity Provinces. Geochemistry, Geophysics, Geosystems, 13(1), n/a- n/a. https://doi.org/10.1029/2011GC003808 Tackley, P. J. (1998). Self-consistent generation of tectonic plates in three-dimensional mantle convection. Earth and Planetary Science Letters, 157(1–2), 9–22. https://doi.org/10.1016/S0012-821X(98)00029-6 Tackley, P. J. (2000). Mantle convection and plate tectonics: Toward an integrated physical and chemical theory. In Science (Vol. 288, Issue 5473, pp. 2002–2007). American Association for the Advancement of Science. https://doi.org/10.1126/science.288.5473.2002 Tackley, P. J. (2012). Dynamics and evolution of the deep mantle resulting from thermal, chemical, phase and melting effects. Earth-Science Reviews, 110(1–4), 1–25. https://doi.org/10.1016/J.EARSCIREV.2011.10.001 Tackley, P. J., & King, S. D. (2003). Testing the tracer ratio method for modeling active compositional fields in mantle convection simulations. Geochemistry, Geophysics, Geosystems, 4(4). https://doi.org/10.1029/2001GC000214 Tarduno, J. A., Cottrell, R. D., Watkeys, M. K., Hofmann, A., Doubrovine, P. V., Mamajek, E. E., Liu, D., Sibeck, D. G., Neukirch, L. P., & Usui, Y. (2010). Geodynamo, Solar Wind, and Magnetopause 3.4 to 3.45 Billion Years Ago. Science, 327(5970), 1238–1240. https://doi.org/10.1126/science.1183445 Thompson, P. E., & Tackley, P. J. (1999). Generation of mega‐plumes from the core‐mantle boundary in a compressible mantle with temperature‐dependent 14 viscosity. In Wiley Online Library (Vol. 25, Issue 11). American Geophysical Union. https://doi.org/10.1029/98GL01228 Yu, S., & Garnero, E. J. (2018). Ultralow Velocity Zone Locations: A Global Assessment. Wiley Online Library, 19(2), 396–414. https://doi.org/10.1002/2017GC007281 Zhang, N., Zhong, S., Leng, W., & Li, Z. X. (2010). A model for the evolution of the Earth’s mantle structure since the Early Paleozoic. Journal of Geophysical Research: Solid Earth, 115(6). https://doi.org/10.1029/2009JB006896 Zhong, S. (2006). Constraints on thermochemical convection of the mantle from plume heat flux, plume excess temperature, and upper mantle temperature. Journal of Geophysical Research: Solid Earth, 111(4), 4409. https://doi.org/10.1029/2005JB003972 15 CHAPTER 2: GOVERNING EQUATIONS AND NUMERICAL METHODS When performing mantle geodynamic calculations, we use numerical tools to solve the equations of conservation of mass, momentum, and energy that govern fluid dynamics. The equations are solved under different approximations for different research scales and purposes, such as the infinite Prandtl number approximation and the Boussinesq approximation. The infinite Prandtl number approximation is usually applied in a buoyancy-driven convection system. The Boussinesq approximation is widely used to simplify incompressible natural convection problems, such as ocean convection, atmosphere convection and mantle convection. It assumes that density is temporally constant, and density variations only affect buoyancy forces in the flow field. It also assumes that there is no viscous dissipation. This approximation is reasonable when studying slabs and plumes traveling through the mantle. The Boussinesq approximation marginally reduces computation costs. In this dissertation, all the calculations are performed under the Boussinesq approximation. But we should note that the Boussinesq approximation is not always valid for mantle dynamic problems. For instance, near boundary layers or phase change locations, where the density gradient is large, volume change due to phase changes can be as large as more than 20%, which cannot be neglected. If the research interest is on, for example, how dehydration phase changes affect the subduction of slabs, other approximations should be applied to (e.g., Gassmöller et al., 2020). The discussion for such problems and approximations is beyond the scope of this dissertation, though. When modeling mantle convection, we apply estimated parameters of mantle material properties to the best of our knowledge. The properties of the Earth’s mantle, especially the lower mantle, however, as described in Chapter 1, have a wide range of uncertainty due to the lack of direct samples from the lowermost mantle. Therefore, it’s important to explore a wide range of 16 parameter space and testify proposed dynamical hypotheses to Better understand the structure and dynamics of the Earth’s mantle. Some proposed dynamical hypotheses directly relate to current seismic observations or mineral physics findings. These hypotheses are motivated to explain observations (i.e., research in Chapter 5). Others are in the “if… then…” style, exploring what observations we might expect if the mantle has certain properties or is in a particular type of dynamical regime. Usually, current observations are not sufficient to test these hypotheses. Sometimes, we do not have a fundamental understanding of these hypotheses at all. Dynamical calculations to test these hypotheses are similar to computational physical experiments. They explore the consequences of the underexplored parameters or dynamical regime, especially those that might be detectable in future seismic studies (i.e., research in Chapters 3 and 4). In this chapter, I first introduce the governing conservation equations and some of the critical parameters I vary in this dissertation (Chapter 2.1). Next, I described the numerical methods and model setups applied to the studies in this dissertation (Chapter 2.2). 2.1 Governing Equations Without approximations, the equation for the conservation of mass is: ∂ρ ∂t + ∇ ∙ (𝜌𝑢⃑ ) = 0……………………………………………………………………Equation (2.1) where 𝜌 is density, t is time, and 𝑢⃑ is velocity. The equation for the conservation of momentum is: 𝜌 𝐷𝑢⃑⃑ 𝐷𝑡 = −∇p + ∇ ∙ 𝜏̿ + 𝜌 𝑔 ………………………………………………………….Equation (2.2) where 𝐷 𝐷𝑡 is the material derivative, p is pressure, 𝜏̿ is deviatoric stress, and 𝑔 is gravitational acceleration. The equation for the conservation of energy is: 17 𝜌𝐶𝑝 𝐷𝑇 𝐷𝑡 = ∇ ∙ 𝑘∇T + α𝑇 𝐷𝑝 𝐷𝑡 + 𝜏̿: ∇𝑢⃑ + 𝜌 𝐻 − 𝑇 𝐷𝑆∗ 𝐷𝑡 ………………………………….Equation (2.3) where 𝐶𝑝 is heat capacity, k is thermal conductivity, α is thermal expansivity, T is temperature, H is heat production rate, and 𝐷𝑆∗ is entropy changes due to additional processes. In this equation, ∇ ∙ 𝑘∇T describes energy produced from heat conduction, α𝑇 𝐷𝑃 𝐷𝑡 represents energy produced from thermal expansion, 𝜏̿: ∇𝑢⃑ describes energy produced from viscous heating, 𝜌 𝐻 describes internal heat production, and −𝑇 𝐷𝑆∗ 𝐷𝑡 represents energy change due to additional entropy changes of other processes. Equation (2.1) can be rewritten as: ∂ρ ∂t + 𝜌∇ ∙ 𝑢⃑ + 𝑢⃑ ∙ ∇𝜌 = 0……………………………………………………………Equation (2.4) For an incompressible fluid, because 𝐷𝜌 𝐷𝑡 = ∂ρ ∂t + 𝑢⃑ ∙ ∇𝜌 = 0, Equation (2.4) is equivalent to: ∇ ∙ 𝑢⃑ = 0…………………………………………………………………………….Equation (2.5) To better understand the role perturbations play in a dynamic process, density, pressure, and gravitational acceleration are often decomposed into a reference value and a perturbation value: 𝜌 = 𝜌𝑟𝑒𝑓 + 𝜌̃ 𝑝 = 𝑝𝑟𝑒𝑓 + 𝑝̃ 𝑔 = 𝑔 𝑟𝑒𝑓 + 𝑔 ̃……………………………………………………………………….Equations (2.6) where variables ending with “ref” are reference values, and variables capped with “~” are perturbation values. Using the decomposition equations (2.6), the conservation equations for incompressible fluid are converted to: 18 ∇ ∙ 𝑢⃑ = 0…………………………………………………………………………….Equation (2.7) 𝜌 𝐷𝑢⃑⃑ 𝐷𝑡 = −∇𝑝̃ + ∇ ∙ 𝜏̿ + 𝜌̃ 𝑔 𝑟𝑒𝑓 + 𝑔 ̃(𝜌𝑟𝑒𝑓 + 𝜌̃)………………………………………Equation (2.8) 𝜌𝐶𝑝 𝐷𝑇 𝐷𝑡 = ∇ ∙ 𝑘∇T + α𝑇𝜌𝑟𝑒𝑓(𝑢⃑ ∙ 𝑔 𝑟𝑒𝑓) + α𝑇 𝐷𝑝̃ 𝐷𝑡 + 𝜏̿: ∇𝑢⃑ + 𝜌 𝐻 − 𝑇 𝐷𝑆∗ 𝐷𝑡 …………….Equation (2.9) In Equation (2.9), the last entropy term is usually neglected or otherwise specified. To better scale the mantle convection system and reduce necessary parameters that characterize the convection process, we first non-dimensionalize the variables in the conservation equations in Equation (2.7) – Equation (2.9) with the following scaling laws: ρ = 𝜌0ρ′ α = α0α′ g = g0g′ κ = κ0κ′ η = η0η′ 𝐶𝑝 = 𝐶𝑝0 𝐶𝑝′ 𝑢⃑ = κ0 ℎ 𝑢⃑ ′ T = ∆𝑇T′ 𝜏̿ = κ0η0 ℎ2 𝜏̿′ 𝐻 = 𝐶𝑝0 κ0∆𝑇 ℎ2 𝐻′ 19 ∇= 1 ℎ ∇′ t = ℎ2 κ0 t′ 𝐷 𝐷𝑡 = κ0 ℎ2 𝐷 𝐷𝑡′ 𝑝̃ = κ0η0 ℎ2 𝑝̃′………………………………………………………………………...Equations (2.10) where 𝜌0, α0, g0, κ0, η0, and 𝐶𝑝0 are reference density, reference thermal expansivity, reference gravity, reference thermal diffusivity, reference viscosity, and reference heat capacity, respectively. Variables with a prime are non-dimensional variables. 𝜂 is viscosity. h is mantle thickness. ∆𝑇 is the temperature contrast between the surface and the CMB. After non- dimensionalization, conservation equations for incompressible fluid are converted to: ∇′ ∙ 𝑢⃑ ′ = 0………………………………………………………………………….Equation (2.11) 1 𝑃𝑟 ρ′ 𝐷𝑢⃑⃑ ′ 𝐷𝑡 = −∇′𝑝̃′ + ∇′ ∙ 𝜏̿′ + 𝑅𝑎𝜌̃𝑔 𝑟𝑒𝑓′ + 𝑅𝑎 𝛼0∆𝑇 𝑔 ̃′(𝜌𝑟𝑒𝑓′ + 𝛼0∆𝑇𝜌̃′)……………..Equation (2.12) 𝜌′𝐶𝑝′ 𝐷𝑇′ 𝐷𝑡′ = ∇′ ∙ 𝑘′∇′T′ + 𝐷𝑖α′𝑇′𝜌𝑟𝑒𝑓′(𝑢⃑ ′ ∙ 𝑔 𝑟𝑒𝑓 ′ ) + 𝐷𝑖𝛼0∆𝑇 𝑅𝑎 α′𝑇′ 𝐷𝑝̃′ 𝐷𝑡′ + 𝐷𝑖 𝑅𝑎 𝜏̿′: ∇′𝑢⃑ ′ + 𝜌′𝐻′………… ……………………………………………………………………………………..Equation (2.13) where the dimensionless parameters Pr, Di, and Ra are introduced to simplify the conservation equations and characterize the convection system more conveniently. Pr is the Prandtl number that describes the ratio of momentum diffusivity to thermal diffusivity. It is defined as: Pr = η0 𝜌0κ0 …………………………………………………………………………...Equation (2.14) 20 It measures the viscous resistance to inertia. After inserting the estimated reference viscosity, density, and thermal diffusivity of the Earth’s mantle, the Prandtl number is around the order of 1023. Thus, 1 𝑃𝑟 is close to zero for the Earth’s mantle. It means inertia is negligible in mantle convection, which is essentially a buoyancy-driven process. This is called the “infinite Prandtl number approximation”. Because of this, in mantle convection models, the left side of the conservation of momentum Equation (2.12) is often assumed to be zero. Di is the dissipation number that describes how effective mechanical energy (i.e., due to self-compression) is converted to heat. It is defined as: 𝐷𝑖 = α0g0ℎ 𝐶𝑝0 …………………………………………………………………………..Equation (2.15) For incompressible fluid, Di is assumed to be zero. Ra is the Rayleigh number that describes how vigorous a convection system is. It is defined as: 𝑅𝑎 = 𝜌0𝑔𝛼0∆𝑇ℎ3 𝜂0𝜅0 …………………………………………………………….….…..Equation (2.16) To further simplify the right of Equation (2.12), we compare the last two force related terms. The 𝑅𝑎𝑔 ̃′𝜌̃′ part in the last term is to the orders of magnitude smaller than 𝑅𝑎𝜌̃𝑔 𝑟𝑒𝑓′ and 𝑅𝑎𝑔 ̃′𝜌𝑟𝑒𝑓′, thus neglected to the first order. Under the Boussinesq approximation, we reserve the density perturbation term 𝑅𝑎𝜌̃𝑔 𝑟𝑒𝑓′ to drive the flow by buoyancy, thus this term is reserved in Equation (2.12). The gravity perturbation term, practically, was assumed to be much smaller than the second term and avoided. But actually, the effect of gravity perturbation could be approximately 10%-50% of the density perturbation (e.g., Ricard, 2015). From perspective of accuracy, if density perturbation term is included, the gravity perturbation term should be included. The gravity perturbation term can be important for particular research topics, such as 21 CMB dynamic topography (e.g., Lassak et al., 2007), but is beyond the scope of this dissertation. After omitting gravitational perturbation and combining the infinite Prandtl number assumption, the conservation of momentum Equation (2.12) is converted into: 0 = −∇′𝑝̃′ + ∇′ ∙ 𝜏̿′ + ℎ3 𝜂0𝜅0 𝜌̃𝑔 𝑟𝑒𝑓′…………………………………………………Equation (2.17) Under Boussinesq approximation, the dissipation number is zero. This converts conservation of energy Equation (2.13) into: 𝐷𝑇′ 𝐷𝑡′ = ∇′2T′ + 𝐻′…………………………………………………………………. Equation (2.18) For isochemical convection calculations, 𝜌̃′ is dependent on temperature: 𝜌̃ = −𝜌0𝛼0∆𝑇𝑇′…………………………………………………………………. Equation (2.19) For thermochemical convection calculations, 𝜌̃′ is dependent on temperature and composition: 𝜌̃ = ∆𝜌𝐶 − 𝜌0𝛼0∆𝑇𝑇′…………………………………………………………….Equation (2.20) where ∆𝜌 is the density difference between different compositions, and C is the compositional coefficient. Besides, the convection system is also constrained by the equation for chemical advection: ∂C ∂t + (𝑢⃑ ∙ ∇)𝐶 = 0………………………………………………………………….Equation (2.21) In addition, for Newtonian materials, which are widely applied to the Earth’s mantle in large-scale mantle dynamics, the deviatoric stress has a linear correlation with strain rate 𝜀̇ ̿: 𝜏̿ = 𝜂𝜀̇ ̿…………………………………………………………………………..….Equation (2.22) In reality, the mantle is more complicated with non-Newtonian rheology. Coupling non- Newtonian and Newtonian rheology is computationally more challenging and not always applied 22 to geodynamic calculations. When studying lithosphere-mantle dynamics and upper mantle dynamics, for example, where dislocation creep is the evident deformation mechanism, non- Newtonian rheology cannot be ignored and this equation should be modified, taking into account non-Newtonian rheology. When studying the lower mantle, because the majority of the lower mantle is seismically isotropic (e.g., F. Niu & Perez, 2004), it is often assumed that the dominating creep mechanism in the lower mantle is diffusion creep. Thus, non-Newtonian rheology is not considered in this dissertation. But we should note that the Earth’s lower mantle might have a combination of diffusion (employed with Newtonian rheology) and dislocation (employed with non-Newtonian rheology) creep mechanisms. In summary, when dropping all the primes and plugging Equation (2.20) into Equation (2.17), for thermochemical mantle convection models, we have the following non-dimensional governing conservation equations under infinite Prandtl number, static gravitational field, and Boussinesq approximation: ∇ ∙ 𝑢⃑ = 0…………………………………………………………………………...Equation (2.23) 0 = −∇𝑝̃ + ∇ ∙ 𝜂𝜀̇ ̿ + 𝑅𝑎(𝐵𝐶 − 𝑇)…………………………………………….…...Equation (2.24) 𝐷𝑇 𝐷𝑡 = ∇2T + 𝐻…………………………………………………………………. ….Equation (2.25) where B is the buoyancy ratio and is defined as: 𝐵 = ∆𝜌 𝜌0𝛼0∆𝑇 …………………………………………………………….….…..Equation (2.26) These are the equations governing models in Chapter 4 and Chapter 5. For isochemical mantle convection models, The C in Equation (2.24) is zero, as applied to in Chapter 3. 23 2.2 Numerical Methods For the studies investigated in this dissertation, we use the finite element code CitcomCU (e.g., Zhong, 2006) to solve the approximated conservation equations as derived in Chapter 2.1. Depending on different research problems to probe, we modify the code’s boundary conditions, rheological formulations, initial conditions, etc. for calculations in this dissertation. Typically, initial temperature condition, velocity and temperature boundary conditions, and viscosity formulations are prescribed before solving the conservation equations computationally. Usually, free-slip and isothermal top and bottom boundary conditions are applied in mantle dynamic calculations. Details of the modifications are described in Chapters 3, 4, and 5. For Chapters 4 and 5, we incorporate composition-dependent rheology into the code to investigate how composition- dependent rheology affects the viscous coupling and interactions among LLSVP materials, ULVZ materials, and the background mantle. All of the calculations in Chapter 3 are isochemical calculations performed in 3D partial spherical domains. All of the calculations in Chapter 4 and Chapter 5 are thermochemical calculations performed in 2.5 D spherical domains. In 2.5 D spherical models, model domains are composed of a layer with one-element thickness. By selecting this model geometry, we make use of the computational efficiency of 3D CitcomCU and obtain a great approximation of 2D annulus model geometries. The 2.5 D geometry enables us to investigate high-resolution cross-sectional shapes of LLSVPs and ULVZs. Compositional fields are tracked using ratio tracer method (e.g., Tackley & King, 2003). In Chapter 4, because the thickness of ULVZ materials is very small, resolving their geometries can be super computationally expensive if imposing uniform ultra-high resolution across depth. To improve computational efficiency, we construct meshes with lower 24 resolution near the surface and gradually increase the resolution of meshes as depth increases, so that we have the highest resolution in the lowermost mantle to investigate ULVZ materials. In Chapter 4 and Chapter 5, LLSVPs are hypothesized to be caused by denser thermochemical piles formed from Early Earth’s differentiation. To formulate such thermochemical piles, we initially prescribed a layer with higher density than the background mantle at the lowermost mantle, assuming this is the layer piles formed from after the Earth’s differentiation. Such layer approximately has a similar volume fraction of the mantle as that of the LLSVPs observed in seismic tomography models. Then, we let the layer evolve in the model to obtain thermochemical piles, assuming that piles eventually form from the mixing process in mantle convection. This method is just the way of formulating thermochemical piles in this dissertation under multiple assumptions and depending on the goal of the two chapters. 25 REFERENCES Gassmöller, R., Dannberg, J., Bangerth, W., Heister, T., & Myhill, R. (2020). On formulations of compressible mantle convection. Geophysical Journal International, 221(2), 1264–1280. https://doi.org/10.1093/gji/ggaa078 Lassak, T. M., McNamara, A. K., & Zhong, S. (2007). Influence of thermochemical piles on topography at Earth’s core–mantle boundary. Earth and Planetary Science Letters, 261(3– 4), 443–455. https://doi.org/10.1016/J.EPSL.2007.07.015 Niu, F., & Perez, A. M. (2004). Seismic anisotropy in the lower mantle: A comparison of waveform splitting of SKS and SKKS. Geophysical Research Letters, 31(24), 1–4. https://doi.org/10.1029/2004GL021196 Ricard, Y. (2015). 7.02 - Physics of Mantle Convection. In G. Schubert (Ed.), Treatise on Geophysics (Second, pp. 23–71). Elsevier. https://doi.org/10.1016/B978-0-444-53802- 4.00127-5 Tackley, P. J., & King, S. D. (2003). Testing the tracer ratio method for modeling active compositional fields in mantle convection simulations. Geochemistry, Geophysics, Geosystems, 4(4). https://doi.org/10.1029/2001GC000214 Zhong, S. (2006). Constraints on thermochemical convection of the mantle from plume heat flux, plume excess temperature, and upper mantle temperature. Journal of Geophysical Research: Solid Earth, 111(4), 4409. https://doi.org/10.1029/2005JB003972 26 CHAPTER 3: LOWERMOST MANTLE FLOW PATTERNS CAUSED BY SLABS OF DIFFERENT STRENGTH AND TRENCH LENGTH IMPINGING ON THE CORE-MANTLE BOUNDARY 3.1 Abstract The geometry of slabs within the upper mantle has been relatively well-imaged by regional tomography and seismic studies; however, slab deformation in the lower mantle remains poorly understood. Although tomography models reveal that the lower mantle beneath paleo-subduction regions are faster-than-average, the resolution is not high enough to image how slabs are actually deforming there. Geodynamical studies of lower mantle slabs have shown a wide range of possible deformational behaviors, ranging from stiff, buckling slabs to more-ductile accumulations of slab material undergoing pure shear (e.g., Loubet et al., 2009). Observations of seismic anisotropy due to crystal preferred orientation (CPO) may ultimately provide valuable insight into the manner of slab deformation in the lower mantle; however, the parameters and processes that cause CPO are not well-constrained. The development of CPO depends upon the intrinsic anisotropy and slip systems of lower mantle minerals, under the influence of the mantle flow field (i.e., strain). As a first step toward better understanding how CPO is generated in lowermost mantle slab regions, we employ 3D spherical geodynamical calculations to investigate how slab strength and trench length influence lowermost mantle lateral flow patterns. We find that strong slabs dominantly lead to stiff, buckling slabs in the deep mantle, resulting in mostly unidirectional lateral flow patterns in the lowermost mantle. In contrast, weaker slabs are dominantly smooshed under pure shear, resulting in more-azimuthally divergent lowermost mantle lateral flow patterns. Furthermore, we find that decreasing trench length is effectively similar to decreasing slab strength in terms of resultant deformational behavior and lowermost mantle lateral flow patterns. In summary, strong 27 and/or wide slabs are expected to have different lowermost mantle flow patterns (more- unidirectional flow) than weak and/or narrow slabs (more azimuthally-divergent flow). This difference will likely result in a characteristically different directionality to the CPO that would develop, possibly allowing observations of seismic anisotropy to provide constraints on slab strength and slab deformational styles in the lower mantle. 3.2 Introduction Seismic tomography models (e.g., Grand, 2002; Koelemeijer et al., 2016; C. Li et al., 2008; Ritsema et al., 2004; Van Der Hilst et al., 1997; D. Zhao, 2004) and regional waveform studies (e.g., Hutko et al., 2006) reveal some continuous faster-than-average seismic velocity anomalies extending into the lower mantle from the upper mantle and high seismic velocity patches at the bottom of the mantle. The observations provide strong evidence that subducted lithosphere slabs can penetrate the Earth’s lower mantle and reach the core-mantle boundary (CMB). A better understanding of how slabs deform in the lower mantle may lead to a clearer view of large-scale mantle convection patterns and dynamics in the deep Earth. Subducted slabs drive plate tectonics and mantle convection (e.g., Conrad & Lithgow-Bertelloni, 2002; Forsyth & Uyedaf, 1975; Zahirovic et al., 2015). Styles of slab deformation are viscously coupled to the background mantle (e.g., Z. H. Li et al., 2019; W. Mao & Zhong, 2018). A deeper understanding of the deformation slabs currently undergo in the Earth’s mantle provides significant information on the mixing efficiency of mantle materials (e.g., Goes et al., 2017). Besides, deformational modes of slabs in the lower mantle will determine characteristics of thermal and chemical heterogeneities in the lowermost mantle. When slabs descend into the lower mantle, they induce large-scale convective currents and introduce thermal and compositional heterogeneities into the deep mantle (e.g., M. Li et al., 2014; Tan et al., 2002). Such heterogeneities can cause several dynamic consequences. As 28 cold slab or slab remnants reach the CMB, they cause lateral temperature variations directly into this thermal boundary layer and promote the cooling of the core. Thermal conditions at the CMB may change the magnetic reversal frequency and will likely affect the geodynamo (e.g., Larson and Olson, 1991). In the meantime, as slabs descend to the CMB, they carry dense basaltic component, whose segregation from slabs are hypothesized to be one of the mechanisms that cause the accumulation of dense thermochemical piles. The accumulation of such compositional reservoir is a candidate for the origin of the Large Low Shear Velocity Provinces (LLSVPs) at the base of the mantle (e.g., Brandenburg & van Keken, 2007; Christensen & Hofmann, 1994; J. Huang & Davies, 2007; M. Li & McNamara, 2013; Tackley, 2011). The composition and origin of such piles affect CMB heat flux, cycling of basaltic materials, and the mantle’s trace element history (e.g., M. Li et al., 2014). Characterizing the deformation of slabs in the lowermost mantle from seismic observations is important but challenging. Firstly, slab geometry in the lower mantle is not well-constrained. Although the advance of seismic tomography provides evidence of slabs penetrating the mantle transition zone and extending to a depth as deep as 1600 km, slabs often appear as broader structures with indistinct boundaries in the Earth’s lower mantle (e.g., Domeier et al., 2016; Fukao et al., 2001; Fukao & Obayashi, 2013; Hosseini et al., 2020; C. Li et al., 2008; Simmons et al., 2012) in the models. The blurring of slab geometries in seismic tomography models could be due to dynamic effects such as slab buckling in the lower mantle (e.g., Ribe et al., 2007) or the limitations in seismic imaging itself (e.g., Ritsema et al., 2007; Shephard et al., 2017). Recently, regional studies analyzing pre- cursors and/or post-cursors of seismic waves provide more detailed information on lower mantle slabs. It is found that some scatters clustering around the depth of mid- and lower mantle (e.g., 29 Bentham & Rost, 2014; Hutko et al., 2006; Kito et al., 2008; Rost et al., 2008) coincide with high seismic velocity anomalies mapped from tomography models or areas of present subduction. These scatters were thus interpreted to be along slab-mantle boundary and related to deep subduction. In Hutko et al., 2006, a step-like structure outlined by seismic scatters at the base of the mantle is interpreted to be caused by the folding of a cold slab. Although scattered signals from basaltic structures in the lower mantle are hard to detect (e.g., Kaneshima & Helffrich, 1999; F. Niu et al., 2003), with more evidence from seismic scattering studies, highly resolved images of lowermost mantle slabs may be obtained. Secondly, although the presence of significant and laterally varying seismic anisotropy at the base of the mantle is generally considered to be related to lower mantle deformation (e.g., Deng et al., 2017; X. He & Long, 2011; Long, 2009; Lynner & Long, 2014), how to correlate D” seismic anisotropy observations to the deformation is unclear (e.g., Lin et al., 2014; Romanowicz & Wenk, 2017). In D” regions with faster than average seismic velocities, horizontally polarized shear waves generally travel faster than vertically polarized ones (e.g., Asplet et al., 2020; Auer et al., 2014; Chang et al., 2014; French & Romanowicz, 2014; Kustowski et al., 2008; Nowacki, 2013; Panning & Romanowicz, 2006; Walpole et al., 2017). Since fast regions are generally caused by subduction or paleo subduction (e.g., Garnero & Lay, 2003), seismic anisotropy in fast regions is naturally thought to be likely related to slab deformation in the lower mantle (e.g., Cottaar et al., 2014; Faccenda & Capitanio, 2012). Unlike how straightforward upper mantle anisotropy reflects deformation in the olivine system (e.g., Becker, 2008; Karato et al., 2008), linking lowermost mantle seismic anisotropy observations to slab deformation is not easy. The difficulty, on the one hand, is due to challenges in lowermost mantle anisotropy observations, such as inadequate sampling and contamination in the upper mantle (e.g., Cottaar et al., 2014; Romanowicz & Wenk, 30 2017). On the other hand, to link seismic anisotropy to deformation, both steps of linking anisotropy to flow and linking flow to mantle deformation are required but not well understood. The step of linking anisotropy to lowermost mantle flow patterns is challenging because of the nonuniqueness of shear wave splitting measurement interpretations (e.g., Wolf et al., 2019), imperfectly understood lowermost mantle mineralogy, and unclear deformation mechanisms at CMB conditions (e.g., Creasy et al., 2020; Jackson & Thomas, 2021; Long, 2013). The deformation mechanism is unclear because such lab experiments cannot be done at realistic CMB temperatures or strain rates. The strong seismic anisotropy in faster than average D” regions related to subduction or paleo subduction is hypothesized to be caused by lattice-preferred orientation (LPO) (e.g., Wenk et al., 2006, 2011), if dislocation creep dominates such regions (e.g., McNamara et al., 2003). The major lower mantle minerals post-perovskite (pPV), magnesiowüstite [(Mg, Fe)O], and bridgmanite are all intrinsically anisotropic. They can therefore contribute to an LPO fabric that causes seismic anisotropy in the deep mantle. The LPO hypothesis raised particular interest with the discovery of post-perovskite (pPV) (e.g., Merkel et al., 2007). The phase transition from perovskite to pPV can occur in cold regions at the D” layer depth (e.g., Miyagi et al., 2010; Murakami et al., 2004; Tsuchiya et al., 2004; Grocholski et al., 2012). pPV can develop LPO in shear flow at the base of the mantle (e.g., Cottaar et al., 2014; Wu et al., 2017), possibly induced by subducted slabs impinging on the CMB. In addition, subduction/paleo subduction regions are generally dominated by dislocation creep (e.g., McNamara et al., 2001), which provides a proper deformation environment for pPV to develop LPO (e.g., Wenk et al., 2011). Besides pPV, magnesiowüstite would develop LPO, leading to D” azimuthal anisotropy (e.g., Long et al., 2006). Therefore, the development of LPO of post-perovskite, (Mg,Fe)O, bridgmanite or a combination 31 of them in the deformation of slabs impinging on CMB is further hypothesized to cause D” anisotropy in fast regions. Under the assumption that the development of LPO in the lowermost mantle causes the D” seismic anisotropy, recent studies (e.g., Cottaar et al., 2014; Wenk et al., 2011) made progress on the first step of establishing the link between D” anisotropy measurements and lowermost mantle flow patterns introduced by slab deformation in the lower mantle. They performed forward modeling that combines dynamic flow models, shear wave splitting measurements, and mineral physics tools in a subducting slab. They also compared synthetic results with seismic observations to find best-fit flow patterns. However, the other step required to eventually correlate slab deformation and seismic anisotropy: linking slab deformation and flow is significant but not fully investigated. Earlier laboratory experiments and geodynamic modeling have proposed a dynamical richness of possible slab deformational styles (e.g., Christensen & Hofmann, 1994; Loubet et al., 2009; Schellart, 2008; Tan et al., 2002), which provides a fundamental understanding of what physically controls lower mantle slab deformational styles. (Loubet et al., 2009; Schellart, 2008) show that how the lower mantle slab deforms depends on the viscosity contrast between the slab and the background mantle (slab strength). Slab strength controls the amount of internal deformation a slab undergoes (e.g., Schellart, 2008). Based on the amount of internal deformation within slabs, deformation modes of slabs in the lower mantle can generally be categorized into three types (e.g., Loubet et al., 2009): 1) Straight descent, (Figure 3.1a) where slabs undergo simple shear and descend through the lower mantle with minimal internal deformation as they reach the CMB. 32 2) Buckling and folding, (Figure 3.1b) where slabs undergo rigid buckling or folding with more internal deformation. 3) Squishing under pure shear, (Figure 3.1c) where slabs undergo intense pure shear deformation and are squished into a glob of subducted material as they reach the CMB. Besides slab strength, the tabular geometry of slabs also affects the deformation modes of slabs in the deep mantle (e.g., Guillaume et al., 2021; Loubet et al., 2009). Calculations in 3D spherical geometry give us more view angles on slab 3D geometry and a more comprehensive understanding of slab morphologies. For instance, in 3D spherical geometries, slabs with very wide trench length can develop a “wrinkle-like” concave curvature in slab edges reaching the CMB (e.g., Morra et al., 2009). Another example is that although plumes are typically triggered at slab edges in 2D models, plumes are more irregularly distributed around slab edges in 3D models (e.g., Tan et al., 2002). To ultimately understand slab deformation in the lower mantle through anisotropy observation, it is critical to link between the lowermost mantle flow patterns and slab deformational styles. In this scope, we hereby investigate how slab strength and trench length (as one of the dimensions composing slab tabular geometry), affect lowermost mantle flow patterns. Incoherent with the slab deformation styles in the lower mantle as categorized in the previous paragraph, we hypothesize that 3 categories of deformation lead to 3 distinctive lowermost mantle lateral flow patterns: 1) Unidirectional, (Figure 3.1a) where the lowermost mantle flow induced by a strong slab is aligned with one direction. Depending on transient dynamical scenarios and rheological parameters, such direction is either aligned with or opposed to the plate motion direction. 33 2) Bilaterally-divergent, (Figure 3.1b) where lowermost mantle flow induced by a moderate strength slab diverges at the projection of slab fold axis near the CMB, and along the plate motion direction. 3) Azimuthally-divergent, (Figure 3.1c) where the lowermost mantle flow induced by a weak slab uniformly spreads out in all directions. Furthermore, we hypothesize that slabs with longer trenches would induce less azimuthally divergent flow at the base of the mantle because the map-view geometry of slabs in the lowermost mantle is constrained by their horizontal cross-sectional shape when subducted. With a longer trench length, the slab is subducted approximately as a tabular sheet, while with a shorter trench length, the slab is subducted approximately as a column. The map-view geometry of slabs constrains the effective directionality of the lowermost mantle flow patterns slabs induce. Here, we perform geodynamical calculations in 3-D partial-sphere geometry to systematically investigate whether and how slab strength and trench length control lateral flow patterns in subduction regions of the lowermost mantle. We design our experiment to generate continuous slabs that straightly descend through the mantle transition zone and impinge on the CMB. We analyze the directional symmetry of lateral flow in subduction regions at the base of the mantle. 34 Figure 3.1 Conceptual models of three categories of slab deformation in the lower mantle with a) low internal deformation, b) moderate internal deformation, and c) large internal deformation and the hypothesized lateral flow pattern beneath them, respectively. Blue in a), b), and c) represents slabs and arrows on the slabs represent the movement direction of the slab. Dashed contours and black arrows in the rectangular boxes in a), b), and c) denote map view of cold regions beneath slabs and lateral flow patterns, respectively. Ideally, most lateral flow would be unidirectional, biaxially perpendicular to the slab axis, and azimuthally divergent for strong a) slab, b) moderate strength slab, and c) weak slab, respectively. Asymmetry parameters of flow patterns in a), b), and c) would be 1, 0, and 0, respectively, with 1 meaning the flow pattern is not symmetrical by the slab axis and 0 meaning the flow pattern is perfectly symmetrical by the slab axis. Azimuthal divergence parameter of flow patterns in a), b), and c) would be 0, 0, and 1, respectively, with 1 meaning the flow pattern is not azimuthally divergent at all and 0 meaning the flow pattern is perfectly azimuthally divergent. Refer to the Method section for the definition of the parameters. 35 3.3 Methods The geodynamic calculations are performed by solving the following non-dimensional equations derived from conservation of mass, momentum, and energy under the Boussinesq approximation: ∇ ∙ 𝑢⃑ = 0…………………………………………………………………………….Equation (3.1) −∇𝑃 + ∇ ∙ (𝜂𝜀̇ ̿) = 𝜆𝑅𝑎𝑇𝑟̂…………………………………………………………..Equation (3.2) 𝜕𝑇 𝜕𝑡 + (𝑢⃑ ∙ ∇)𝑇 = ∇2𝑇 + 𝑄…………………………………………………………..Equation (3.3) Here, 𝑢⃑ , P, η, 𝜀̇ ̿, T, 𝑟̂ , t and Q represent the velocity, dynamic pressure, viscosity, strain rate, temperature, the unit vector in the radial direction, time, and the volumetric internal heating rate. 𝜆 = (𝑅/ℎ)3 with R as the radius of the Earth and h as the mantle thickness. Ra is the Rayleigh number, which is defined as: 𝑅𝑎 = 𝜌0𝑔𝛼0∆𝑇ℎ3 𝜂0𝜅0 …………………………………………………………….….……Equation (3.4) where 𝜌0 , 𝛼0 , 𝜂0 , 𝜅0 are the dimensional reference values of density, thermal expansivity, viscosity at non-dimensional temperature T = 0.5, and thermal diffusivity. ∆𝑇 is the dimensional reference temperature contrast between the surface and the CMB. g is the constant of gravitational acceleration. Although there is large uncertainty in these parameters, we apply values typically employed in mantle convection problems, as listed in Table 3.1. We use a Rayleigh number of Ra = 1.0 × 107 in the reference case (Case 1) and most of the cases we are interested in. The variant parameters used in all the cases are listed in Table 3.2. To focus on the fundamental deformational styles of slabs in the lower mantle and better facilitate vertically sinking slabs in the upper mantle, we do not include phase transition at 660 km depth and chemical variation here. The complexity 36 added from phase transition and chemical variation to the mantle flow pattern near the CMB is therefore excluded, so that we can clearly observe how lowermost mantle flow patterns respond to different deformation mode of slabs in the lower mantle. Because we do not include chemical variation, when we refer to “slab”, we mean “cold anomaly” with -0.1 non-dimensional residual temperature in this study. The viscosity 𝜂 is temperature- and depth-dependent: 𝜂(𝑇, 𝑧) = 𝜂𝑧(𝑧)exp [𝐴(0.5 − 𝑇)]………………………………………………….Equation (3.5) Here 𝜂𝑧(𝑧) is 1 and 30 for the upper and lower mantle, leading to a 30x viscosity increase from the upper mantle to the lower mantle. A and T are the non-dimensional activation coefficient and temperature, respectively. A controls the viscosity contrast between the coldest (at T = 0 on the surface) and the hottest (at T = 0.5 on the CMB) regions in the model. We investigate cases with activation coefficients of 18.4207, 13.8155, 9.2103 and 4.6052, which lead to a 104 ×, 103 × , 102 ×, 𝑎𝑛𝑑 10 × viscosity contrast due the temperature differences between the hottest and the coldest regions, respectively. As defined by equation (3.5), greater viscosity contrast (i.e., greater A) leads to greater viscosity in the coolest regions in the model. If we consider slabs as cold anomalies, greater A effectively leads to stronger and stiffer subducted slabs. 37 Table 3.1 Parameters used in this study. Parameter Description ∆𝑇 𝛼0 𝜌0 𝑔 R h 𝜂0 𝜅0 𝜂𝑧 Q potential temperature drop across mantle reference thermal expansivity reference density gravitational constant radius of the Earth mantle thickness reference viscosity reference thermal diffusivity viscosity contrast between the upper and lower mantle non-dimensional (scaled by Earth radius) internal heating rate a Vplate plate velocity Value 2600 3 × 10−5 3340 9.8 6370 2870 Units K K-1 kg m-3 m s-2 km km 6.03 × 1021 Pa s 1 × 10−6 m2 s-1 30 60 4 - - cm yr-1 a The plate velocity listed in the table is the plate velocity at the equator. For a place at colatitude ϴ, the plate velocity should be computed as v = Vplate × 𝑠𝑖𝑛𝜃. 38 Table 3.2 Cases investigated in this study. Case Ra A trench length (km) Viscosity contrasta 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 1.0 × 107 1.0 × 107 1.0 × 107 5.0 × 107 1.0 × 107 1.0 × 107 1.0 × 107 5.0 × 107 1.0 × 107 1.0 × 107 1.0 × 107 5.0 × 107 9.2103 3,336 4.6052 3,336 13.8155 3,336 18.4207 3,336 4.6052 1,112 9.2103 1,112 13.8155 1,112 18.4207 1,112 4.6052 7,784 9.2103 7,784 13.8155 7,784 18.4207 7,784 5.0 × 106 4.6052 1.0 × 107 9.2103 445 445 1.0 × 107 13.8155 445 100 10 1,000 10,000 10 100 1,000 10,000 10 100 1,000 10,000 10 100 1,000 a The viscosity contrast between the coldest region (surface) and the hottest region. The conservation equations are solved using our modified version of CitcomCU (e.g., Zhong, 2006). The models are formulated in a regional 3-D spherical domain, which encompasses 1/6 of a whole spherical shell. The longitude ranges from 30° to 150°; the colatitude ranges from −90° to 90°. The model consists of 192, 128 and 64 elements in longitude, colatitude, and radial directions, respectively. Temperature boundary conditions are insulating on the side, and isothermal on the surface and bottom. The model is heated from both below and internally, with a non-dimensional internal heating rate Q = 60 (scaled by mantle radius). To generate a descending slab, we apply regional kinematic boundary condition on the surface. A northward constant velocity of 4 cm/year is imposed on surface region from −90° to 0° longitude and 75° to 105° 39 colatitude (Figure 3.2), so that slabs with trench length of 3,336 km can sink along 0° longitude. The remaining surface area is stationary, having zero velocity. Except for the surface, the other velocity boundary conditions are free slip. Slabs with trench lengths of 445, 1,112, and 7,784 km are also investigated. For each calculation, we initially prescribe a cold plate along the surface with uniform T = 0. The thickness of the plate is determined from the half space cooling model of oceanic plates. Because the goal of study is to investigate how slabs deform in the lower mantle and the lowermost mantle flow patterns induced by slabs, we want to isolate slab deformation, so that slab morphologies and flow patterns are not complicated by the presence of plumes. One way to achieve this is by equating the initial background mantle temperature with the temperature of the CMB. Here, we apply an initial temperature of T = 0.5 to the background mantle (Figure 3.3a). We then analyze lowermost mantle flow patterns after the system reaches thermal equilibrium. We fix the plate velocity at 4 cm/year, which is a typical present-day plate motion value. For slabs with different width, we tune the Rayleigh number separately. For slabs with a certain width but at different activation parameter A, we change the Rayleigh number to formulate slabs that penetrate directly through the transition zone and extend to the CMB as continuously and steadily as possible. We adjust the Rayleigh number to match the intrinsic convection speed with the subduction speed of slabs. If the two speeds do not match, slab deformation will be affected in different ways. When the convection speed is faster the slab will subduct in a manner of flat-slab subduction (Figure 3.8(a)). On the contrary, when the convection speed is slower, slabs will accompany large blob of droplets beneath them (Figure 3.8(b)). In reality, we do not know whether the intrinsic convection speed will match perfectly with the subduction speed of slabs and how will this influence slab deformation and flow patterns in the lower mantle. As complementary experiments, we systematically investigate for different As, the influence of Rayleigh number on 40 Figure 3.2 A slab model (blue isosurface) from Case1. The left panel shows the map view of the model. The black outline delineates the plate and the blue arrow denotes that the plate has a fixed angular velocity of 4 cm/year from south to north. The core is orange. The right panel shows temperature cross section at colatitude 90⁰ (orange line a-a’ on the left). Color denotes non- dimensional temperature. Figure 3.3 Cross-sections at a-a’ in Figure 3.2 of (a) initial temperature field and (b) temperature field after 500 Ma since the model initiates from (a). 41 slab deformation in the lowermost mantle and flow patterns near the CMB. For each different A value selected, we fix A at that value, and firstly run the case with a low Rayleigh number around 1.0 x 106. Then we keep doubling the Rayleigh number and run the case, until the Rayleigh number is as high as around 1.0 x 108. For different A value, starting low Ra value and ending high Ra value are not always the same, since different A leads to a different temperature-dependent viscosity due to Equation (3.5), thus a different rheology. To estimate the lateral flow symmetry beneath slabs in the lowermost mantle, we focus on the lateral flow distribution within the colder-than-average region at the depth of 2800 km (called “the region” in the context). Assuming that the net flow direction can represent the direction where most flow is directed, we can categorize the symmetry of the lateral flow within the region based on it. To estimate the symmetry character of the flow patterns by the slab axis, we define the following asymmetry parameter: 𝛼 = ∑(𝑣𝜏⃑⃑⃑⃑ ∙ 𝑢⃑⃑ ) ∑(𝑣𝜏⃑⃑⃑⃑ ∙ 𝑢⃑⃑ ) 𝑤ℎ𝑒𝑛 (𝑣𝜏⃑⃑⃑⃑ ∙ 𝑢⃑⃑ )> =0.0 ……………………………………………….……..….Equation (3.6) where 𝑣𝜏⃑⃑⃑ is the lateral velocity at each element node in the region and 𝑢⃑ is the unit vector (also called “maximum direction”) of the net flow within the region. The numerator of the parameter is the sum of flow towards the maximum direction minus the sum of flow against the maximum direction, while the denominator of the parameter is the sum of flow towards the maximum direction. If all the flow is towards the maximum direction, the flow is extremely asymmetric by the slab axis and we will have α = 1.0. On the other hand, if the flow is evenly distributed by the slab axis, we will have α = 0.0. To estimate the azimuthal divergence of the flow patterns, we define the following azimuthal asymmetry parameter: 𝛽 = | ∑|𝑣𝜏⃑⃑⃑⃑ ∙ 𝑢⃑⃑ |− ∑ |𝑣𝜏⃑⃑⃑⃑ ∙ 𝑢⊥⃑⃑⃑⃑⃑ || ∑ |𝑣𝜏⃑⃑⃑⃑ ∙ 𝑢⃑⃑ | ………………………………………….………….……..….Equation (3.7) 42 where 𝑢⊥⃑⃑⃑⃑ is the unit vector that is orthogonal to the maximum direction. The numerator of the parameter is the difference between the sum of the flow magnitude towards or against the maximum direction and the sum of the flow magnitude orthogonal to the maximum direction. The denominator is the sum of the flow magnitude towards or against the maximum direction. If the flow is symmetric azimuthally, we will have β = 0.0. However, if the flow is biased parallel to the maximum direction, we will have β = 1.0. To estimate the time average of the asymmetry character by the slab axis and the azimuthal asymmetry, we integrate these parameters by time and average them by time for each calculation so that the integral is normalized to the range of 0 to 1. 3.4 Results 3.4.1 Lowermost mantle flow patterns induced by slabs as a function of slab strength Taking slabs of a trench length of 3,336 km as a reference, we can see that weak slabs are smooshed under pure shear in the lowermost mantle (Figure 3.4(b)), creating three colder-than- average patches near the CMB; moderate strength slabs are buckling in the lower mantle before it reaches the CMB (Figure 3.4(e)); strong slabs bend in the lower mantle and shear along the CMB in one direction (Figure 3.4(i) and 3.4(g)). The lowermost mantle lateral flow velocities induced by weak slabs (Figure 3.4(a)) and moderate strength slabs (Figure 3.4(d)) look not significantly different and both are symmetric by the slab axis (colatitude 90°), while those induced by strong slabs are directed towards one side of the slab axis (Figure 3.4(i)). Weak and moderate strength slabs also induce flow spreading out in all directions, including parallel to the slab axis, while the flow strong slabs induce is dominantly perpendicular to the slab axis (Figure 3.4(a), 3.4(d) and 3.4(g)). As we quantify the time average of the asymmetry character, we consistently find that the asymmetry character by the slab axis (𝛼) and azimuthal asymmetry character (β) is small and close 43 for the lowermost mantle lateral velocity induced by weak and moderate strength slabs and increase by more than two times for that induced by strong slabs (Figure 3.5). For slabs with trench length of 1,112 km, 3, 336 km and 7,784 km, we find as slab strength increases, the dynamical behavior of slabs in the lower mantle changes from being smooshed under pure shear-buckling-bending and shearing along the CMB in one direction, and the lowermost mantle lateral flow pattern becomes less symmetrical by the slab axis and less azimuthally mantle lateral flow pattern becomes less symmetrical by the slab axis and less azimuthally divergent in general (Figure 3.9(a) and Figure 3.9(b)). α increases abruptly as A increases from 9.2103 to 13.8155 (Figure 3.9(c)), while α varies smaller than 0.1 as A increases from 4.6052 to 9.2103 except that α is very high for weak slabs with 1,112 km trench length and drops as A increases to 9.2103. α increases to over 0.7 for all the slabs as A increases to equal to or larger than 13.8155 (Figure 3.9(c)), meaning that the directionality of the lowermost mantle lateral flow pattern is biased towards one side of the slab axis and has a unidirectional pattern beneath strong slabs. We can see from Figure 3.9(d) that the increase rate of β as a function of A decreases as the trench length increases. The absolute value of β ranges between 0.6 and 0.9 and only increases by 5% as the viscosity contrast between the CMB and the surface plate (“viscosity contrast”) increases by 1,000 times if the trench length of slabs is 7,784 km, while it has a larger range, between 0.16 and 0.71 and increases by 430% as the viscosity contrast increases by 1,000 times if the trench length of slabs is 1,112 km. These indicate that the degree of azimuthal-divergence decreases as slab strength increases, and the degree of azimuthal divergence becomes steadier and is less influenced by slab strength as the trench length increases. Combining the characteristics of α and β as a function of slab strength, we can conclude that, in general, the lateral flow directionality becomes 44 less symmetric by the slab axis and less azimuthally symmetric, and the influence of slab strength on the azimuthal symmetry of the flow pattern becomes less as the trench length increases. 45 Figure 3.4 Representative snapshots of lowermost mantle (120 km above the CMB) lateral flow velocity (black arrows) in cold regions and slabs as a function of slab strength (manipulated by A). a-c, snapshots at time = 2.92 billion years after it reaches a thermal equilibrium for a weak slab model (Case 2), where the surface plate has a 10x viscosity contrast with the CMB. d-f, snapshots at time = 3.33 billion years after it reaches a thermal equilibrium for a moderate strength slab model (Case 1), where the surface plate has a 100x viscosity contrast with the CMB. g-i, snapshots at time = 1.17 billion years after it reaches a thermal equilibrium for a strong slab model (Case 3), where the surface plate has a 1000x viscosity contrast with the CMB. Colder and hotter temperatures are represented by blue and red colors, respectively. Left column shows the map view of the lowermost mantle lateral flow velocity, with a background of temperature field. Top right graphs on each row shows the slab (blue isosurface) that induces the flow pattern on the left. Bottom right graphs on each row are temperature cross-sections at colatitude 90⁰ (blue line in (a), (d), and (h)). 46 Figure 3.5 Time average of asymmetry characters of lowermost mantle lateral velocity induced by slabs with a trench length of 3,336 km (30-degree arc angle) as a function of slab strength. Low 𝛼 means the lateral flow is symmetric by the slab axis, and low β indicates the lateral flow is symmetric azimuthally. Models with A of 4.61, 9.21, and 13.82 lead to 10x, 100x, and 1000x viscosity contrast between the surface plate and the CMB in the models 47 3.4.2 Lowermost mantle flow patterns induced by slabs as a function of trench length Taking slabs with A of 9.2103 as a reference, we discover that regardless of the different trench lengths, the slabs consistently buckle in the lower mantle (Figure 3.6(c), (e), and (h)), which means trench length does not evidently affect the dynamical behavior of slabs in the lower mantle. However, the increase of trench length makes the lowermost mantle lateral velocity induced by slabs to be less azimuthally symmetric (Figure 3.7). The increase in trench length does not change the symmetry of the lowermost mantle lateral velocity by the slab axis, though (Figure 3.7). When we look at slabs with greater (A of 13.85) and smaller internal deformation (A of 4.61), this pattern is true as well. As we can see in Figure 3.10 (a) and 3.10 (c), for the symmetry of slab lowermost mantle lateral velocity by the slab axis, strong slabs always have an α larger than 0.8 and close to 1 regardless of the change in trench length, meaning that the lateral velocity is always biased towards one side of the slab axis; moderate slabs have a steady α regardless of the change in trench length as well, except for showing an α=1.0 for slabs with a trench length of 445 km; weak slabs does not show a clear correlation between α and trench length, meaning that trench length might not be a dominant factor on the symmetry of slab lowermost mantle lateral velocity by the slab axis. Thus, we can conclude that trench length is not a dominant factor that affects the symmetry of the slab lowermost mantle lateral velocity by the slab axis. On the other hand, for the azimuthal symmetry of the lateral velocity, except for slabs of 445 km trench length, we see an increase in β, which means a decrease in azimuthal symmetry as the trench length increases in weak, moderate strength, and strong slabs (Figure 3.10 (b) and 3.10 (d)). Combining the characteristics of α and β as a function of trench length, we can draw to a conclusion that in general, the slab lowermost mantle lateral flow becomes less azimuthally symmetric as the trench 48 length increases, but the symmetry of lateral flow by the slab axis does not have a clear correlation with 49 Figure 3.6 Representative snapshots of lowermost mantle (120 km above the CMB) lateral flow velocity (black arrows) in cold regions and slabs as a function of trench length. a-c, snapshots at time = 1.54 billion years after it reaches a thermal equilibrium for a narrow slab model (Case 6). d-f, snapshots at time = 3.33 billion years after it reaches a thermal equilibrium for a moderate width slab model (Case 1). g-i, snapshots at time = 1.71 billion years after it reaches a thermal equilibrium for a wide slab model (Case 10). These three cases are identical, and all formulate moderate strength slabs except that they have a trench length of 1,112 km (Case 6), 3,336 km (Case1), and 7,778 km (Case 10). Colder and hotter temperatures are represented by blue and red colors, respectively. Left column shows the map view of the lowermost mantle lateral flow velocity, with a background of temperature field. Top right graphs on each row shows the slab (blue isosurface) that induces the flow pattern on the left. Bottom right graphs on each row are temperature cross-sections at colatitude 90⁰ (blue line in (a), (d), and (h)). 50 Figure 3.7 Time average of asymmetry characters of lowermost mantle lateral velocity induced by slabs with a 100x viscosity contrast between the surface plate and the CMB as a function of trench length. Low 𝛼 means the lateral flow is symmetric by the slab axis, and low β indicates the lateral flow is symmetric azimuthally. The lateral velocity is less azimuthal symmetric as trench length increases. the trench length. 3.5 Discussions We investigated a wide range of slabs with various temperature-dependent rheology and trench lengths. In 3D partial spherical geometry, we observed a richness of slab dynamical behaviors in the lower mantle. In general, slab dynamical behaviors fall into the three categories as a function of slab strength as we hypothesized (Figure 3.1), regardless of the trench length they have on the surface: weak slabs with low internal deformation are dominantly smooshed under shear in the lowermost mantle; moderate strength slabs with moderate internal deformation dominantly buckle in the lower mantle; strong slabs with great internal deformation prominently bend in the lower mantle and glide across the CMB in one direction. The trench length dimension, though, does add complexity to the slab geometry in the lower mantle. For weak and moderate strength slabs with greater than or equal to 3,336 km trench length, we frequently, if not dominantly, observed slabs having more than or equal to 2 rounded-shape cold-patched cores (i.e., Figure 3.4(a)) at the base of the mantle. Those cold patches follow the cylindrical downwelling bulges on the trench side of the slabs (Figure 3.4 (b) and Figure 3.11 (b)). If we track those cylindrical bulges all the way up to the surface plate, we will observe that they are 51 caused by cold anomalies in the surface plate (Figure 3.11 (c)-(f)). The cold anomalies are not applied at the start of the calculation but are evolved and occur dominantly after the model reaches thermal equilibrium. Those cold anomalies in the surface plate are sustained down the CMB within the slabs, making the slab regions beneath those cylindrical bulges colder than other parts of the slabs. As for narrower slabs, the anomalies might overlap or barely separate, so we do not expect those cold patches in the lowermost mantle within narrower slabs. In addition, we found that weak and moderate strength slabs with a 445 km trench length sometimes have a simple shearing behavior on slab edges along the CMB in the lower mantle (Figure 3.12). Because narrow slabs generate a small volume of cold anomalies in the mantle, they are warmed up by the surrounding mantle easily. When they are not strong enough, they will sometimes be driven by the surrounding flow. The simple shear along the CMB behaviors is driven by the surrounding flow. The 3D spherical slab models here show a large variety of slab geometries and slab dynamical behaviors in the lower mantle that cannot be illustrated with 2D slab models and are much more complicated than the 2D conceptual models that are widely constructed in earth science (i.e., Dolejš, 2016; Jiang et al., 2015). When we try to interpret seismological observations or geochemistry observations, we need to consider the richness of slab dynamical behavior. We found that we can roughly distinguish slab strength from the directional distribution of lateral flow patterns along the CMB. However, the effects of slab trench length on the flow patterns make it difficult to infer slab strength and deformational styles of slabs in the lower mantle from lowermost mantle flow patterns. For example, stronger slabs induce more unidirectional lowermost mantle flow patterns, but the trench length effect on flow patterns might counteract the slab strength effect if they have shorter trench lengths. In reality, subduction 52 zones have a large variety of trench lengths, ranging from 400 km to 7700 km according to the near surface subduction zone geometry depicted in slab 2.0 model (Hayes et al., 2018) using seismicity. In this study, we covered a wide range of trench lengths, from 300 km to 7700 km, and took slabs with a trench length of 3,336 km, which is close to the trench length of Kermadec subduction zone and central America subduction zone, as a reference. Therefore, the parameter space investigated here is sufficient to determine if the trench length effect would prohibit us from distinguishing slab strength by lowermost mantle flow patterns. To infer slab strength and, thus slab deformational styles in the lower mantle from seismic anisotropy observations, linking lowermost mantle flow patterns with slab strength and deformational styles in the lower mantle is only the first step. The second step of linking seismic anisotropy patterns and flow patterns requires assumptions of mineral phase assemblages and matching observation with synthetic anisotropy patterns, which involves a collaboration with mineral physicists and seismologists. Use the phase system from a previous study that built the second step (Cottaar et al., 2014) as an example to infer what seismic anisotropy pattern we might expect. Suppose the lowermost mantle seismic anisotropy is caused by CPO of post- perovskite. They found in seismic models that the fast axis direction of the split shear waves was linked to the direction of horizontal flow around slab edges in the lower mantle. Therefore, a unidirectional flow pattern can be linked to a seismic anisotropy pattern with fast axis aligned towards one direction. However, because the link is complicated with another factor, the trench length, it’s not worthwhile to link the flow patterns with seismic anisotropy patterns and eventually relate slab deformational styles in the lower mantle to seismic anisotropy patterns. 53 3.6 Conclusions Our findings show that we cannot directly infer slab strength or slab deformational styles in the lower mantle from lowermost mantle flow patterns, because of the complications slab trench length added to the lowermost mantle flow patterns. It is shown that slab strength controls slab deformational styles in the lower mantle and lowermost flow patterns, just as hypothesized. However, slabs with shorter trench lengths and weaker slabs can non-uniquely produce azimuthally divergent flow patterns, making it difficult to determine the dominant cause of a particular flow pattern beneath a subduction region. Therefore, although we established the link between lowermost mantle flow patterns and slab deformational styles in the lower mantle, because of the complications of the cause of lowermost mantle flow patterns, we should be cautious about inferring slab deformational styles in the lowermost mantle from the flow patterns and seismic anisotropy observations. 54 REFERENCES Asplet, J., Wookey, J., & Kendall, M. (2020). A potential post-perovskite province in D″ beneath the Eastern Pacific: evidence from new analysis of discrepant SKS–SKKS shear-wave splitting. Geophysical Journal International, 221(3), 2075–2090. https://doi.org/10.1093/GJI/GGAA114 Auer, L., Boschi, L., Becker, T. W., Nissen-Meyer, T., & Giardini, D. (2014). Savani: A variable resolution whole-mantle model of anisotropic shear velocity variations based on multiple data sets. Journal of Geophysical Research: Solid Earth, 119(4), 3006–3034. https://doi.org/10.1002/2013JB010773 Becker, T. W. (2008). Azimuthal seismic anisotropy constrains net rotation of the lithosphere. Geophysical Research Letters, 35(5). https://doi.org/10.1029/2007GL032928 Bentham, H. L. M., & Rost, S. (2014). Scattering beneath Western Pacific subduction zones: Evidence for oceanic crust in the mid-mantle. Geophysical Journal International, 197(3), 1627–1641. https://doi.org/10.1093/gji/ggu043 Brandenburg, J. P., & van Keken, P. E. (2007). Deep storage of oceanic crust in a vigorously convecting mantle. Journal of Geophysical Research: Solid Earth, 112(6). https://doi.org/10.1029/2006JB004813 Chang, S. J., Ferreira, A. M. G., Ritsema, J., Van Heijst, H. J., & Woodhouse, J. H. (2014). Global radially anisotropic mantle structure from multiple datasets: A review, current challenges, and outlook. Tectonophysics, 617, 1–19. https://doi.org/10.1016/J.TECTO.2014.01.033 Christensen, U. R., & Hofmann, A. W. (1994). Segregation of subducted oceanic crust in the convecting mantle. In JOURNAL OF GEOPHYSICAL RESEARCH (Vol. 99, Issue B10). https://doi.org/10.1029/93JB03403 Conrad, C. P., & Lithgow-Bertelloni, C. (2002). How mantle slabs drive plate tectonics. Science, 298(5591), 207–209. https://doi.org/10.1126/science.1074161 Cottaar, S., Li, M., McNamara, A. K., Romanowicz, B., & Wenk, H. R. (2014). Synthetic seismic anisotropy models within a slab impinging on the core–mantle boundary. Geophysical Journal International, 199(1), 164–177. https://doi.org/10.1093/GJI/GGU244 Creasy, N., Miyagi, L., & Long, M. D. (2020). A Library of Elastic Tensors for Lowermost Mantle Seismic Anisotropy Studies and Comparison With Seismic Observations. Geochemistry, Geophysics, Geosystems, 21(4). https://doi.org/10.1029/2019GC008883 Deng, J., Long, M. D., Creasy, N., Wagner, L., Beck, S., Zandt, G., Tavera, H., & Minaya, E. (2017). Lowermost mantle anisotropy near the eastern edge of the Pacific LLSVP: constraints from SKS–SKKS splitting intensity measurements. Geophysical Journal International, 210(2), 774–786. https://doi.org/10.1093/GJI/GGX190 55 Domeier, M., Doubrovine, P. V., Torsvik, T. H., Spakman, W., & Bull, A. L. (2016). Global correlation of lower mantle structure and past subduction. Geophysical Research Letters, 43(10), 4945–4953. https://doi.org/10.1002/2016GL068827 Faccenda, M., & Capitanio, F. A. (2012). Development of mantle seismic anisotropy during subduction-induced 3-D flow. Geophysical Research Letters, 39(11), n/a-n/a. https://doi.org/10.1029/2012GL051988 Forsyth, D., & Uyedaf, S. (1975). On the Relative Importance of the Driving Forces of Plate Motion. Geophysical Journal of the Royal Astronomical Society, 43(1), 163–200. https://doi.org/10.1111/j.1365-246X.1975.tb00631.x French, S. W., & Romanowicz, B. A. (2014). Whole-mantle radially anisotropic shear velocity structure from spectral-element waveform tomography. Geophysical Journal International, 199(3), 1303–1327. https://doi.org/10.1093/gji/ggu334 Fukao, Y., & Obayashi, M. (2013). Subducted slabs stagnant above, penetrating through, and trapped below the 660 km discontinuity. Journal of Geophysical Research: Solid Earth, 118(11), 5920–5938. https://doi.org/10.1002/2013JB010466 Fukao, Y., Widiyantoro, S., & Obayashi, M. (2001). Stagnant slabs in the upper and lower mantle transition region. Reviews of Geophysics, 39(3), 291–323. https://doi.org/10.1029/1999RG000068 Garnero, E. J., & Lay, T. (2003). D″ shear velocity heterogeneity, anisotropy and discontinuity structure beneath the Caribbean and Central America. Physics of the Earth and Planetary Interiors, 140(1–3), 219–242. https://doi.org/10.1016/J.PEPI.2003.07.014 Goes, S., Agrusta, R., van Hunen, J., & Garel, F. (2017). Subduction-transition zone interaction: A review. Geosphere, 13(3), 644–664. https://doi.org/10.1130/GES01476.1 Grand, S. P. (2002). Mantle shear–wave tomography and the fate of subducted slabs. Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 360(1800), 2475–2491. https://doi.org/10.1098/rsta.2002.1077 Grocholski, B., Catalli, K., Shim, S. H., & Prakapenka, V. (2012). Mineralogical effects on the detectability of the postperovskite boundary. Proceedings of the National academy of Sciences, 109(7), 2275-2279. Guillaume, B., Funiciello, F., & Faccenna, C. (2021). Interplays Between Mantle Flow and Slab Pull at Subduction Zones in 3D. Journal of Geophysical Research: Solid Earth, 126(5), e2020JB021574. https://doi.org/10.1029/2020JB021574 He, X., & Long, M. D. (2011). Lowermost mantle anisotropy beneath the northwestern Pacific: Evidence from PcS, ScS, SKS, and SKKS phases. Geochemistry, Geophysics, Geosystems, 12(12), 12012. https://doi.org/10.1029/2011GC003779 56 Hosseini, K., Sigloch, K., Tsekhmistrenko, M., Zaheri, A., Nissen-Meyer, T., & Igel, H. (2020). Global mantle structure from multifrequency tomography using P, PP and P-diffracted waves. Geophysical Journal International, 220(1), 96–141. https://doi.org/10.1093/gji/ggz394 Huang, J., & Davies, G. F. (2007). Stirring in three-dimensional mantle convection models and implications for geochemistry: 2. Heavy tracers. Geochemistry, Geophysics, Geosystems, 8(7). https://doi.org/10.1029/2007GC001621 Hutko, A. R., Lay, T., Garnero, E. J., & Revenaugh, J. (2006). Seismic detection of folded, subducted lithosphere at the core-mantle boundary. Nature, 441(7091), 333–336. https://doi.org/10.1038/nature04757 Jackson, J. M., & Thomas, C. (2021). Seismic and Mineral Physics Constraints on the D″ Layer. Mantle Convection and Surface Expressions, 193–227. https://doi.org/10.1002/9781119528609.CH8 Kaneshima, S., & Helffrich, G. (1999). Dipping low-velocity layer in the mid-lower mantle: Evidence for geochemical heterogeneity. Science, 283(5409), 1888–1891. https://doi.org/10.1126/SCIENCE.283.5409.1888 Karato, S. I., Jung, H., Katayama, I., & Skemer, P. (2008). Geodynamic Significance of Seismic Anisotropy of the Upper Mantle: New Insights from Laboratory Studies. Https://Doi.Org/10.1146/Annurev.Earth.36.031207.124120, 36, 59–95. https://doi.org/10.1146/ANNUREV.EARTH.36.031207.124120 Kito, T., Thomas, C., Rietbrock, A., Garnero, E. J., Nippress, S. E. J., & Heath, A. E. (2008). Seismic evidence for a sharp lithospheric base persisting to the lowermost mantle beneath the Caribbean. Geophysical Journal International, 174(3), 1019–1028. https://doi.org/10.1111/J.1365-246X.2008.03880.X/3/174-3-1019-FIG007.JPEG Koelemeijer, P., Ritsema, J., Deuss, A., & van Heijst, H. J. (2016). SP12RTS: a degree-12 model of shear- and compressional-wave velocity for Earth’s mantle. Geophysical Journal International, 204(2), 1024–1039. https://doi.org/10.1093/GJI/GGV481 Kustowski, B., Ekström, G., & Dziewoński, A. M. (2008). Anisotropic shear-wave velocity structure of the Earth’s mantle: A global model. Journal of Geophysical Research: Solid Earth, 113(B6). https://doi.org/10.1029/2007JB005169 Li, C., Van Der Hilst, R. D., Engdahl, E. R., & Burdick, S. (2008). A new global model for P wave speed variations in Earth’s mantle. Geochemistry, Geophysics, Geosystems, 9(5), 1213. https://doi.org/10.1029/2007GC001806 Li, M., & McNamara, A. K. (2013). The difficulty for subducted oceanic crust to accumulate at the Earth’s core-mantle boundary. Journal of Geophysical Research: Solid Earth, 118(4), 1807–1816. https://doi.org/10.1002/JGRB.50156 Li, M., McNamara, A. K., & Garnero, E. J. (2014). Chemical complexity of hotspots caused by 57 cycling oceanic crust through mantle reservoirs. Nature Geoscience, 7(5), 366–370. https://doi.org/10.1038/ngeo2120 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- Science Reviews, 196, 102874. https://doi.org/https://doi.org/10.1016/j.earscirev.2019.05.018 Lin, Y. P., Zhao, L., & Hung, S. H. (2014). Full-wave effects on shear wave splitting. Geophysical Research Letters, 41(3), 799–804. https://doi.org/10.1002/2013GL058742 Long, M. D. (2009). Complex anisotropy in D″ beneath the eastern Pacific from SKS-SKKS splitting discrepancies. Earth and Planetary Science Letters, 283(1–4), 181–189. https://doi.org/10.1016/j.epsl.2009.04.019 Long, M. D. (2013). Constraints on subduction geodynamics from seismic anisotropy. Reviews of Geophysics, 51(1), 76–112. https://doi.org/10.1002/rog.20008 Long, M. D., Xiao, X., Jiang, Z., Evans, B., & Karato, S. ichiro. (2006). Lattice preferred orientation in deformed polycrystalline (Mg,Fe)O and implications for seismic anisotropy in D″. Physics of the Earth and Planetary Interiors, 156(1–2), 75–88. https://doi.org/10.1016/J.PEPI.2006.02.006 Loubet, N., Ribe, N. M., & Gamblin, Y. (2009). Deformation modes of subducted lithosphere at the core-mantle boundary: An experimental investigation. Geochemistry, Geophysics, Geosystems, 10(10). https://doi.org/10.1029/2009GC002492 Lynner, C., & Long, M. D. (2014). Lowermost mantle anisotropy and deformation along the boundary of the African LLSVP. Geophysical Research Letters, 41(10), 3447–3454. https://doi.org/10.1002/2014GL059875 Mao, W., & Zhong, S. (2018). Slab stagnation due to a reduced viscosity layer beneath the mantle transition zone. Nature Geoscience, 11(11), 876–881. https://doi.org/10.1038/s41561-018-0225-2 McNamara, A. K., Karato, S. I., & Van Keken, P. E. (2001). Localization of dislocation creep in the lower mantle: implications for the origin of seismic anisotropy. Earth and Planetary Science Letters, 191(1–2), 85–99. https://doi.org/10.1016/S0012-821X(01)00405-8 McNamara, A. K., van Keken, P. E., & Karato, S.-I. (2003). Development of finite strain in the convecting lower mantle and its implications for seismic anisotropy. Journal of Geophysical Research: Solid Earth, 108(B5), 2230. https://doi.org/10.1029/2002JB001970 Merkel, S., McNamara, A. K., Kubo, A., Speziale, S., Miyagi, L., Meng, Y., Duffy, T. S., & Wenk, H. R. (2007). Deformation of (Mg,Fe)SiO3 post-perovskite and D″ anisotropy. Science, 316(5832), 1729–1732. https://doi.org/10.1126/SCIENCE.1140609 Miyagi, L., Kanitpanyacharoen, W., Kaercher, P., Lee, K. K. M., & Wenk, H. R. (2010). Slip 58 systems in MgSiO3 post-perovskite: Implications for D″ anisotropy. Science, 329(5999), 1639–1641. https://doi.org/10.1126/SCIENCE.1192465 Morra, G., Chatelain, P., Tackley, P., & Koumoutsakos, P. (2009). Earth curvature effects on subduction morphology: Modeling subduction in a spherical setting. Acta Geotechnica, 4(2), 95–105. https://doi.org/10.1007/s11440-008-0060-5 Murakami, M., Hirose, K., Kawamura, K., Sata, N., & Ohishi, Y. (2004). Post-Perovskite Phase Transition in MgSiO3. Science, 304(5672), 855–858. https://doi.org/10.1126/SCIENCE.1095932 Niu, F., Kawakatsu, H., & Fukao, Y. (2003). Seismic evidence for a chemical heterogeneity in the midmantle: A strong and slightly dipping seismic reflector beneath the Mariana subduction zone. Journal of Geophysical Research: Solid Earth, 108(B9), 2419. https://doi.org/10.1029/2002jb002384 Nowacki, A. (2013). Deformation of the Lowermost Mantle from Seismic Anisotropy. In Plate Deformation from Cradle to Grave: Seismic Anisotropy and Deformation at Mid-Ocean Ridges and in the Lowermost Mantle (pp. 99–122). https://doi.org/10.1007/978-3-642- 34842-6_4 Panning, M., & Romanowicz, B. (2006). A three-dimensional radially anisotropic model of shear velocity in the whole mantle. Geophysical Journal International, 167(1), 361–379. https://doi.org/10.1111/J.1365-246X.2006.03100.X/3/167-1-361-FIG017.GIF Ribe, N. M., Stutzmann, E., Ren, Y., & van der Hilst, R. (2007). Buckling instabilities of subducted lithosphere beneath the transition zone. Earth and Planetary Science Letters, 254(1–2), 173–179. https://doi.org/10.1016/j.epsl.2006.11.028 Ritsema, J., McNamara, A. K., & Bull, A. L. (2007). Tomographic filtering of geodynamic models: Implications for models interpretation and large-scale mantle structure. Journal of Geophysical Research: Solid Earth, 112(1). https://doi.org/10.1029/2006JB004566 Ritsema, J., van Heijst, H. J., & Woodhouse, J. H. (2004). Global transition zone tomography. Journal of Geophysical Research: Solid Earth, 109(B2). https://doi.org/10.1029/2003JB002610 Romanowicz, B., & Wenk, H. R. (2017). Anisotropy in the deep Earth. Physics of the Earth and Planetary Interiors, 269, 58–90. https://doi.org/10.1016/J.PEPI.2017.05.005 Rost, S., Garnero, E. J., & Williams, Q. (2008). Seismic array detection of subducted oceanic crust in the lower mantle. Journal of Geophysical Research: Solid Earth, 113(6). https://doi.org/10.1029/2007JB005263 Schellart, W. P. (2008). Kinematics and flow patterns in deep mantle and upper mantle subduction models: Influence of the mantle depth and slab to mantle viscosity ratio. Geochemistry, Geophysics, Geosystems, 9(3), n/a-n/a. https://doi.org/10.1029/2007GC001656 59 Shephard, G. E., Matthews, K. J., Hosseini, K., & Domeier, M. (2017). On the consistency of seismically imaged lower mantle slabs. Scientific Reports, 7(1), 1–17. https://doi.org/10.1038/s41598-017-11039-w Simmons, N. A., Myers, S. C., Johannesson, G., & Matzel, E. (2012). LLNL-G3Dv3: Global P wave tomography model for improved regional and teleseismic travel time prediction. Journal of Geophysical Research: Solid Earth, 117(B10), 10302. https://doi.org/10.1029/2012JB009525 Tackley, P. J. (2011). Living dead slabs in 3-D: The dynamics of compositionally-stratified slabs entering a “slab graveyard” above the core-mantle boundary. Physics of the Earth and Planetary Interiors, 188(3–4), 150–162. https://doi.org/10.1016/j.pepi.2011.04.013 Tan, E., Gurnis, M., & Han, L. (2002). Slabs in the lower mantle and their modulation of plume formation. Geochemistry, Geophysics, Geosystems, 3(11), 1–24. https://doi.org/10.1029/2001gc000238 Tsuchiya, T., Tsuchiya, J., Umemoto, K., & Wentzcovitch, R. M. (2004). Phase transition in MgSiO3 perovskite in the earth’s lower mantle. Earth and Planetary Science Letters, 224(3–4), 241–248. https://doi.org/10.1016/J.EPSL.2004.05.017 Van Der Hilst, R. D., Widiyantoro, S., & Engdahl, E. R. (1997). Evidence for deep mantle circulation from global tomography. Nature, 386(6625), 578–584. https://doi.org/10.1038/386578a0 Walpole, J., Wookey, J., Kendall, J. M., & Masters, T. G. (2017). Seismic anisotropy and mantle flow below subducting slabs. Earth and Planetary Science Letters, 465, 155–167. https://doi.org/10.1016/J.EPSL.2017.02.023 Wenk, H. R., Cottaar, S., Tomé, C. N., McNamara, A., & Romanowicz, B. (2011). Deformation in the lowermost mantle: From polycrystal plasticity to seismic anisotropy. Earth and Planetary Science Letters, 306(1–2), 33–45. https://doi.org/10.1016/j.epsl.2011.03.021 Wenk, H. R., Speziale, S., McNamara, A. K., & Garnero, E. J. (2006). Modeling lower mantle anisotropy development in a subducting slab. Earth and Planetary Science Letters, 245(1– 2), 302–314. https://doi.org/10.1016/j.epsl.2006.02.028 Wolf, J., Creasy, N., Pisconti, A., Long, M. D., & Thomas, C. (2019). An investigation of seismic anisotropy in the lowermost mantle beneath Iceland. Geophysical Journal International, 219, S152–S166. https://doi.org/10.1093/gji/ggz312 Wu, X., Lin, J.-F., Kaercher, P., Mao, Z., Liu, J., Wenk, H.-R., & Prakapenka, V. B. (2017). Seismic anisotropy of the D″ layer induced by (001) deformation of post-perovskite. Nature Communications, 8(1), 14669. https://doi.org/10.1038/ncomms14669 Zahirovic, S., Müller, R. D., Seton, M., & Flament, N. (2015). Tectonic speed limits from plate kinematic reconstructions. Earth and Planetary Science Letters, 418, 40–52. https://doi.org/10.1016/j.epsl.2015.02.037 60 Zhao, D. (2004). Global tomographic images of mantle plumes and subducting slabs: Insight into deep Earth dynamics. Physics of the Earth and Planetary Interiors, 146(1–2), 3–34. https://doi.org/10.1016/j.pepi.2003.07.032 Zhong, S. (2006). Constraints on thermochemical convection of the mantle from plume heat flux, plume excess temperature, and upper mantle temperature. Journal of Geophysical Research: Solid Earth, 111(4), 4409. https://doi.org/10.1029/2005JB003972 61 APPENDIX Figure 3.8 Zoom-in snapshots of temperature cross-sections at 0⁰ longitude of (a) Ra = 5.0 x 107 and (b) Ra = 1.0 x 106 while other parameters are the same as in Case 13. Cold regions (blue) are flat or horizontal in (a) and are seemingly thickened when large blob of cold droplets drops from the plate before they reach the trench. 62 Figure 3.9 Time average of asymmetry character (α) and azimuthal asymmetry character (β) of lowermost mantle lateral flow patterns within slabs of 1,112 km (10 degree arc angle), 3,336 km (30 degree arc angle) and 7,784 km (70 degree arc angle) as a function of slab strength (parameter A). (a) and (b), time average of α and β as a function of trench length and slab strength are color- coded in (a) and (b), respectively. Blue means low value and pink means high value. Vertical axes are A values and horizontal axes are trench length in arc angle unit. (c) and (d), time average of α and β as a function of slab strength, with orange markers and lines representing slabs with a trench length of 10 degree arc angle, green symbols representing slabs with a trench length of 30 degree arc angle, and yellow symbols representing slabs with a trench length of 70 degree arc angle. Models with A of 4.61, 9.21, and 13.82 lead to 10x, 100x, 1,000x, and 10,000x viscosity contrast between the surface plate and the CMB in the models. Rhombuses are data points; solid lines connect the data points. 63 Figure 3.10 Time average of asymmetry character (α) and azimuthal asymmetry character (β) of lowermost mantle lateral flow patterns within weak slabs (with A of 4.60517), moderate strength slabs (with A of 9.21034) and strong slabs (with a of 13.815511) as a function of trench length (in arc angle). (a) and (b), time average of α and β as a function of trench length and slab strength are color-coded in (a) and (b), respectively. Blue means low value and pink means high value. Vertical axes are A values and horizontal axes are trench length in arc angle unit. (c) and (d), time average of α and β as a function of trench length, with blue markers and lines representing weak slabs, orange symbols representing moderate strength slabs, and grey symbols representing strong slabs. Models with A of 4.61, 9.21, and 13.82 lead to 10x, 100x, 1,000x, and 10,000x viscosity contrast between the surface plate and the CMB in the models. Squares are data points; solid lines connect the data points. 64 Figure 3.11 Representative snapshots for a wide slab (7,778 km of trench length) with moderate strength (100x viscosity contrast between the surface plate and the CMB) at t = 3.15 billion years after it reaches thermal equilibrium (Case 10). Similar to Case 2 (Figure 5 a-c), in the lowermost mantle, slabs frequently deform into seemingly separate cold patches. a, residual temperature cross-section of the slab along the slab axis, with the blue contour being the colder residual temperature contour showing the slab. Colder regions are blue and warmer regions are red. b, 3D residual temperature isosurface showing the slab geometry. c, map view of radial velocity at 320 km depth, with positive meaning the material is moving upward and negative meaning the material is moving downward. The radial velocity is represented by the color scale beneath. Black arrows are lateral velocity. d, map view of the magnitude of tangential velocity at 320 km depth, with red meaning higher tangential velocity and blue meaning lower tangential velocity. e,f, map view of residual temperature at 60 km depth at t = 357 Ma and t = 3.15 billion years after it reaches thermal equilibrium, respectively. The color scale in (e) and (f) are the same with (a). The orange contour is radial velocity contour showing the slab geometry. (a) images that there are two cylindrical colder cores in the lowermost mantle slab. Here in (b) we can clearly track the cold patches all the way up to the surface. In c and d, we can see that the cylindrical structure in the back of the slab has a higher radial velocity going downward and higher tangential velocity near the surface. 65 Figure 3.12 A representative snapshot of a narrow slab with a trench length of 445 km (Case 14) producing shearing in the lowermost mantle driven by the surrounding flow. Colder and hotter temperatures are represented by blue and red colors, respectively. zoom-in figure on the bottom right corner shows the map view of the lowermost mantle lateral flow velocity (black arrows), with a background of temperature field. Blue isosurface is the slab that induces the flow pattern. 66 CHAPTER 4: ULVZ CROSS-SECTIONAL SHAPE GIVES CLUE TO THE ORIGIN AND INTRINSIC RHEOLOGICAL NATURE OF THEMSELVES AND LLSVPS 4.1 Abstract Small-scale ultralow velocity zones (ULVZs) have been detected as disparate patches on the core-mantle boundary, oftentimes at the margins of or within the two large low-shear-velocity provinces (LLSVPs) in the Earth’s lowermost mantle. LLSVPs and ULVZs are hypothesized to be caused by intrinsically denser large-scale primordial compositional reservoirs (thermochemical piles). Because the compositions of the LLSVPs and ULVZs cannot be determined directly, the origin and the rheological properties of the reservoirs remain elusive. The piles and ultra-dense materials may have different intrinsic viscosity from the background mantle due to their distinct composition (e.g., Ballmer et al., 2017), and the intrinsic rheological property change of ultra- dense materials could have a dynamical feedback to the shape of ultra-dense materials (e.g., M. Li et al., 2017; McNamara et al., 2010). With more advanced seismic techniques, it is hopeful to resolve the cross-sectional shape of ULVZs in the future. However, it is unclear how the cross- sectional shape of ultra-dense materials could be related to the intrinsic rheological properties of large compositional reservoirs. Here we perform high-resolution geodynamical calculations to investigate if and how the cross-sectional shapes of ULVZs are related to the rheology of the large- scale lower mantle compositional reservoirs. We find that ultra-dense materials accumulate into symmetrical patches within the interior of thermochemical piles and accumulate into asymmetrical patches along the margin of thermochemical piles. Furthermore, the type of viscous coupling between piles and the background mantle controls the symmetry of cross-sectional shapes of ultra- dense materials accumulated at pile edges. Large compositional reservoirs with higher intrinsic viscosity than the background mantle decrease the differential viscous coupling between the pile 67 and the background mantle. Therefore, the cross-sectional shape of ultra-dense materials becomes more symmetric or exhibits opposite asymmetry at the margins of more viscous piles. Thus, the lateral morphology of the patches provide insight into the type of viscous coupling between thermochemical piles and the background mantle, and furthermore, the origin and the rheological properties of the piles. With high-resolution seismic imaging technique of ULVZ cross-sectional shapes (e.g., Z. Li et al., 2022), the shapes of ULVZs can therefore help to constrain the rheological nature and the origin of the large-scale lower mantle compositional reservoirs. 4.2 Introduction The two large low-shear-velocity provinces (LLSVPs) in the Earth’s lowermost mantle beneath the Pacific and Africa are hypothesized to be formed by the large accumulation of primordial remnants of Earth’s early differentiation process (e.g., Kellogg et al., 1999; Lee et al., 2010). Travel-time and waveform studies infer sharp seismic-wave velocity gradients along LLSVP boundaries (e.g., Bréger & Romanowicz, 1998; Ford et al., 2006; Ni & Helmberger, 2003b; X. Sun et al., 2007; To et al., 2005; Van Keken, 1997). Normal mode studies (e.g., Davaille et al., 2002; Ishii & Tromp, 2001, 2004) and body tide modeling (e.g., Lau et al., 2017) suggest that LLSVPs may be anomalous in density. Furthermore, geochemical observations from the Earth's surface imply a more primordial compositional reservoir in the Earth’s lower mantle (e.g., Tackley, 2000). The origin of LLSVPs may play a critical role in understanding the large-scale dynamics in the Earth’s mantle and, thus, the evolution of the Earth. The observations, though, do not exclude other hypotheses (e.g., Christensen & Hofmann, 1994) (i.e., LLSVPs are caused by the accumulation of subducted oceanic crust (e.g., Ballmer et al., 2017; C. Huang et al., 2020; M. Li et al., 2014; Nakagawa et al., 2010); super plume (e.g., D. R. Davies, 2015); a cluster of thermal plumes (e.g., Davaille & Romanowicz, 2020); thermochemical plumes (e.g., Thompson & 68 Tackley, 1999)) for the origin of LLSVPs. Furthermore, the observations provide little clue to the rheological nature of LLSVPs. At a significantly smaller scale than LLSVPs, ultralow velocity zones (ULVZs) have been mapped as thin (5–40 km) and small-scale (100–1000 km) (e.g., Hutko et al., 2009; Rost et al., 2005; Thorne et al., 2013) patches often located near or within the LLSVPs (e.g., McNamara et al., 2010) on the core-mantle boundary (CMB). Recent seismic studies (e.g., Thorne et al., 2021; Yu & Garnero, 2018) have proposed distribution maps for ULVZs on which some possible ULVZs are located outside of the LLSVPs but we focus on those associated with LLSVPs. ULVZs have a substantial seismic wave velocity reduction of 10% to 30% (e.g., Rost et al., 2005; Yu & Garnero, 2018), which is much slower than LLSVPs. ULVZs are hypothesized to be caused by solely partial melt (e.g., Y. He & Wen, 2009) or compositional heterogeneity (i.e., iron-enriched (Mg,Fe)O (e.g., Wicks et al., 2010, 2017), iron-enriched post-perovskite (e.g., W. L. Mao et al., 2006), subducted banded iron formation (e.g., Dobson & Brodholt, 2005), etc.) which may or may not contain partial melting. Although ULVZs located within LLSVPs may dominantly be caused by partial melting (e.g., M. Li et al., 2017), most ULVZ materials accumulated at the edge of compositional reservoirs must surpass a 5% density increase than the background mantle to survive entrainment from the background mantle (e.g., M. Li et al., 2017; McNamara et al., 2010). It is found in previous geodynamics studies (e.g., M. Li et al., 2017; McNamara et al., 2010) that ULVZ morphologies could contain information about the dynamics of compositional reservoirs. ULVZs at the compositional reservoir boundaries have an asymmetrical cross-sectional shape (e.g., M. Li et al., 2017; McNamara et al., 2010) that is consistent with the Hawaiian ULVZ morphology observed in a recent seismic study (e.g., Z. Li et al., 2022). They are found to be thicker on the outside of the thermochemical pile due to differential viscous coupling between the 69 background mantle and compositional reservoirs (e.g., M. Li et al., 2017; McNamara et al., 2010). The previous studies (e.g., M. Li et al., 2017; McNamara et al., 2010) assume that piles and the background mantle have the same intrinsic rheology nature, which, however, may not be true as they are hypothesized to have different compositions. For instance, if pile materials have a larger grain size than the background mantle, they will have a higher intrinsic viscosity (e.g., Karato & Wu, 1993; McNamara et al., 2003) (Table 4.1). Previous studies (e.g., Davaille et al., 2002, 2003; McNamara & Zhong, 2004) showed that an addition in intrinsic viscosity of denser materials could lead to an increase in the amount of thermochemical “domes” (superplumes) and the formation of a high viscosity rind along the interface between the dense component and the background mantle. Although these studies did not implement ULVZ material, the possible change of the dynamics of the thermochemical piles due to the intrinsic viscosity change could influence how the ultra-dense material interacts with piles and the background mantle. It is unclear, therefore, whether we will always expect asymmetrical accumulations of ultra-dense materials at thermochemical pile edges and how the viscous coupling along the pile boundaries controls ULVZ morphologies and symmetry, in particular. Combined with more high-resolution seismic observations on ULVZ morphologies (e.g., Cottaar & Romanowicz, 2012; Z. Li et al., 2022; Thorne et al., 2013; Yu & Garnero, 2018; K. Yuan & Romanowicz, 2017) in the future, understanding how ULVZ morphologies are controlled by the viscous coupling along compositional reservoir boundaries provide meaningful insights into the rheological nature and origin of LLSVPs. Thus, we perform high-resolution 2.5 D annulus thermochemical numerical calculation and hypothesize that ULVZs are caused by ultra-dense material accumulations (also referred to as “UDMAs” thereafter). We study how UDMA shapes, especially in terms of symmetry, respond to thermochemical piles with different intrinsic rheological properties. 70 4.3 Method We modified the Spherical convection code CitcomCU (e.g., Zhong, 2006) to include multiple species and compositional rheology to solve the conservation equations of mass, momentum and energy in Boussinesq approximation. We employ a Rayleigh number Ra = 1 × 107 for all the cases. The viscosity is depth-, temperature-, and composition-dependent: 𝜂(𝑇, 𝑧, 𝐶) = 𝜂𝑧(𝑧)𝜂𝑐(𝐶)exp [𝐴(0.35 − 𝑇)]…………………………………………Equation (4) Here ηz(z) = 1 and 30 for the upper and lower mantle, leading to a 30x viscosity increase from the upper mantle to the lower mantle. ηC(C) is a viscosity coefficient depending on a given compositional species. A and T are non-dimensional activation coefficient and temperature, respectively. We use A = 9.21, which leads to a viscosity range of 104 due to temperature differences between the coldest and hottest regions of the model. We employ an annulus geometry in a cylindrical coordinate (r, ϕ), where the dimensionless radius ranges from 0.55 to 1.0 (thus, corresponds to the CMB and the Earth’s surface). The domain is divided into 2560 elements in the ϕ direction and 192 elements in the radial direction. All boundaries are free-slip, isothermal on the top and bottom and insulating on both sides. The models are heated from below and internally with non-dimensional heat production of 60. To establish an initial condition, we first perform a lower resolution (1,600 × 160 elements) calculation with two compositional components (thermochemical pile material and background mantle) until the model reaches quasi-steady thermal state and thermochemical piles are developed. Pile material has a buoyancy number of B = 0.8 (about 1.5-4.0% denser than the background mantle). About 24.5 million tracers are introduced to model the compositional field using ratio tracer method (e.g., Tackley & King, 2003). Thereafter, we interpolate the temperature and composition field to higher resolution (2560 × 192 elements) and introduce a uniform layer 71 composed of ultra-dense ULVZ material with a buoyancy number of B = 2.0 (about 5.0-10.0% intrinsically denser than the background mantle) in the lowermost 7.6 km of the mantle. The new mesh is refined in the lowermost 38 km of the mantle, providing a resolution of 2.5 km radially and ~8.5 km laterally near the CMB. We define a reference case (Case 1) in which the viscosity of thermochemical piles and the background mantle are not dependent on composition with ηC(C) = 1, meaning that they have the same intrinsic viscosity. The parameters used in the reference case are shown in Table 4.2. Figure 4.1 shows a representative time snapshot of the reference case at 1,135 million years after the start of the calculation. 4.4 Results Figure 4.1a is the temperature field of the upper semi-sphere of the model. Figure 4.1b is the compositional field of the same semi-sphere in Figure 4.1a. At the time, piles have formed and are framed by downwellings. Plumes are rising up on top of the piles, entraining a small amount of pile materials. 5 discontinuous patches of UDMAs (dark red region in Figure 4.1b) have accumulated in the lowermost 50 km of the mantle: 4 at pile edges, and one within the pile interior. The five boxed UDMA regions in Figure 4.1b are enlarged and displayed in Figure 4.1c-g. Figures 4.1c-g illustrate the representative cross-sectional shape of the UDMAs (yellow contour) and viscosity field surrounding the accumulations. All the accumulations displayed are in convergent flow regions and are in a quasi-stable state. Figures 4.1c, 4.1d, 4.1e, and 4.1g demonstrate that ultra-dense materials typically accumulate into asymmetrical thin patches along the boundaries of piles, thicker on the side outboard of the pile, 72 Figure 4.1 ULVZ cross-sectional shapes in the reference case, shown at 1,135 million years after the start of the calculation. a Temperature field (non-dimensional) of one of the semi spheres of the model. b Compositional field of cross-section in a. Pink denotes thermochemical pile materials, black shows the background mantle, and red regions are UDMAs. The pile contains 5 UDMAs, zoomed-in in Figure 1c-g. c-g Zoom-ins of the regions around UDMAs in b, from left to right. Color represents the logarithm of non-dimensional viscosity. Gray lines are temperature contours with an interval of non-dimensional value 0.1. Yellow lines outline UDMAs. Red lines delineate pile material. which is consistent with previous studies (e.g., M. Li et al., 2017; McNamara et al., 2010). The accumulation produced at pile edges can survive and maintain its shape for as long as 1,400 million years. Although the lateral width of the accumulations at pile margins varies from place to place, and from time to time (Figure 4.3 and Figure 4.5), from ~80 km to ~400 km, the asymmetry feature of the cross-sectional shape of the accumulations remains. Along the boundary of the pile (red contour), the viscosity in the background mantle is nearly 4x higher than the inside of the reservoir due to the temperature-dependent rheology. Considering that viscosity dominantly controls stress over strain rate at pile edges, the stress acting on the ULVZ material from the side in contact with 73 the background mantle is thus greater than that from the side within the primordial reservoir. Therefore, the differential stress acting on the two sides of ultra-dense material shapes the accumulation into asymmetrical patches. Figure 4.1f displays that the UDMAs are symmetrical thin patches within the pile interior. Within the interior of thermochemical reservoirs, the accumulations are quasi-stable under upwelling regions within piles for as long as 178 million years but will eventually be transported to pile boundaries by pile internal flow. The viscosity is uniform, surrounding the accumulation within the pile interior, leading to a uniform viscous coupling between UDMAs and their surroundings. Case 2 is identical to the reference case except that the intrinsic viscosity of compositional reservoirs is 100x greater (ηC(C) = 100) than that of the background mantle. Figure 4.2 shows a representative time snapshot of Case 2 at 580 million years, demonstrating the cross-sectional shapes of UDMAs. ULVZ material accumulates at the edge of the reservoir (Figure 4.2c and 4.2f) and within the interior of the reservoir (Figure 4.2d and 4.2e), forming UDMAs in a quasi- equilibrium state. The pile on the right side has just been separated from the pile on its left side by the downwelling on the left and is being swept rightwards. Interestingly, in this case, a high viscosity rind also forms along the interface between the compositional reservoir and the background mantle (Figure 4.10), which is consistent with previous geodynamic models (e.g., Davaille et al., 2003). Similar to Case 1, the viscosity is uniform, surrounding the accumulation within the interior of the reservoir. Consistent with Case 1, the accumulation within the pile interior has a symmetrical cross-sectional shape (Figure 4.2d and 4.2e). At the edge of the pile (red contour), in a convergent flow region, UDMAs are longer-lived and accumulate into a more symmetrical patch than that in Case 1 (Figure 4.2c and 4.2f). The UDMAs even slightly exhibit an opposite asymmetry, with the side in contact with the background mantle ~10% thinner than the 74 Figure 4.2 UDMAs cross-sectional shapes in the intrinsically more viscous pile (Case 2), shown at 580 million years after the start of the calculation. The case is identical to the reference case except that the piles are intrinsically 100x more viscous than the background mantle (ηC(C) = 100). Panels represent the same fields and colors as in Figure 4.1. side inboard of the pile. The size of UDMAs at pile edges are of variety, with lateral lengths from ~100 km to ~300 km, but the symmetry to slightly opposite asymmetry characteristic of its cross- sectional shape is consistent. Surrounding the UDMAs (yellow contour), the side within the reservoir is ~6x more viscous than the side in contact with the background mantle. Compared to Case 1, the opposite differential viscous coupling on the two sides of UDMAs causes the stress force acting on the accumulation from the pile interior to be larger than that from the background mantle. Thus, the opposite differential stress shapes the accumulation that exhibits an opposite asymmetry. 75 Figure 4.3 Summary chart of characteristic cross-sectional shapes of UDMA (red, yellow, cyan and blue contours) at pile edges (orange-shaded column) and within pile interiors (green-shaded column), in different pile intrinsic rheology regimes (Case 1-4). The UDMA contours are compositional field contours at C=1.0, zoomed-in at the same scale, and representing the cross- sectional shapes of UDMAs. The UDMA shapes of each case are randomly picked and representative from time to time, and place to place. The shapes are all oriented so that the left direction is the side in contact with the background mantle and the right direction is the side inboard of the reservoir. We examined other cases with 10x greater (ηC(C) = 10) and less (ηC(C) = 0.1) pile intrinsic viscosity (Case 3-4, Figures 4.6 and 4.7) and we found that symmetry of the cross-sectional shape of the accumulations formed along pile boundaries followed the same trend; as pile intrinsic viscosity increases, the accumulation shape changes from asymmetrical, thicker on the side in contact with the background mantle, to more symmetrical, and even oppositely asymmetrical, thinner on the side outboard of the reservoir. Figure 4.3 displays characteristic ULVZ shapes in different pile rheology regimes (Case 1-4) and demonstrates the trend that ULVZ shape symmetry 76 changes as pile intrinsic viscosity changes, as summarized above. In the intermediate case (Case 3) between Case 1 and Case 2, the differential viscous coupling of the side in contact with the background mantle to the side within the pile interior is ~ 1-3x. Thus, the differential stress is reduced, leading to a more symmetrical-shaped ULVZ material accumulation (Figure 4.3). In Case 4, the viscous coupling is ~ 10-30x. Thus, the differential stress is increased compared with Case 1, shaping the UDMAs into an even more asymmetrical and “arrow-like” shape (Figure 4.3). It is also noted that as pile intrinsic viscosity decreases, UDMAs produced at pile edges become wider and thinner (Case 4, Figure 4.3 and 4.7), which is consistent with the phenomena when only the intrinsic compositional viscosity of UDMAs is decreased (e.g., McNamara et al., 2010). For all the cases, the quasi-stable UDMAs produced within the pile interior are symmetrical. In addition, we explored the case (Case 5) in which the pile is intrinsically 1,000x more viscous (ηC(C) = 1,000) than the background mantle. We found that the compositional reservoir is more stable compared to Cases 1-4, giving a slower response to the change of downwelling patterns. Moreover, the convection pattern within piles is different from that in Cases 1-4. Plumes could form away from piles after the piles have developed (Figure 4.8). UDMAs could be dragged towards plume roots outside the piles, becoming flattened out and maintaining its position and shape for at least 800 million years (Figure 4.9). In addition, we found that the morphology and position of thermochemical piles could remain stable for more than 800 million years (Figure 4.8). Because the dynamic behavior of piles and plumes is different from the cases mainly discussed here, Case 5 is in a different dynamic regime from Case 1-4, and although the trends in UDMA shapes may not apply, it is beyond the scope of this study here. We also investigated other combinations of parameters such as Rayleigh number, the intrinsic density of compositional reservoirs, the intrinsic density of UDMAs, the intrinsic 77 viscosity of UDMAs, and the temperature dependence of viscosity. The varying parameters used in all the cases are listed in Table 4.3. Although varying these parameters leads to second-order differences in UDMA shapes from Case 1 (Figure 4.11-4.15), the fundamental conclusions about the symmetry of cross-sectional shapes of the UDMAs remained unchanged. 4.5 Discussion and Conclusion In this work, we assume that LLSVPs are caused by thermochemical piles and ULVZs have a compositional origin, which may or may not include partial melting. Such assumptions are motivated by the seismic observations (e.g., Hutko et al., 2009; McNamara et al., 2010; Rost et al., 2005; Thorne et al., 2013, 2021; Yu & Garnero, 2018) indicating that ULVZs could possibly be located outside of the LLSVPs. If ULVZs are solely caused by partial melting, they are unlikely to be located outside of the LLSVPs (e.g., McNamara et al., 2010), because their origin then rely on the higher temperature LLSVPs have than the background mantle. Furthermore, a recent study (e.g., Jenkins et al., 2021) suggested that a few resolved ULVZs are asymmetrical at the margin of LLSVPs. We investigated if cross-sectional shapes of ULVZ could provide insights or constraints on the rheological nature of LLSVPs. The Rheological nature of LLSVPs is hard to measure directly by seismic observations and lab experiments, while cross-sectional shapes of ULVZs are easier to observe with more seismic data coverage and better techniques. We find that the viscous coupling between the interface of the UDMAs and its surrounding material controls the cross-sectional shape of UDMAs at first order. Figure 4.4 summarizes and abstracts the typical UDMA shapes when UDMAs are in a quasi-equilibrium state. Differential viscous coupling at the edge of thermochemical piles (LLSVPs in the model) leads to differential stress acting on the two sides of UDMAs and causes asymmetrical UDMA shapes: uniform viscous coupling within the pile interior leads to uniform stress acting on the two sides of UDMAs and 78 Figure 4.4 Conceptual cartoon of possible cross-sectional shapes of ULVZs (red) in intrinsically less viscous thermochemical piles (green) and more viscous piles (orange). Ultra-dense materials (ULVZs) could accumulate within the interior and at the margins of thermochemical piles. A sinking slab is shown in blue, pushing the more primordial materials into two piles on the two sides. Rising plumes on top of the piles are shown in pink. Arrows indicate mantle flow directions. causes symmetrical UDMA shapes. Furthermore, if compositional reservoirs and the background mantle have the same rheology, UDMAs along pile edges are asymmetric (thicker outboard and thinner inboard of the pile); if reservoirs have a higher intrinsic viscosity than background mantle, UDMAs along pile edges become more symmetric or exhibit opposite asymmetry (thinner outboard and thicker inboard of the pile). UDMAs are typically found at pile edges and within the pile interior in convergent flow regions. There are some exceptions to the symmetry of UDMAs. It is found in our models that when UDMAs are in a transition mode, meaning that they are being transported to pile edges, they could either inherit the shape symmetry from when they were previously in a stable state (Figure 4.12, region 2) or transform into a very thin and long symmetrical cross-sectional shape, which is consistent with other geodynamic models (e.g., McNamara et al., 2010). In addition, when 79 UDMAs are in regions where two piles are merging or separating, their cross-sectional shape could be deformed by the downwelling flow sweeping the more primordial material and is symmetrical (e.g., Figure 4.13, region 2). But the above cases are temporary in the models and commonly do not last long. It is also found that when the intrinsic compositional viscosity of thermochemical piles is 1,000x more than that of the background mantle (equivalent to 16 x of grain size increase), the piles become much more stable compared to previous cases, and the mantle convection goes to a different dynamic regime, which is interesting to investigate in the future but beyond the scope of this work. In conclusion, our work shows that we could determine the rheological nature of LLSVPs by resolving the symmetry of ULVZ cross-sectional shapes. Although high-resolution ULVZ images are under-explored, the recent work (e.g., Jenkins et al., 2021; Z. Li et al., 2022) towards resolving ULVZ detailed structures cast light on it. With more high-resolution seismological observations on ULVZ structures, it is promising for us to understand the rheological nature of the LLSVPs and ULVZs. 80 REFERENCES Ballmer, M. D., Houser, C., Hernlund, J. W., Wentzcovitch, R. M., & Hirose, K. (2017). Persistence of strong silica-enriched domains in the Earth’s lower mantle. Nature Geoscience, 10(3), 236–240. https://doi.org/10.1038/ngeo2898 Bréger, L., & Romanowicz, B. (1998). Three-Dimensional Structure at the Base of the Mantle Beneath the Central Pacific. Science, 282(5389), 718–720. https://doi.org/10.1126/science.282.5389.718 Christensen, U. R., & Hofmann, A. W. (1994). Segregation of subducted oceanic crust in the convecting mantle. In JOURNAL OF GEOPHYSICAL RESEARCH (Vol. 99, Issue B10). https://doi.org/10.1029/93JB03403 Cottaar, S., & Romanowicz, B. (2012). An unsually large ULVZ at the base of the mantle near Hawaii. Earth and Planetary Science Letters, 355–356, 213–222. https://doi.org/10.1016/j.epsl.2012.09.005 Davaille, A., Girard, F., & Le Bars, M. (2002). How to anchor hotspots in a convecting mantle? Earth and Planetary Science Letters, 203(2), 621–634. https://doi.org/10.1016/S0012- 821X(02)00897-X Davaille, A., Le Bars, M., & Carbonne, C. (2003). Thermal convection in a heterogeneous mantle. Comptes Rendus - Geoscience, 335(1), 141–156. https://doi.org/10.1016/S1631- 0713(03)00003-8 Davaille, A., & Romanowicz, B. (2020). Deflating the LLSVPs: Bundles of Mantle Thermochemical Plumes Rather Than Thick Stagnant “Piles.” Tectonics, 39(10), e2020TC006265. https://doi.org/10.1029/2020TC006265 Davies, D. R. (2015). Thermally dominated deep mantle LLSVPs: A review Antilles (VoiLA) View project VoiLA: Volatiles Pathways in the Lesser Antilles View project. Springer, 441–477. https://doi.org/10.1007/978-3-319-15627-9_14 Dobson, D. P., & Brodholt, J. P. (2005). Subducted banded iron formations as a source of ultralow-velocity zones at the core-mantle boundary. Nature, 434(7031), 371–374. https://doi.org/10.1038/nature03430 Ford, S. R., Garnero, E. J., & McNamara, A. K. (2006). A strong lateral shear velocity gradient and anisotropy heterogeneity in the lowermost mantle beneath the southern Pacific. Journal of Geophysical Research: Solid Earth, 111(B3), n/a-n/a. https://doi.org/10.1029/2004JB003574 He, Y., & Wen, L. (2009). Structural features and shear‐velocity structure of the “Pacific Anomaly” J. Geophys. Res, 114(2), 2309. https://doi.org/10.1029/2008JB005814 Huang, C., Leng, W., & Wu, Z. (2020). The Continually Stable Subduction, Iron-Spin 81 Transition, and the Formation of LLSVPs From Subducted Oceanic Crust. Wiley Online Library, 125(1). https://doi.org/10.1029/2019JB018262 Hutko, A. R., Lay, T., & Revenaugh, J. (2009). Localized double-array stacking analysis of PcP: D″ and ULVZ structure beneath the Cocos plate, Mexico, central Pacific, and north Pacific. Physics of the Earth and Planetary Interiors, 173(1–2), 60–74. https://doi.org/10.1016/j.pepi.2008.11.003 Ishii, M., & Tromp, J. (2001). Even-degree lateral variations in the Earth’s mantle constrained by free oscillations and the free-air gravity anomaly. Geophysical Journal International, 145(1), 77–96. https://doi.org/10.1111/j.1365-246X.2001.00385.x Ishii, M., & Tromp, J. (2004). Constraining large-scale mantle heterogeneity using mantle and inner-core sensitive normal modes. Physics of the Earth and Planetary Interiors, 146(1–2), 113–124. https://doi.org/10.1016/j.pepi.2003.06.012 Jenkins, J., Mousavi, S., Li, Z., & Cottaar, S. (2021). A high-resolution map of Hawaiian ULVZ morphology from ScS phases. Earth and Planetary Science Letters, 563, 116885. https://doi.org/10.1016/j.epsl.2021.116885 Karato, S. I., & Wu, P. (1993). Rheology of the upper mantle: A synthesis. Science, 260(5109), 771–778. https://doi.org/10.1126/science.260.5109.771 Kellogg, L. H., Hager, B. H., & Van Der Hilst, R. D. (1999). Compositional stratification in the deep mantle. Science, 283(5409), 1881–1884. https://doi.org/10.1126/science.283.5409.1881 Lau, H. C. P., Mitrovica, J. X., Davis, J. L., Tromp, J., Yang, H. Y., & Al-Attar, D. (2017). Tidal tomography constrains Earth’s deep-mantle buoyancy. Nature, 551(7680), 321–326. https://doi.org/10.1038/nature24452 Lee, C. T. A., Luffi, P., Höink, T., Li, J., Dasgupta, R., & Hernlund, J. (2010). Upside-down differentiation and generation of a primordial lower mantle. Nature, 463(7283), 930–933. https://www.nature.com/articles/nature08824 Li, M., McNamara, A. K., & Garnero, E. J. (2014). Chemical complexity of hotspots caused by cycling oceanic crust through mantle reservoirs. Nature Geoscience, 7(5), 366–370. https://doi.org/10.1038/ngeo2120 Li, M., McNamara, A. K., Garnero, E. J., & Yu, S. (2017). Compositionally-distinct ultra-low velocity zones on Earth’s core-mantle boundary. Nature Communications, 8(1), 1–9. https://doi.org/10.1038/s41467-017-00219-x Li, Z., Leng, K., Jenkins, J., & Cottaar, S. (2022). Kilometer-scale structure on the core–mantle boundary near Hawaii. Nature Communications 2022 13:1, 13(1), 1–8. https://doi.org/10.1038/s41467-022-30502-5 Mao, W. L., Mao, H. K., Sturhahn, W., Zhao, J., Prakapenka, V. B., Meng, Y., Shu, J., Fei, Y., & 82 Hemley, R. J. (2006). Iron-rich post-perovskite and the origin of ultralow-velocity zones. Science, 312(5773), 564–565. https://doi.org/10.1126/science.1123442 McNamara, A. K., Garnero, E. J., & Rost, S. (2010). Tracking deep mantle reservoirs with ultra- low velocity zones. Earth and Planetary Science Letters, 299(1–2), 1–9. https://doi.org/10.1016/j.epsl.2010.07.042 McNamara, A. K., van Keken, P. E., & Karato, S.-I. (2003). Development of finite strain in the convecting lower mantle and its implications for seismic anisotropy. Journal of Geophysical Research: Solid Earth, 108(B5), 2230. https://doi.org/10.1029/2002JB001970 McNamara, A. K., & Zhong, S. (2004). Thermochemical structures within a spherical mantle: Superplumes or piles? Journal of Geophysical Research: Solid Earth, 109. https://doi.org/10.1029/2003JB002847 Nakagawa, T., Tackley, P. J., Deschamps, F., & Connolly, J. A. D. (2010). The influence of MORB and harzburgite composition on thermo-chemical mantle convection in a 3-D spherical shell with self-consistently calculated mineral physics. Earth and Planetary Science Letters, 296(3–4), 403–412. https://doi.org/10.1016/j.epsl.2010.05.026 Ni, S., & Helmberger, D. V. (2003). Ridge-like lower mantle structure beneath South Africa. Journal of Geophysical Research: Solid Earth, 108(B2). Rost, S., Garnero, E. J., Williams, Q., & Manga, M. (2005). Seismological constraints on a possible plume root at the core-mantle boundary. Nature, 435(7042), 666–669. https://www.nature.com/articles/nature03620 Sun, X., Song, X., Zheng, S., & Helmberger, D. V. (2007). Evidence for a chemical-thermal structure at base of mantle from sharp lateral P-wave variations beneath Central America. Proceedings of the National Academy of Sciences, 104(1), 26–30. https://doi.org/10.1073/pnas.0609143103 Tackley, P. J. (2000). Mantle convection and plate tectonics: Toward an integrated physical and chemical theory. In Science (Vol. 288, Issue 5473, pp. 2002–2007). American Association for the Advancement of Science. https://doi.org/10.1126/science.288.5473.2002 Tackley, P. J., & King, S. D. (2003). Testing the tracer ratio method for modeling active compositional fields in mantle convection simulations. Geochemistry, Geophysics, Geosystems, 4(4). https://doi.org/10.1029/2001GC000214 Thompson, P. E., & Tackley, P. J. (1999). Generation of mega‐plumes from the core‐mantle boundary in a compressible mantle with temperature‐dependent viscosity. In Wiley Online Library (Vol. 25, Issue 11). American Geophysical Union. https://doi.org/10.1029/98GL01228 Thorne, M. S., Garnero, E. J., Jahnke, G., Igel, H., & McNamara, A. K. (2013). Mega ultra low velocity zone and mantle flow. Earth and Planetary Science Letters, 364, 59–67. https://doi.org/10.1016/j.epsl.2012.12.034 83 Thorne, M. S., Leng, K., Pachhai, S., Rost, S., Wicks, J., & Nissen-Meyer, T. (2021). The Most Parsimonious Ultralow-Velocity Zone Distribution From Highly Anomalous SPdKS Waveforms. Geochemistry, Geophysics, Geosystems, 22(1). https://doi.org/10.1029/2020GC009467 To, A., Romanowicz, B., Capdeville, Y., & Takeuchi, N. (2005). 3D effects of sharp boundaries at the borders of the African and Pacific Superplumes: Observation and modeling. Earth and Planetary Science Letters, 233(1–2), 137–153. https://doi.org/10.1016/j.epsl.2005.01.037 Van Keken, P. (1997). EPSL Evolution of starting mantle plumes: a numerical and laboratory. In Earth and Planetary Science Letters (Vol. 148). https://www.sciencedirect.com/science/article/pii/S0012821X97000423 Wicks, J. K., Jackson, J. M., & Sturhahn, W. (2010). Very low sound velocities in iron-rich (Mg,Fe)O: Implications for the core-mantle boundary region. Geophysical Research Letters, 37(15). https://doi.org/10.1029/2010GL043689 Wicks, J. K., Jackson, J. M., Sturhahn, W., & Zhang, D. (2017). Geophysical Research Letters Sound velocity and density of magnesiowüstites: Implications for ultralow-velocity zone topography. Wiley Online Library, 44(5), 2148–2158. https://doi.org/10.1002/2016GL071225 Yu, S., & Garnero, E. J. (2018). Ultralow Velocity Zone Locations: A Global Assessment. Wiley Online Library, 19(2), 396–414. https://doi.org/10.1002/2017GC007281 Yuan, K., & Romanowicz, B. (2017). Seismic evidence for partial melting at the root of major hot spot plumes. Science, 357(6349), 393–397. https://doi.org/10.1126/science.aan0760 Zhong, S. (2006). Constraints on thermochemical convection of the mantle from plume heat flux, plume excess temperature, and upper mantle temperature. Journal of Geophysical Research: Solid Earth, 111(4), 4409. https://doi.org/10.1029/2005JB003972 84 APPENDIX Table 4.1 Effects of grain size on diffusion creep viscosity* Non-dimensional Diffusion creep viscosity of dense material1 Non-dimensional grain size of dense material2 0.1 1 10 100 1,000 0.40 1.00 2.51 6.31 15.85 𝑚 * The table shows how the intrinsic diffusion creep viscosity of dense material (thermochemical pile and ULVZ material) is related to the grain size of it. The calculated values are based on the power-law relationship between grain-size and viscosity in diffusion creep as stated in McNamara et al., 2003: 𝜂𝑑𝑖𝑓𝑓 = 𝜇 𝐴′ 𝑑𝑖𝑓𝑓 ( 𝑑 𝑏 ) exp( 𝑔𝑑𝑖𝑓𝑓𝑇𝑚 𝑇𝑑𝑖𝑚 ), where 𝑇𝑚is the dimensional melting temperature, = ( 𝑑𝑑 𝑑𝑚 𝑇𝑑𝑖𝑚 is the dimensional temperature, d is the grain size, 𝑔𝑑𝑖𝑓𝑓 is an empirical parameter related to activation free energy, 𝜇 is the reference rigidity, b is the reference Burgers vector, and 𝐴′ 𝑑𝑖𝑓𝑓 is a prefactor. m is a constant and is taken as 2.5 according to Karato and Wu, 1993. Note here that 𝜂𝑑 )2.5 , where 𝜂𝑑 and 𝜂𝑚 are the intrinsic viscosity of the dense material and the 𝜂𝑚 background mantle, and 𝑑𝑑 and 𝑑𝑚 are the grain size of the dense material and the background mantle. 1. Non-dimensional diffusion creep viscosity of dense material, taking the diffusion creep viscosity of the background mantle as 1. 2. Non-dimensional grain size, taking the grain size of the background mantle as 1.00. Note: Although grain size affects only the diffusion creep viscosity, because thermochemical piles are low stress regions and are dominated by diffusion creep14, it is a reasonable approximation to multiply the viscosity formula by a composition-dependent prefactor 𝜂𝑐(𝐶) in the model. 85 Table 4.2 Parameters used in this study. Parameters Mantle density ρ0 Reference value 3,500 kg/m3 Gravitational acceleration g 9.8 m/s2 Thermal expansitivity α0 1 × 10-5 K-1 Temperature difference between the core and the surface ∆T 3000 K Mantle thickness h 2891 km Reference viscosity η0 2.5 × 1019 Pa s Thermal diffusivity κ0 1 × 10-6 m2/s 86 Table 4.3 Model parameters for all cases used in this study. Case Ra 1 2 3 4 5 6 7 8 9 1 × 107 1 × 107 1 × 107 1 × 107 1 × 107 5 × 107 1 × 107 5 × 107 1 × 107 10 1 × 107 11 1 × 107 12 1 × 107 A 9.21 9.21 9.21 9.21 9.21 9.21 9.21 9.21 9.21 9.21 9.21 4.61 Bu 2.0 2.0 2.0 2.0 2.0 2.0 2.0 4.0 2.0 2.0 2.0 2.0 Bp 0.8 0.8 0.8 0.8 0.8 0.8 1.2 0.8 0.8 0.8 0.8 0.8 𝜂𝑐(𝐶)u 1.0 𝜂𝑐(𝐶)p 1.0 Figure 4.1 100.0 100.0 4.2, 4.10 10.0 0.1 10.0 0.1 4.6 4.7 1,000.0 1,000.0 4.8, 4.9 1.0 1.0 100.0 100.0 1.0 0.1 0.01 10.0 1.0 1.0 1.0 1.0 1.0 1.0 4.11 4.12 4.13 4.14 4.14 4.14 4.15 Ra: Rayleigh number; A: Activation coefficient; Bu: Buoyancy number of ultra-dense material; Bp: Buoyancy number of compositional reservoirs; 𝜂𝑐(𝐶)u and 𝜂𝑐(𝐶)p: intrinsic compositional viscosity coefficient of ultra-dense material and thermochemical piles. 87 Figure 4.5 Snapshots of cross-sectional shapes of the UDMAs at different time and location than Figure 4.1 in the reference case (Case 1). The snapshot is at 2,305 million years after the start of the calculation, which is 1,170 million years after the snapshot in Figure 1. (a) Superposition of non-dimensional temperature field and compositional field, with light blue means cool and light red means warm. Dark red represents UDMAs (effective buoyancy number B greater than 1.0). Panel (i) – (vii) zoom-in of the regions with UDMAs, as marked in (a). Color in panel (i) –(vii) represents logarithm of non-dimensional viscosity. Red lines outline UDMAs. White lines delineate pile material. Black arrows are unit vectors of mantle flow velocity, indicating flow directions. Grey lines are temperature contours at a step of T = 0.1. 88 Figure 4.6 UDMA cross-sectional shapes in the intrinsically 10x more viscous pile (Case 3), shown at 1482 million years after the start of the calculation. b is zoom-in view of the black box in a. c and d are zoom-ins of the regions outlined with black lines in b. The panels represent the same fields and colors as Figure 4.5. The two UDMAs at pile edges are in quasi-equilibrium. 89 Figure 4.7 UDMA cross-sectional shapes in the intrinsically 10x less viscous pile (Case 4), shown at 515 million years after the start of the calculation. The panels represent the same fields and colors as Figure 4.6 except that the length of the zoom-in panel c-e is ~1.35x longer than that in Figure 1. UDMAs are clearly more flattened compared to other cases. UDMAs at pile edges are highly asymmetrical (c and e), while symmetrical within the pile interior (d). There are 5 discontinuous UDMAs found in the zoomed pile, which is more than 2-3 patches that are found in other cases investigated. Only 1 UDMA within the pile interior is zoomed in, because the other two exhibit similar symmetry to (d). 90 Figure 4.8 Plumes continuously forming (pointed by the yellow arrow) away from thermochemical piles after piles have developed, where piles are intrinsically 1,000 x more viscous than the background mantle (Case5). The left column is the non-dimensional temperature field and the right column is showing the compositional field. (a) At 72 million years, the plume away from piles on its both sides just start to form. On its right side, UDMAs have accumulated at the pile edge adjacent to the plume. Yet neither UDMAs nor pile material are observed beneath the plume. (b) At 877 million years, two plumes in the highlighted region is generating, yet both the position and morphology of the piles on both sides of the region changed very little during the past ~800 million years. UDMA has just been swept to the plume root on the right by the mantle flow. Yet neither UDMA nor pile material are observed beneath the plume on the left. Note that in Case 2, at the similar time (Figure 4.10), more primordial materials have been swept to this region, filling up the root of the plumes so that no plumes are generating directly from the lowermost mantle. 91 Figure 4.9 UDMA could be flattened out at the edge of the intrinsically 1,000x more viscous pile (Case 5), shown at 877 million years after the start of the calculation. The transparent background is temperature field, with blue meaning cool and red meaning warm. The yellow region is the UDMAs. The red contour outlines a primordial pile. Arrows show mantle flow direction, with the length of them scaled to the magnitude of flow velocity. The UDMA is dragged towards the root of the upwelling plumes, but the movement of the majority of the highly viscous pile is slower and shows a delay response to the background mantle flow acting on it, leading to the UDMA being flattened out. 92 Figure 4.10 High viscosity rind formed along the interface between the background mantle and the intrinsically 100x more viscous pile (Case 2), shown at 875 million years after the start of the calculation. Left is the logarithm of non-dimensional viscosity field and right is the compositional field. The yellow arrow points to a region where an upwelling plume is forming on top of the compositional reservoir. Note here that the plume is not rooted directly from the core-mantle boundary. Furthermore, more primordial material has been swept beneath the highlighted upwelling region, as part of the thermochemical pile on its right side. 93 Figure 4.11 Effects of Rayleigh number on the cross-sectional morphology of UDMAs (Case 6), shown at 173 million years after the start of the calculation. The figure is a snapshot of the compositional field of the semi-sphere of the model. The case is identical to Case 1 except that the Rayleigh number of the case is 5.0 x 107. Red represents UDMAs. Pink represents pile material. Black represents the background mantle. The 7 UDMAs are zoomed-in in the white boxes. Regions 4 and 5 are in an area where piles are approaching each other. Region 7 is in a pile merging area. UDMAs in other regions are in a quasi-stable state. 94 Figure 4.12 Effects of intrinsic density for thermochemical piles on the cross-sectional morphology of UDMAs (Case 7), shown at 432 million years after the start of the calculation. The figure is a snapshot of the compositional field of the semi-sphere of the model. The case is identical to Case 2 except that the effective buoyancy number for the thermochemical pile is Bp = 1.2. Red represents UDMAs. Pink represents pile material. Black represents the background mantle. The 4 UDMAs are zoomed-in in the white boxes. UDMAs in region 1 is in quasi-stable state. Region 2 is where two piles just merged and formed into one pile. The UDMAs in regions 3 and 4 are in transition mode to pile edge. 95 Figure 4.13 Effects of intrinsic density of UDMAs on the cross-sectional morphology of UDMAs (Case 8), shown at 337 million years after the start of the calculation. The figure is a snapshot of the compositional field of the semi-sphere of the model. The case is identical to Case 6 except that the effective buoyancy number for the UDMAs is Bu = 4.0. Red represents UDMAs. Pink represents pile material. Black represents the background mantle. The 5 UDMAs are zoomed-in in the white boxes. UDMAs in regions 1 and 5 are in quasi-stable state. Regions 2 and 4 are where two piles just merged and formed into one pile. The UDMA in region 3 is in transition mode to pile edge. 96 Figure 4.14 Effects of intrinsic compositional viscosity for the UDMAs on the cross-sectional morphology of UDMAs (Case 9-11), shown at around the same time, after 873 million years of the start of the calculation. (a), (b), and (c) show the same zoomed-in evolving pile and are snapshots at around the same time after the start of the calculation, so that the comparison here is least affected by the variation of UDMA shapes from different time and location. The three cases are identical to each other except that the intrinsic compositional viscosity of the UDMAs is 10 times less (Case 9), 100 times less (Case 10), and 10 times more (Case 11) than that of the background mantle (𝜂𝑐(𝐶)u = 1) and thermochemical piles (𝜂𝑐(𝐶)u = 1). The pile in the 3 cases all have 2 UDMAs in quasi-equilibrium state at two sides of the pile edges at the time shown. The UDMAs are zoomed-in in white boxes. Red represents UDMAs (effective buoyancy number greater than 1.0), pink represents thermochemical piles (C is between 0.25 to 1.0), and black represents the background mantle (C is less than 0.25). 97 Figure 4.15 Effects of temperature-dependent rheology on the cross-sectional morphology of UDMAs (Case 12), shown at 463 million years after the start of the calculation. The figure is a snapshot of the compositional field of the semi-sphere of the model. The case is identical to Case 1 except it has a different temperature-dependent rheology. In this case, the activation parameter for the temperature dependent viscosity is A = 4.61, leading to 10x viscosity contrast between the coolest and hottest area. The background is compositional field. UDMAs are shown by red (effective buoyancy number greater than 1.5) and are zoomed-in in the white boxes. At the time, the two piles on the sides are moving away from the pile in the center, driven by dowelling flow. UDMA 1 is located where two piles are about to merge and will merge with another UDMA soon. UDMA 5 and 6 were just split up from one UDMA and will reach their equilibrium 273 million years later. Other UDMAs were in quasi-equilibrium state at the time. 98 CHAPTER 5: ARE INTRINSICALLY MORE VISCOUS THERMOCHEMICAL PILES STATIONARY AT THE CORE-MANTLE BOUNDARY? 5.1 Abstract Global seismic tomography models reveal two Large Low Shear Velocity Provinces (LLSVPs) that are hypothesized to be caused by large-scale, compositionally distinct reservoirs (e.g., thermochemical piles) in the Earth’s lowermost mantle. Whether the LLSVPs are fixed or mobile is important in paleogeographic reconstructions and our understanding of the dynamic nature of the Earth’s lower mantle. Recently, paleomagnetic studies (e.g., Burke et al., 2008; Torsvik et al., 2010) have implied that the LLSVPs may remain stationary for several hundred million years, while another study (e.g., Flament et al., 2022) showed that the paleomagnetic data could also be explained by mobile LLSVPs. Although geodynamical studies reveal that thermochemical piles are easily swept around by changing subduction patterns, there are other parameters such as compositional viscosity contrast between piles and the ambient mantle that have not been sufficiently explored in these models. Conventionally, thermochemical calculations employ the same rheological formulation for piles and the background mantle; however, piles may have a higher intrinsic viscosity due to a larger grain size or other compositional differences. In this study, we investigate whether and how composition-dependent rheology can explain the possible lateral fixity of thermochemical piles against changing subduction patterns. We performed 2D annulus geodynamical calculations to explore how the change of intrinsic viscosity of thermochemical piles affects the morphology and lateral fixity of LLSVPs in response to slabs and plumes. We introduce isolated thermal anomalies in the Earth’s mantle to generate changes in convective flow that interact with the thermochemical piles. Cold and hot anomalies are introduced in the upper mantle and lowermost mantle to represent changing subduction and upwelling flow, 99 respectively. As pile viscosity increases, we find that the pile morphology is more stable against changing upwelling flow. However, piles remain laterally mobile in all cases, even if they are up to 5,000x more intrinsically viscous than the ambient mantle. In summary, we find that increased thermochemical pile viscosity stabilizes pile morphology yet does not make them less laterally mobile. If LLSVPs are laterally fixed, another mechanism must be found. 5.2 Introduction Basalts formed at oceanic islands (OIBs) are derived from rocks that are less depleted in incompatible elements than those formed mid-ocean ridge basalts (MORBs) (e.g., Hofmann, 1997; Y. Niu & O’Hara, 2003). display higher and more variable 3He/4He ratios than depleted mid-ocean ridge basalts (MORBs) (e.g., Hofmann, 1997; Houser et al., 2008; Porcelli & Elliott, 2008; White, 2015; Williams et al., 2015). Assuming that OIBs are related to plumes rooted in a deep thermal boundary layer, this geochemical observation suggests the existence of a more primordial chemical reservoir in the lower mantle (e.g., Labrosse et al., 2007; Richard et al., 1976; Tackley, 2000; Tucker & Mukhopadhyay, 2014; Wasserburg & DePaolo, 1979). The structure and dynamics of the reservoir control the large-scale mantle convective flow and core-mantle boundary (CMB) heat flux (e.g., M. Li et al., 2018; Nakagawa et al., 2010), which bring insight into the chemical and thermal evolution of the Earth (e.g., Garnero & McNamara, 2008a; Tackley, 2012). Understanding the dynamic nature of the compositional heterogeneity in the Earth’s lower mantle remains one of the challenging problems in Earth Science. Global seismic tomography models reveal two Large Low Shear Velocity Provinces (LLSVPs) located beneath Africa and the Pacific at the base of the mantle (e.g., Dziewonski et al., 2010; French & Romanowicz, 2014; Grand, 2002; Hosseini et al., 2020; Houser et al., 2008; Lekic et al., 2012; Megnin & Romanowicz, 2000; Ritsema et al., 2011; Schuberth et al., 2009; Simmons 100 et al., 2010). Although the hypotheses that LLSVPs are caused by thermal structures cannot be entirely excluded (e.g., D. R. Davies et al., 2012, 2015), there is seismic evidence suggesting that LLSVPs may be compositionally distinct (e.g., Deschamps & Tackley, 2008; Deschamps & Trampert, 2003; French & Romanowicz, 2014; Ishii & Tromp, 2004; Lekic et al., 2012; Masters et al., 2000; Su & Dziewonski, 1997; Trampert et al., 2004). For instance, waveform studies infer sharp seismic velocity gradients along LLSVP margins (e.g., Ford et al., 2006; Y. He & Wen, 2009, 2012; Ni & Helmberger, 2003a, 2003b; D. Sun et al., 2007; X. Sun et al., 2007; To et al., 2005; Wang & Wen, 2007; Wen et al., 2001). LLSVP regions are also characterized by an anticorrelation of bulk sound velocity to shear-wave velocity (e.g., Antolik et al., 2003; Koelemeijer et al., 2016; Ritsema & van Heijst, 2002; Su & Dziewonski, 1997). Additionally, normal mode studies (e.g., Ishii & Tromp, 1999, 2001, 2004; Koelemeijer et al., 2017) and body tide modeling (e.g., Lau et al., 2017) infer that provinces anomalous in density exist in the lowermost few hundreds of kilometers of the Earth’s mantle. Combining the seismic evidence with the observation that many, although not all, surface hot spot locations overlie the LLSVPs (e.g., Thorne et al., 2004; Torsvik et al., 2006; C. Zhao et al., 2015), the LLSVPs are hypothesized to be caused by thermochemical reservoirs at the base of the Earth’s mantle. Furthermore, because the OIB geochemical signatures imply the presence of a more primordial chemical reservoir at the depth (e.g., Boyet & Carlson, 2005; Deschamps et al., 2011; Lee et al., 2010; Samuel & Farnetani, 2003; Tackley, 2000), one of the hypotheses regarding the origin of LLSVPs proposes that the LLSVPs originated from accumulation of more primordial materials (thermochemical piles) formed early in Earth’s history (e.g., McNamara & Zhong, 2005; Tackley, 1998, 2002, 2012). Thermochemical geodynamic models show that the primordial materials have to be intrinsically 101 slightly denser than the background mantle (e.g., McNamara, 2019; Tackley, 2012), to accumulate at such scale along the CMB. There has been a debate about whether the LLSVPs are stable and laterally fixed at their current locations for up to several hundred million years (e.g., Dziewonski et al., 2010; McNamara, 2019; Torsvik, 2019). In the past decade, some studies used plate reconstruction to date back the past eruption locations of kimberlites and paleo-locations of Large Igneous Provinces (LIPs) up to a few hundred million years ago (e.g., Burke et al., 2008; Burke & Torsvik, 2004; Torsvik et al., 2008, 2010). They found that the reconstructed LIP locations overlay the edges of present-day LLSVPs (Figure 5.1). Assuming that deep plumes form on the peripheries of the LLSVPs and generate surface hotspots and LIPs (e.g., Burke et al., 2008; Heyn et al., 2020; Torsvik et al., 2006), they suggested that past eruption locations of LIPs could be footprints of ancient LLSVP locations. Based on the assumptions, they further hypothesize that LLSVPs have remained at their current places for at least 200 million years. The inference has the potential to provide a powerful tool to greatly refine paleogeographic reconstructions and paleolongitude beyond the past few hundred million years (e.g., Burke et al., 2008; Burke & Torsvik, 2004; Torsvik et al., 2008, 2010, 2014). On the contrary, a recent study by (Flament et al., 2022) indicated that mobile LLSVPs could non- uniquely explain the paleomagnetic data. The research brings the debate to the table again. In contrast to the previous paleomagnetic studies implying fixed LLSVPs, several geodynamic studies indicate that lateral locations and morphologies of LLSVPs are transient configurations. Geodynamic studies that assume LLSVPs are caused by thermochemical piles and employ tectonic plate motion history as surface velocity boundary conditions (e.g., Bower et al., 2013; Flament et al., 2022; Langemeyer et al., 2020, 2022; McNamara & Zhong, 2005; Tan et al., 2011; Zhang et al., 2010) revealed that 102 Figure 5.1 Hotspot locations, past eruption locations of kimberlites, and paleo-locations of large igneous provinces (LIPs) superimposed on a map of SMEAN seismic shear wave tomography model at 2800 km depth. Margins of the two Large Low Shear Velocity Provinces (LLSVPs) are shown by -1% shear wave velocity perturbation contour, outlined with red lines. This figure shows that present-day hotspot locations, paleo-locations of LIPs, and most original locations of kimberlites generally overlap the LLSVP margins. thermochemical piles are passive structures easily swept around by changing subduction patterns. In addition, it is found that as mantle convection converges at upwelling centers, pile materials accommodate their shapes to upwelling currents, causing plumes to form on top of piles (e.g., Heyn et al., 2018; M. Li & Zhong, 2017; Y. Li et al., 2014, 2018; McNamara & Zhong, 2004; Q. Yuan & Li, 2022). At first sight, the lateral fixity of LLSVPs is contrary to geodynamic calculations. However, other parameters, such as compositional viscosity contrast between piles and the background mantle, have not been sufficiently explored in these models. Therefore, it is necessary and significant to study the thermochemical parameters underexplored in traditional 103 geodynamic models and examine whether a feasible physical mechanism exists for the possible lateral fixity of thermochemical piles against changing subduction patterns. Conventionally, thermochemical geodynamic calculations employ the same temperature- dependent rheological formulation for thermochemical piles and the background mantle (e.g., McNamara & Zhong, 2005; Tan et al., 2011; Zhang et al., 2010). However, if LLSVPs are hypothesized to be compositionally distinct from the background mantle, they might have different intrinsic viscosity than the background mantle, leading to different rheological formulations for piles and the background mantle. In particular, the LLSVPs may have an intrinsically higher viscosity than the background mantle. From one aspect, LLSVPs may have a higher intrinsic viscosity due to a larger grain size than the background mantle. It is currently impossible to constrain the grain size in the lower mantle from real Earth samples. However, in some grain size evolution computations (e.g., Dannberg et al., 2017; Solomatov & Reese, 2008) that combine mineral physics studies, it is suggested that the grain size in hotter and more primordial thermochemical piles may have the condition to grow larger than that in the background mantle. In addition, a recent attenuation tomography model (e.g., Deuss et al., 2021) suggested a larger grain size of LLSVPs to explain their weak attenuation and higher density features. A larger grain size will lead to a higher diffusion creep viscosity because it’s sensitive to grain size (e.g., Karato & Wu, 1993; Mcnamara et al., 2003) (Table 5.3). If the LLSVPs have a larger grain size and are assumed to be caused by thermochemical piles, their viscosity reduction due to temperature dependence rheology is reduced. Thus, piles have a higher effective viscosity than the background mantle in some cases. From another aspect, LLSVPs may be enriched in bridgmanite and depleted in ferropericlase (e.g., Trampert et al., 2004), which causes them to be stronger and more viscous than the surrounding mantle (e.g., Ballmer et al., 2017). 104 The importance of composition-dependent rheology in lower mantle dynamics was realized a long time ago (e.g., Karato, 1984; McNamara & Zhong, 2004; Solomatov, 2001; Solomatov & Reese, 2008), but how it affects the stability and shapes of thermochemical piles have not been emphasized until recently (e.g., Langemeyer et al., 2022; Y. Li et al., 2019; Schierjott et al., 2020; Q. Yuan & Li, 2022). With the increased intrinsic pile viscosity, the increased compositional viscosity can counteract or surpass the viscosity reduction due to the temperature dependence of viscosity in the interior of piles. It is found that with raised pile intrinsic viscosity, because the interface between the pile and the background mantle is cooler than the pile interior, a high viscosity rind developed along pile surfaces (e.g., McNamara & Zhong, 2004). This rind may decouple the rheology of piles and the background mantle, leading to different dynamical behavior between them. For instance, it may cause piles to be more resistant to changing convection patterns. Therefore, it is crucial to investigate whether the increased intrinsic viscosity of thermochemical piles could cause lateral fixity and the stability of the cross-sectional morphology of piles against changing convection currents. In this study, we perform 2.5 D annulus thermochemical geodynamic experiments to systematically investigate how dense thermochemical piles with increased pile intrinsic viscosity respond to changing convective flow patterns in terms of lateral mobility and cross-sectional morphology. 5.3 The Experiment Setup The aim of our experiment is to examine whether higher viscosity thermochemical piles can resist changing convection flow patterns. In other words, can higher viscosity piles become more spatially stable? We performed 2 types of experiments to explore how the intrinsic viscosity of thermochemical piles affects the ability for piles to respond to changing mantle flow patterns. 105 In the first type, we examine how piles respond to changing downwelling currents. In the second type, we examine how piles respond to changing upwelling currents. For each of the two types of experiment, we explore 2 different Rayleigh numbers, leading to 4 sets of experiments in total. We have 2 initial conditions (temperature and composition), one for each Rayleigh number (Figure 5.2). These 2 initial conditions were generated by performing 2-component thermochemical convection calculations. The calculations for initial conditions started from a linear temperature profile and a 382 km (~7.2% of the volume of the mantle) thick layer of more-dense material atop the CMB, and reached thermal quasi-equilibrium to form the initial conditions applied later. Details such as Buoyancy and Rayleigh numbers used are provided in Table 5.2. That first calculations (to generate the initial conditions) utilized only temperature- dependent rheology (no composition-dependence). For each set of experiments, we explore five different viscosity pre-factors for the pile material: 1x (no composition-dependence), 100x, 1,000x, 2,000x, and 5,000x. This viscosity pre- factor is applied to rheological formulation of pile material only. When combined with temperature-dependence, the 5 pre-factors effectively lead to the following viscosity contrasts with the background mantle: 0.03x (no composition-dependence), 1.0x, 3.2x, 10 x, and 31.6x (these are approximate values). Figure 5.3 illustrates the effective viscosities of the initial condition when Rayleigh number is 9.34 × 106 , taking into account both temperature- and composition- dependence. Figure 5.10 illustrates the effective viscosities of the initial condition when Rayleigh number is 9.34 × 107, taking into account both temperature- and composition-dependence. For each set of experiments, we perturb the flow by introducing spherical thermal anomalies at the beginning of the calculation and again, at a later, intermediate time during the calculation. Half of the experiments incorporate cold spheres (to induce downwelling flow) and 106 Figure 5.2 Initial temperature and composition conditions for (a) experiment set with Rayleigh number of 9.34 × 106 and (b) 9.34 × 107 . Here, the background of the figures shows non- dimensional temperature fields. Yellow thermochemical pile materials outline compositional field with a buoyancy number no more than 0.4 and are superimposed on the temperature field. half incorporate hot spheres (to induce upwelling). For the cold sphere experiments, we examine how effective the newly induced downwellings are at moving thermochemical piles laterally on the CMB. For the hot sphere experiments, we examine how effective the newly induced upwellings are at changing the cross-sectional morphology of the piles. In summary, our main collection of experiments consists of 20 calculations to investigate both Rayleigh numbers for the 5 different pile viscosity pre-factors. In addition, we explore cases with a different pile buoyancy number, a different temperature-dependent rheology, and different resolutions (Table 5.2). 107 Figure 5.3 The logarithm of non-dimensional initial effective viscosity fields used for the cases in experiment Set A and B in this study. The intrinsic viscosity contrasts between thermochemical piles and the background mantle are none in Case A1 and B1 (a), 100x in Case A2 and B2 (b), 1,000x in case A3 and B3 (c), 2,000x in Case A4 and B4 (d), and 5,000x in Case A5 and B5 (e), respectively. Rayleigh number applied in these sets are 9.34 x 106. 5.4 Methods The geodynamic calculations are performed by solving the following non-dimensional equations derived from the conservation of mass, momentum, and energy under the Boussinesq approximation: ∇ ∙ 𝑢⃑ = 0…………………………………………………………………………….Equation (5.1) −∇𝑃 + ∇ ∙ (𝜂𝜀̇ ̿) = 𝜆𝑅𝑎(𝑇 − 𝐵𝐶)𝑟̂………………………………………………….Equation (5.2) 𝜕𝑇 𝜕𝑡 + (𝑢⃑ ∙ ∇)𝑇 = ∇2𝑇 + 𝑄…………………………………………………………..Equation (5.3) 108 where 𝑢⃑ , P, η, 𝜀̇ ̿, T, C, 𝑟̂, t and Q represent the velocity, dynamic pressure, viscosity, strain rate tensor, temperature, composition, the unit vector in the radial direction, time, and the volumetric internal heating rate, all in non-dimensional units. 𝜆 = (𝑅/ℎ)3 with R as the radius of the Earth and h as the mantle thickness. The Rayleigh number is defined as: 𝑅𝑎 = 𝜌0𝑔𝛼0∆𝑇ℎ3 𝜂0𝜅0 …………………………………………………………….….……Equation (5.4) where 𝜌0 , 𝛼0 , 𝜂0 , 𝜅0 are the dimensional reference values of density, thermal expansivity, viscosity at non-dimensional temperature T = 0.35, and thermal diffusivity. ∆𝑇 is the dimensional temperature contrast between the surface and the CMB. g is the constant of gravitational acceleration. We use values typically employed in mantle convection problems. We employ a Rayleigh number of 𝑅𝑎 = 9.34 × 106 for most cases and 𝑅𝑎 = 9.34 × 107 for a group of complementary cases. The values of the physical parameters in this study are listed in Table 5.1. The density contrast between the thermochemical pile and background mantle materials is controlled by the buoyancy ratio B, which is defined as: 𝐵 = ∆𝜌 𝜌0𝛼0∆𝑇 …………………………………………………………….….….……..Equation (5.5) where ∆𝜌 is the intrinsic density difference between the pile material and the background mantle material. The buoyancy ratio is set at 0.8 (about 2.4% intrinsic density difference) for most cases. A buoyancy number of 0.8 leads to the formation of thermochemical piles. We also explore a case set with higher buoyancy number of 1.2 that includes the reference case (without composition- dependent viscosity) and another case with 2,000x viscosity pre-factor for the pile materials. In this set, we aim to investigate if the increased buoyancy number affects how piles are reshaped by changing upwelling currents. The viscosity is depth-, temperature-, and composition-dependent: 109 Table 5.1 Fixed physical parameters and their values. Parameter Reference temperature Reference thermal expansivity Reference density Acceleration of gravity Radius of the Earth Mantle thickness Reference viscosity Reference thermal diffusivity Viscosity contrast between the upper and lower mantle Non-dimensional internal heating rate Symbol Value ∆𝑇 𝛼0 𝜌0 𝑔 R h 𝜂0 𝜅0 𝜂𝑧 Q 3000 K 1 × 10−5 K-1 3300 kg m-3 9.8 m s-2 6371 km 2891 km 2.34 × 1021 Pa s 1 × 10−6 m2 s-1 50 60 𝜂(𝑇, 𝑧, 𝐶) = 𝜂𝑧(𝑧)𝜂𝑐(𝐶)exp [𝐴(0.35 − 𝑇)]……………………………………….Equation (5.6) where ηz(z) = 1 and 50 for the upper and lower mantle, leading to a 50x viscosity increase from the upper mantle to the lower mantle (e.g., Bertelloni & Gurnis, 1997). 𝜂𝑐(𝐶) is a viscosity pre- factor applied only to the more-dense thermochemical pile material. A is a non-dimensional activation coefficient, and T is the non-dimensional temperature. We use A = 9.21 for most cases so that there is ~104x temperature-dependent viscosity contrast between the hottest and the coldest regions in the model. To perform our numerical experiments, we modified CitcomCU (e.g., Zhong, 2006) to include compositional rheology. More than 9 million tracers are introduced to the models to track the compositional field, using tracer-ratio method (e.g., Tackley & King, 2003). Our models include 2 compositional components: intrinsically more-dense thermochemical pile material and background mantle material. We run calculations in an effectively 2.5-D spherical geometry where the colatitude direction is only composed of a one-element layer. The model parameters are the 110 same along colatitudinal direction and we use this geometry to represent spherical annulus geometry. The side boundary conditions are reflective, meaning that flow and heat do not go through side boundaries. The model domain is divided into 1920 elements in the longitude direction and 192 elements in the radial direction, respectively. The inner and outer radius of the model is 0.55 and 1.00, respectively. The grid provides ~15km resolution near the CMB. All velocity boundary conditions are free-slip. Temperature boundary conditions are isothermal on the top (T = 0) and bottom (T = 1). The models are heated from both below and internally, with a non- dimensional internal heating rate of Q = 60. As described in the previous section (the Experiment Setup), we performed 2 different types of experiments. In the cold sphere experiment set, in order to examine how piles respond to changing dowelling currents, we incorporate cold spheres to induce downwelling flow. The cold spheres introduced into the models have a non-dimensional temperature of T = 0.0, with a radius of ~400 km (Figure 5.4a and 5.4b). We introduce the cold spheres centering at the depth of 445 km near the surface to induce downwelling currents initiated from the upper mantle. In the hot sphere experiment set, in order to examine how piles respond to changing upwelling currents, we incorporate hot spheres to induce upwelling flow. The hot spheres introduced into the models have a non-dimensional temperature of T = 1.0, with a radius of ~573 km (Figure 5.4c and 5.4d). We introduce the hot spheres centering at the depth of 2230 km near the CMB to induce upwelling currents originating from the lowermost mantle. 5.5 Results Table 5.2 lists the parameters used in all the cases studied. Cases starting with the same letter (A-D) are of the same case set. Within each set, the cases have the same parameters except for composition-dependent viscosity pre-factor in pile materials (ηC). Sets A and C are cold sphere 111 experiment sets B and D are hot sphere experiment sets. The time when thermal anomalies are incorporated at the beginning of the calculations is defined as t = 0. Sets A and B have a Rayleigh number of 9.34 × 106, while Sets C and D have the same Rayleigh number of 9.34 × 107. The cases ending with “1” is the reference case for that set, meaning it has no composition-dependent viscosity applied to thermochemical piles. We first show the cold (Set A) and hot (Set B) anomaly experiment sets with Rayleigh number of 9.34 × 106 . Sets A and B utilized the same initial condition and contained five thermochemical piles (Figure 5.2a). Immediately before the introduction of the anomalies, Pile I was in the process of separating into two piles by a downwelling sinking on top of it. Piles II and III each have two relatively stable downwellings sandwiching them from both sides. Pile IV was being pushed clockwise by the downwelling next to it in the counter-clockwise direction. Pile V was being pushed counter-clock wisely by the downwelling in its clockwise direction. Pile IV and Pile V were moving towards each other and would merge into a combined larger pile shortly after. Figure 5.5 shows the results of Set A. The panel on the left shows the initial thermal and compositional condition. The five panels on the right display representative snapshots of Cases A1-A5. On the right panels, each panel box contains snapshots from a case. The first panel shows the reference Case A1, where piles have no composition-dependent viscosity. The other panels show the cases with increased pile intrinsic viscosity. Snapshots on the same column were taken at the same model time for all the cases. Four cold anomalies were introduced into the cases at the start and again at 376 million years. Two anomalies were near Pile I, and the other two were near Pile V. Within the first 376 million years (Figure 5.5), the first set of cold anomalies reached the 112 Table 5.2 Cases studied and varied parameters in the cases. Case Ra 𝜂𝐶 9.34 × 106 1 9.34 × 106 500 A 9.21 9.21 𝐵𝐶 Anomaly typea 0.8 Cold sphere 0.8 Cold sphere Note A1 A2 A3 A4 A5 B1 B2 B3 B4 B5 C1 C2 C3 C4 C5 D1 D2 D3 D4 D5 21 22 23 24 25 9.34 × 106 1,000 9.21 0.8 Cold sphere 9.34 × 106 2,000 9.21 0.8 Cold sphere 9.34 × 106 5,000 9.21 0.8 Cold sphere 9.34 × 106 1 9.34 × 106 500 9.21 9.21 0.8 Hot sphere 0.8 Hot sphere 9.34 × 106 1,000 9.21 0.8 Hot sphere 9.34 × 106 2,000 9.21 0.8 Hot sphere 9.34 × 106 5,000 9.21 0.8 Hot sphere 9.34 × 107 1 9.34 × 107 500 9.21 9.21 0.8 Cold sphere 0.8 Cold sphere 9.34 × 107 1,000 9.21 0.8 Cold sphere 9.34 × 107 2,000 9.21 0.8 Cold sphere 9.34 × 107 5,000 9.21 0.8 Cold sphere 9.34 × 107 1 9.34 × 107 500 9.21 9.21 0.8 Hot sphere 0.8 Hot sphere 9.34 × 107 1,000 9.21 0.8 Hot sphere 9.34 × 107 2,000 9.21 0.8 Hot sphere 9.34 × 107 5,000 9.21 0.8 Hot sphere 9.34 × 107 1 13.82 0.8 Hot sphere 9.34 × 107 2,000 13.82 0.8 Hot sphere 9.34 × 107 1 9.21 1.2 Hot sphere 9.34 × 107 2,000 9.21 1.2 Hot sphere 9.34 × 107 1 9.21 0.8 Hot sphere Spatial resolution is 4x higher than Case D1. 113 Table 5.2 (cont’d) 26 9.34 × 107 1 9.21 0.8 Hot sphere Density of tracers is 4x more than Case D1. a The initial location and temperature of the anomalies introduced are identical in Cases A1-A5, Cases B1-B5, Cases C1-C5, Cases D1-D5, Cases 21 and 22, and Cases 23 and 24, respectively. lowermost mantle, creating a strong downwelling flow that pushed Pile I and Pile V laterally away from each other. By 376 million years, Pile IV and Pile V completely merged into one pile in Case A1, while in Cases A2-A5, they have not merged yet. At 376 million years, Pile I stopped diverging. At the time, we put in the second set of cold anomalies. Afterward, the second set of cold anomalies created a downwelling flow that pushed Pile I and Pile V even further apart. Over the entire calculation, the separation distance between Pile I and Pile V increased by 97, 116, 110, 105, and 66 epicentral degrees in Cases A1-A5, respectively. In all cases (including cases with increased pile viscosity), piles moved laterally away from the induced downwellings substantially. 114 Figure 5.4 Experiment setup for cold anomaly experiment sets (a and b) and hot anomaly experiment sets (c and d). The figure shows snapshots of non-dimensional temperature fields when thermal anomalies were first introduced into the models (t = 0) for experiment Set A (a), B (c), C (b), and D (d), respectively. Compositional fields with buoyancy number BC of no more than 0.4 are colored yellow and imposed on temperature fields. Here, yellow regions represent denser (thermochemical pile) materials. Cold anomalies and hot anomalies implemented are shown by blue circles and red circles outlined by black lines, respectively. a and c share the same initial model, where thermochemical piles are enumerated as I-V in a. b and d share the same initial model, where thermochemical piles are enumerated as I-V in b. 115 Figure 5.5 Snapshots in a time sequence for Cases A1-A5 with 1x, 500x, 1,000x, 2,000x, and 5,000x pile intrinsic viscosity, respectively. Time t is scaled to geologic time after the introduction of cold anomalies and is listed on top of the snapshots. The figure on the left shows the initial model of Set A, which is shared through Case A1-A5. Each panel shows snapshots for the same case, while each column compares snapshots for the 5 cases at the same time. Color scales here have the same meaning as that in Figure 5.4. 116 Results from Set B are shown in Figure 5.6. Set B is identical to Set A, except we introduce hot anomalies instead of cold ones. Three hot anomalies were introduced into the cases in the beginning and again at 62 million years. The anomalies were located near Pile I, between Pile I and Pile II, and near Pile V, respectively. In the reference Case B1, the anomalies induced strong upwelling flows. In response, lateral return flow swept the surrounding lower mantle toward the newly formed upwelling regions. Nearby piles swiftly changed shape to accommodate the new upwelling flow pattern. Piles repositioned themselves to be at the base of the upwelling regions, forming mantle plumes atop them. Unlike Case B1, Cases B2-5 employ an intrinsic viscosity increase within pile material. Two noticeable differences between them and Case B1 occur within the first 60 million years. Firstly, the upwellings induced by the hot anomalies do not rise as quickly as in Case B1, hinting that they may have been slowed down by a more sluggish lower mantle lateral return flow. Furthermore, the piles themselves appear relatively undeformed. Only after introducing the second set of anomalies, reinforcing the upward flows, did the piles in Cases B2 and B3 begin to reshape toward the upwelling centers. The piles in Cases B4 and B5 appear to be only minimally affected (if at all) by the anomaly-induced upwellings. To further examine how pile morphology is affected by anomaly-induced upwelling, zoom-ins are shown in Figure 5.7. We compare reference Case B1 (with no intrinsic viscosity increase) with Case B4 (with 2,000 x intrinsic viscosity increase). In Case B1 (Figure 5.7b-d), pile material is effectively 2-3 orders of magnitude less viscous than background mantle. In Case B4 (Figure 5.7e-g), the piles are effectively ~1.5 more viscous than background mantle. The anomaly induced stronger upwelling flow and, therefore, caused a greater magnitude of lateral return flow in reference Case B1 than in Case B4. As a result, within 42.7 million years, in Case B1, the 117 Figure 5.6 Snapshots in a time sequence for Cases B1-B5 with 1x, 500x, 1,000x, 2,000x, and 5,000x pile intrinsic viscosity, respectively. Time t is scaled to geologic time after incorporating hot anomalies and is listed on top of the snapshots. The figure on the left shows the initial model of set B, which is shared through Case B1-B5. Each panel shows snapshots for the same case, while each column compares snapshots for the 5 cases at the same time. Color scales here have the same meaning as that in Figure 5.4. 118 Figure 5.7 Comparison between the zoomed-in views of sequential snapshots around one of the hot anomalies for Case B1 and B4. a shows the snapshot of non-dimensional temperature field superimposed with compositional field at the initial time for both cases. Black lines in a outline the zoomed-in region. The top (b-d) and bottom (e-g) panels display snapshots of the logarithm of non-dimensional effective viscosity field with mantle flow velocity superimposed (green arrows) for Case B1 and B4, respectively. Denser materials (thermochemical piles) are outlined with black lines in b-g. Snapshots at t = 7.7 million years (b and e), t = 15.5 million years, and t = 42.7 million years are compared here between the reference case (case B1) and the case with intrinsically 2,000 x more viscous piles (Case B4). adjacent pile was quickly swept to the base of the upwelling region and accommodated its shape to the newly-formed upwelling flow. However, in Case B4, the pile remains undeformed and barely laterally moved. We repeat these experiments with increased convective rigor (Sets C and D), employing a Rayleigh number of 9.34 × 107 . These sets have the same initial condition and contain five thermochemical piles (Figure 5.2b). Immediately before the anomalies were introduced, Piles I and II were pushed apart by the downwelling flow between them. Pile III was sandwiched between 2 downwellings, with upwellings forming from both sides of the edges. Pile IV was swept around by the two downwellings on both sides. A small portion of Pile IV was edging toward Pile V in such a process. Pile V had a flat top with two upwellings forming atop it. The two upwellings were moving towards each other and would reshape the top of Pile V to be sharper soon after. 119 Figure 5.8 shows the results of Set C. Three cold anomalies were introduced into the cases at the beginning and again at 271 million years. They were incorporated between Piles III and IV. Within the first 271 million years, the cold anomalies induced a strong downwelling flow that pushed Pile III and IV away from each other. In all cases, as Pile IV was being pushed towards Pile V by the downwellings, the two piles merged into one pile (“new Pile IV”) by 271 million years. Then, at 271 million years, we put in the second set of cold anomalies to reinforce the downwelling current. Afterward, Pile III and new Pile IV were pushed further apart by the lateral flow induced by the reinforced downwellings. Throughout the calculation, the lateral separation distance between Pile III and Pile IV increased by 75-85 epicentral degrees in all cases, including cases with significantly increased pile viscosity. Results from Set D are shown in Figure 5.9. Set D is identical to Set C, except we introduce hot anomalies instead of cold ones. Three hot anomalies were introduced at the start and again at t = 73.9 million years. The three anomalies were positioned between Piles II and III, Piles III and IV, and Piles IV and V, respectively. In the reference Case D1, the anomalies induced a strong upwelling flow. Nearby pile materials quickly changed their shapes to accommodate the new upwelling current. By 81.6 million years, two small new piles with cusp-like tops were formed, separating themselves from Pile II and Pile IV. In addition, Pile V changed its top shape of from flattened to higher and fluffier on the side near the anomaly. In the meantime, lateral return flow swept the pile materials toward the base of newly formed upwelling regions. Unlike Case D1, in Cases D2-D5, piles are intrinsically more viscous than the background mantle. Our results clearly show that when introducing the second set of anomalies, piles in all cases moved laterally toward the anomalies by a similar amount of distance. However, similar to the results from Set B, the two 120 Figure 5.8 Snapshots in a time sequence for Cases C1-C5 with 1x, 500x, 1,000x, 2,000x, and 5,000x pile intrinsic viscosity, respectively. Time t is scaled to geologic time after the introduction of cold anomalies and is listed on top of the snapshots. The figure on the left shows the initial condition of Set C, which is shared through Case C1-C5. Each panel shows snapshots for the same case, while each column compares snapshots for the 5 cases at the same time. Color scales here have the same meaning as that in Figure 4. Rayleigh number applied in these sets is 9.34 x 107. 121 Figure 5.9 Snapshots in a time sequence for Cases D1-D5 with 1x, 500x, 1,000x, 2,000x, and 5,000x pile intrinsic viscosity, respectively. Time t is scaled to geologic time after the incorporation of hot anomalies and is listed on top of the snapshots. The figure on the left shows the initial condition of Set B, which is shared through Case D1-D5. Each panel shows snapshots for the same case, while each column compares snapshots for the 5 cases at the same time. Color scales here have the same meaning as that in Figure 5.4. Rayleigh number applied in these sets is 9.34 x 107. 122 significant differences between Case D1 and Cases D2-D5 exist. Firstly, the upwellings induced by the hot anomalies reach the surface as quickly as within 7.7 million years in the reference case, while they only slightly rise in other cases by the time. Secondly, pile cross-sectional morphologies remain relatively unchanged in Cases D2-D5. The piles were slightly reshaped toward the upwelling center in Cases D2 and D3. The pile shapes were only affected neglectably in Cases D4 and D5 by the anomaly-induced upwellings. We explore increasing the temperature dependence of viscosity in Case 21 and Case 22. Figure 5.12 shows the results of the two cases. Case 21 and Case 22 use the same parameters as that in Case D1 and Case D4, respectively, except they have a higher activation coefficient (A = 13.82), which leads to a 10,000x viscosity contrast due to temperature across the models. The initial thermal and compositional condition that Case 21 and Case 22 shared differs from Case D1 and D4, as shown in Figure 5.11a. We employ no composition-dependent viscosity in Case 21 and 2,000x intrinsic viscosity in piles in Case 22. When combined with temperature dependence, the two pile viscosity pre-factors effectively lead to 0.001x and 0.31x viscosity contrast with the background mantle. Figures 5.11b and 5.11c illustrate the effective viscosities of the initial condition for Case 21 and Case 22, respectively. In Case 21, piles are drastically reshaped by the upwellings newly induced by the anomalies, with an evident amount of pile materials entrained into the newly formed upwellings. In Case 22, with higher pile intrinsic viscosity, unlike Case 21, pile materials only reshape mildly to accommodate changing upwelling currents or are moved laterally by the return flow undergoing minimal deformation. As pile intrinsic viscosity increases, piles maintain their shapes more against changing upwelling currents. Increased temperature dependence of viscosity did not affect how effective the newly induced upwellings were at changing the cross-sectional morphology of the piles. 123 We also explore increasing the buoyancy number of pile materials in Case 23 and Case 24. Results from these two cases are shown in Figure 5.13. Case 23 and Case 24 utilized the same initial thermal and compositional conditions as in Case 21 and Case 22 (Figure 5.11a). The cases use the same parameters as that in Case D1 and Case D4, respectively, except they have a higher buoyancy number (B = 1.2), leading to a ~4% density increase in piles. We employ no composition-dependent viscosity in Case 23 and 2,000x intrinsic viscosity in piles in Case 24. When combined with temperature dependence, the two pile viscosity pre-factors effectively lead to 0.03x and 10x viscosity contrast with the background mantle, similar to Case D1 and D4. Figures 5.11b and 5.11c illustrate the effective viscosities of the initial condition for Case 23 and Case 24, respectively. In Case 23, the cross-sectional morphology of piles was reshaped by newly formed upwellings as swiftly as that in Case D1. In Case 24, as pile intrinsic viscosity is higher, similar to that in Case D4, the pile shapes are minimally affected (if any) by the increased pile intrinsic viscosity. Increased buoyancy number did not influence how piles with varied intrinsic viscosities reshaped themselves to accommodate newly formed upwellings. We find that increased intrinsic viscosity of thermochemical piles does not cause them to be more resistant to changing downwelling patterns. In all the cases of Sets A and C, piles are pushed laterally away from the downwelling flow induced by the cold anomalies. In all cases of both sets, the increase of separation distance between piles is as substantial as around 90 degrees of epicentral distance during a few hundred million years, regardless of the intrinsic viscosity of piles and the Rayleigh number applied in the cases. Furthermore, the increase of separation distance only has less than 2x of difference among cases in the same set. Piles with higher intrinsic viscosity are not less mobile along the CMB because they are being pushed along the free-slip boundary condition. The increased pile intrinsic viscosity does not change this boundary condition 124 and, thus, how piles are being moved along the boundary. In addition, piles with higher intrinsic viscosity can better resist changing their cross-sectional morphology in response to changing upwelling flows in the background mantle. In Sets B and D, in the reference cases with no composition-dependent viscosity, piles are swiftly deformed by the upwelling flow newly induced by the hot anomalies. However, in higher pile intrinsic viscosity cases, pile cross-sectional morphologies are only slightly to minimally (if any) affected by the newly formed upwelling currents. Increasing the intrinsic viscosity of piles decreases the effective Rayleigh number within the piles. The decrease of effective Rayleigh number within the pile interior leads to greater contrast between effective Rayleigh numbers of piles and the background mantle, as the effective Rayleigh number of the background mantle remains unchanged. Therefore, as the background mantle flow interacts with pile materials, piles with higher intrinsic viscosity are less vigorous than the background mantle and, thus, do not respond to the induced upwelling flow as promptly as those in the reference cases. Effectively, increased pile intrinsic viscosity causes cross-sectional morphologies of piles to be more resistant to changing upwelling flow. In the hot anomaly experiment Sets B and D, we find that in cases where piles are intrinsically more viscous, newly induced upwellings can generate away from piles and remain outside of the piles for around 90 million years. However, new upwellings swiftly migrate to the top of piles in reference cases. This is because piles with higher intrinsic viscosity do not respond to the induced upwelling flow as promptly as those in the reference cases. In the lowermost mantle, plume initiate location remains largely controlled by thermal boundary layer instability (e.g., Heyn et al., 2020; M. Li & Zhong, 2017; Olson et al., 1987; Tan et al., 2011). It is very likely that some plumes initiated outside of the LLSVPs, just like the upwellings induced by hot anomalies introduced outside of piles in our experiments. Some studies showed that if LLSVPs are caused 125 by thermochemical piles, plumes are typically rooted on top of piles (e.g., Deschamps et al., 2011; M. Li et al., 2014). Our results provide a possible mechanism for those plumes initiated outside of piles to remain outside of the piles. If the surface hotspots outside the LLSVPs are formed from deep mantle plumes, we could still expect the LLSVPs to be caused by thermochemical piles but perhaps intrinsically more viscous piles. Since plumes outside of the LLSVPs do not sample LLSVPs, this could also explain why some OIBs outside of LLSVP margins have low 3He/4He ratios. 5.6 Discussion and Conclusion In this study, we performed 2.5D annulus thermochemical convection calculations to investigate how thermochemical pile stability is affected by their increased intrinsic viscosity, ηC. We perturbed the convective flow patterns by introducing thermal anomalies and examined pile lateral mobility and cross-sectional morphology. Hereby, we focus on the effects of increased intrinsic viscosity of piles because we envision it to be more exciting and underexplored from geodynamics perspectives than that of decreased intrinsic viscosity of piles. Firstly, due to the temperature-dependence rheology of mantle materials and a higher temperature in thermochemical piles, in conventional geodynamic models, piles are effectively less viscous than the background mantle and appear passively mobile (e.g., McNamara et al., 2010). Thus, if we decrease the intrinsic viscosity of piles, the effective viscosity of piles will be much less than that of the background mantle, which might slightly change the viscosity structure of piles but will not bring new insights into how piles can be more resistant to changing mantle flow patterns. However, suppose we increase the intrinsic viscosity of piles, instead, until the effective viscosity of piles exceeds that of the background mantle. In this case, the effective Rayleigh number within pile materials is higher than that of the background mantle, causing piles 126 to be less vigorous and have a different viscosity structure than piles with lower intrinsic viscosity. Secondly, a previous study (e.g., McNamara & Zhong, 2004) found that piles with increased intrinsic viscosity of 100x will develop a highly viscous rind along the interface between pile materials and the background mantle because the interface is cooler than the pile interior. This viscous envelope has the potential to lead to interesting dynamical scenarios such as less entrainment of pile materials into mantle plumes. The entrainment of pile materials in mantle plumes may carry primordial isotopes (i.e., 3He) and travel to the Earth’s surface (e.g., Deschamps et al., 2011; M. Li et al., 2014), which is important to our understanding of the geochemical isotopic signature in hot spots, especially those overlying the LLSVPs. Although it is beyond the scope of our study to quantify and fully resolve the entrainment of pile materials into upwellings, we do observe a significantly decreased amount of entrainment into newly induced upwellings as pile intrinsic viscosity increases in both of our hot anomaly experimental sets. The observation is consistent with the conclusion from an earlier lab experiment (e.g., Davaille et al., 2002, 2003), arguing that thermochemical plumes with higher viscosity will have a lower entrainment rate of denser materials in the lowermost mantle. In the future, it is worth quantifying and investigating how the entrainment of pile materials into mantle plumes is affected by the increased pile intrinsic viscosity. If we assume that thermochemical piles are dominated by diffusion creep (e.g., Karato et al., 1995; Solomatov et al., 2002) and the composition-dependent viscosity contrast is caused by grain size variance, the five viscosity prefactors for the pile material explored in this study are approximately caused by: 1x, 12x, 16x, 21x, and 30x of grain size increase in pile materials (Table 5.3). Such relevance involves the discussion of two questions. Are the assumptions here reliable? Are the prefactor values applied here reasonable for thermochemical pile materials? The honest 127 answer is that the deformation mechanism, rheological nature, and grain size of LLSVP materials are poorly constrained because there is no direct evidence of lower mantle rock samples (e.g., Fei et al., 2021; Ranalli & Fischer, 1984). The deformation mechanism and grain size of LLSVPs are highly uncertain, but there are some hypotheses. The Earth’s lower mantle is suggested to be dominated by diffusion creep because of the lack of strong seismic anisotropy observations in most lower mantle (e.g., Meade et al., 1995; F. Niu & Perez, 2004). It is thus appropriate to perform calculations assuming that the LLSVPs are dominated by diffusion creep. In addition, an increase of grain size in LLSVPs is possible, as the high temperature within LLSVPs potentially gives mineral grains in LLSVPs the condition to grow larger than that of the background mantle (e.g., Dannberg et al., 2017). An earlier simulation (e.g., Solomatov & Reese, 2008) showed that it was possible for grain size contrasts between primordial grains and recrystallized grains in the Earth’s mantle to be 10 – 100 times and even more in extreme cases. Geodynamical calculations also show that hot regions in the lower mantle could be more viscous due to a larger grain size (e.g., Solomatov et al., 2002). A recent seismic attenuation study (e.g., Deuss et al., 2021) suggested that the LLSVPs might have a larger grain size because temperature variation alone cannot explain the observed weak attenuation in LLSVP regions. Therefore, if we hypothesize that more primordial reservoirs cause the LLSVPs, the parameter space of their increased grain size, as investigated in this study, is within the reasonable range. However, we should note that we do not rule out other possibilities that cause compositional viscosity contrast between LLSVPs and the background mantle but rather explore under such appropriate assumptions. It might be possible to perform this experiment with naturally evolved upwellings and downwellings. However, through plenty of trials and errors, we found our method of systematically introducing thermal anomalies to be most viable for the experiment. On real Earth, 128 plate tectonics leads to a subduction system with changing subduction patterns (e.g., Bower et al., 2013; Stadler et al., 2010; Zhang et al., 2010). Therefore, to investigate the mechanism that possibly causes the LLSVPs to be stationary for a few hundred million years, we should systematically develop a configuration to change mantle flow patterns. This configuration aims to perturb the existing large-scale flow patterns in the background mantle instead of generating upwellings and downwellings as naturally as possible. Furthermore, when comparing cases with varied pile intrinsic viscosities, we benefit from our thermal anomaly method because it induces changing flow patterns that are more predictable and consistent than those through “more natural methods”. In our main sets of 20 calculations, we employ an activation coefficient of A = 9.21, which leads to 1,000x of viscosity contrast in the model due to temperature-dependent viscosity. We do not explore a wide range of values for A, as we do not expect the temperature dependence of viscosity to change the conclusions of our study (Y. Li et al., 2014). Temperature-dependent rheology effectively changes the viscosity difference between thermochemical piles and the background mantle. However, due to the higher temperature in piles, the temperature dependence of viscosity itself will not lead to a higher effective viscosity in piles than in the background mantle, which is the main aspect to investigate in this study. Moreover, how the increased pile intrinsic viscosity affects the effective viscosity contrast between piles and the background mantle is only slightly influenced by the temperature dependence of viscosity. To verify this, we test increasing the temperature dependence of viscosity to 10,000x in a reference case (Case 21) and a case with 2,000x pile intrinsic viscosity (Case 22). As shown in Figures 5.11b and 5.11c, we found that as we increase the pile intrinsic viscosity prefactor by 2,000 times, the effective viscosity contrast between piles and the background mantle increases by around 300 times, which is similar to that 129 of a lower temperature dependence of viscosity (also around 300 times from Case D1 to Case D4). Therefore, when comparing Case 21 and Case 22, we came to the same conclusion that the cross- sectional morphology of piles is more resistant to changing upwelling currents as pile intrinsic viscosity increases. We employ reflective side boundary conditions in all of the calculations. It is possible to employ periodic boundary conditions. However, it is more efficient for us to observe the responses of piles to changing flow patterns if we employ reflective boundary conditions. Pile materials near both side boundaries are not mixed when using reflective boundary conditions and prohibiting mantle flow from going across the side boundaries. Therefore, the two piles near the side boundaries will never be merged into one pile and will always be considered two separate piles, which gives us an additional pile sample to observe. Furthermore, estimating and comparing lateral moving distance and cross-sectional morphology changes of piles in cases with different pile intrinsic viscosities is more straightforward with this boundary condition. To ensure that we capture the correct deformation style in the calculations, we ran two resolution tests: one with four times higher tracer density (Case 26) and the other with four times higher grid resolution (Case 25). The two cases use the same parameters as Case D1. Although there are expected second-order differences, Case 25, Case 26, and Case D1 all lead to a swift change of pile cross-sectional shapes in response to changing upwelling flow patterns and entrainment of pile materials into newly-formed upwellings. It is under debate whether the LLSVPs have been fixed at their current positions for up to a few hundred million years (e.g., McNamara, 2019). In studies by (e.g., Torsvik et al., 2010), they found that the original eruption locations of LIPs determined by paleomagnetic plate reconstruction overlie the margins of the present-day LLSVPs and hypothesize that the LLSVPs 130 are stationary for at least 200 and up to 500 million years. Based on this hypothesis, LLSVPs can be used as a reference frame in plate reconstruction, potentially providing a powerful tool for paleolongitude reconstruction, as the inclination recorded in magnetite rocks only gives clues to paleolatitude. The hypothesis itself seems inconsistent with what we conventionally observe in geodynamic models because it is revealed in many geodynamic models (e.g., McNamara et al., 2010; Tackley, 2002) that thermochemical piles appear as a passive structure being easily swept around by changing subduction patterns. (e.g., Conrad et al., 2013) performed geodynamical calculations with the buoyancy field derived from shear-wave tomography to model instantaneous mantle flow patterns and provided some dynamical correlation to this “stationary LLSVPs” hypothesis. They found that present-day plate motions and mantle convection patterns share similar dipole and quadrupole divergence/convergence characteristics, with LLSVP regions generally correlated with divergence poles. They then show that the divergence quadrupoles have remained over the LLSVP regions by examining the divergence poles for plate motions over the past 250 million years. They thus suggested that this may provide evidence for LLSVPs being stationary for the past few hundred million years. Although intriguing, this study did not provide a dynamical mechanism for LLSVPs to be at their current positions for as long as 250 million years. When discussing and utilizing the hypothesis of “stationary LLSVPs”, we should be aware of the multiple steps of assumptions it is based upon. Two of the important assumptions are that plumes mostly form at pile edges and the surface expression of plumes (hotspots) is the vertical projection of the plume generation location in the lowermost mantle (e.g., Burke et al., 2008). These assumptions are controversial from dynamical perspectives (e.g., Austermann et al., 2014; Davies, Goes, & Sambridge, 2015; Torsvik et al., 2010). Just as other dynamic studies (e.g., M. Li & Zhong, 2017; McNamara & Zhong, 2004; Tan, Leng, Zhong, & Gurnis, 2011) have shown, our 131 models show that the transient location of plumes can be outside, on top, and at edges of thermochemical piles. It is dynamically very possible that the locations of LIPs in the past do not always reflect the projection of the past LLSVP margins. Recently, (e.g., Flament et al., 2022) reconstructed mantle flow models from the past with various input parameters and compared the predicted basal lower mantle structure (assumed to be the cause of the LLSVPs) and present-day LLSVPs in multiple seismic tomographic models. They found that statistically, the past locations of hotspots would reconcile with mobile basal lower mantle structures, just like they would with fixed ones. Their results reignite the interest in the debate of whether “the stationary LLSVPs” are necessary to explain reconstructions of the history of volcanism. In conclusion, our study shows that even thermochemical piles with ultra-high intrinsic viscosity (5,000x) are easily moved by changing subduction patterns. The increased intrinsic pile viscosity does not cause piles to be laterally stationary for a few hundred million years. Therefore, if LLSVPs are fixed, another mechanism must be found to explain their lateral spatial stability. In addition, we find that the increased intrinsic pile viscosity stabilizes the cross-sectional morphology of piles and makes them more resistant to getting entrained into upwellings. This effect should be considered in future studies on the morphology evolution of LLSVPs (e.g., the morphologic difference between the African and Pacific LLSVP) and the cycling of more primitive materials. 132 REFERENCES Antolik, M., Gu, Y. J., Ekström, G., & Dziewonski, A. M. (2003). J362D28: A new joint model of compressional and shear velocity in the Earth’s mantle. Geophysical Journal International, 153(2), 443–466. https://doi.org/10.1046/J.1365-246X.2003.01910.X/2/153- 2-443-FIG018.JPEG Austermann, J., Kaye, B. T., Mitrovica, J. X., & Huybers, P. (2014). A statistical analysis of the correlation between large igneous provinces and lower mantle seismic structure. Geophysical Journal International, 197(1), 1–9. https://doi.org/10.1093/GJI/GGT500 Ballmer, M. D., Houser, C., Hernlund, J. W., Wentzcovitch, R. M., & Hirose, K. (2017). Persistence of strong silica-enriched domains in the Earth’s lower mantle. Nature Geoscience, 10(3), 236–240. https://doi.org/10.1038/ngeo2898 Bertelloni, C. L., & Gurnis, M. (1997). Cenozoic subsidence and uplift of continents from time- varying dynamic topography. Geology, 25(8), 735. https://doi.org/10.1130/0091- 7613(1997)025<0735:CSAUOC>2.3.CO;2 Bower, D. J., Gurnis, M., Seton, M., Bower, D. J., Gurnis, M., & Seton, M. (2013). Lower mantle structure from paleogeographically constrained dynamic Earth models. Geochemistry, Geophysics, Geosystems, 14(1), 44–63. https://doi.org/10.1029/2012GC004267 Boyet, M., & Carlson, R. W. (2005). Geochemistry: 142Nd evidence for early (>4.53 Ga) global differentiation of the silicate earth. Science, 309(5734), 576–581. https://doi.org/10.1126/science.1113634 Burke, K., Steinberger, B., Torsvik, T. H., & Smethurst, M. A. (2008). Plume Generation Zones at the margins of Large Low Shear Velocity Provinces on the core–mantle boundary. Earth and Planetary Science Letters, 265(1–2), 49–60. https://doi.org/10.1016/J.EPSL.2007.09.042 Burke, K., & Torsvik, T. H. (2004). Derivation of Large Igneous Provinces of the past 200 million years from long-term heterogeneities in the deep mantle. Earth and Planetary Science Letters, 227(3–4), 531–538. https://doi.org/10.1016/j.epsl.2004.09.015 Conrad, C. P., Steinberger, B., & Torsvik, T. H. (2013). Stability of active mantle upwelling revealed by net characteristics of plate tectonics. Nature, 498(7455), 479–482. https://doi.org/10.1038/NATURE12203 Dannberg, J., Eilon, Z., Faul, U., Gassmöller, R., Moulik, P., & Myhill, R. (2017). The importance of grain size to mantle dynamics and seismological observations. Geochemistry, Geophysics, Geosystems, 18(8), 3034–3061. https://doi.org/10.1002/2017GC006944 Davaille, A., Girard, F., & Le Bars, M. (2002). How to anchor hotspots in a convecting mantle? Earth and Planetary Science Letters, 203(2), 621–634. https://doi.org/10.1016/S0012- 821X(02)00897-X 133 Davaille, A., Le Bars, M., & Carbonne, C. (2003). Thermal convection in a heterogeneous mantle. Comptes Rendus - Geoscience, 335(1), 141–156. https://doi.org/10.1016/S1631- 0713(03)00003-8 Davies, D. R., Goes, S., Davies, J. H., Schuberth, B. S. A., Bunge, H. P., & Ritsema, J. (2012). Reconciling dynamic and seismic models of Earth’s lower mantle: The dominant role of thermal heterogeneity. Earth and Planetary Science Letters, 353–354, 253–269. https://doi.org/10.1016/J.EPSL.2012.08.016 Davies, D. R., Goes, S., & Lau, H. C. P. (2015). Thermally Dominated Deep Mantle LLSVPs: A Review. The Earth’s Heterogeneous Mantle: A Geophysical, Geodynamical, and Geochemical Perspective, 441–477. https://doi.org/10.1007/978-3-319-15627-9_14 Davies, D. R., Goes, S., & Sambridge, M. (2015). On the relationship between volcanic hotspot locations, the reconstructed eruption sites of large igneous provinces and deep mantle seismic structure. Earth and Planetary Science Letters, 411, 121–130. https://doi.org/10.1016/j.epsl.2014.11.052 Deschamps, F., Kaminski, E., & Tackley, P. J. (2011). A deep mantle origin for the primitive signature of ocean island basalt. Nature Geoscience 2011 4:12, 4(12), 879–882. https://doi.org/10.1038/ngeo1295 Deschamps, F., & Tackley, P. J. (2008). Searching for models of thermo-chemical convection that explain probabilistic tomography: I. Principles and influence of rheological parameters. Physics of the Earth and Planetary Interiors, 171(1–4), 357–373. https://doi.org/10.1016/J.PEPI.2008.04.016 Deschamps, F., & Trampert, J. (2003). Mantle tomography and its relation to temperature and composition. Physics of the Earth and Planetary Interiors, 140(4), 277–291. https://doi.org/10.1016/J.PEPI.2003.09.004 Deuss, A., Jagt, L., Schneider, S., van Tent, R., Talavera-Soza, S., Deuss, A., Jagt, L., Schneider, S., van Tent, R., & Talavera-Soza, S. (2021). Constraining the nature of the LLSVP’s using whole Earth oscillations. AGUFM, 2021, DI14A-02. https://ui.adsabs.harvard.edu/abs/2021AGUFMDI14A..02D/abstract Dziewonski, A. M., Lekic, V., & Romanowicz, B. A. (2010). Mantle Anchor Structure: An argument for bottom up tectonics. Earth and Planetary Science Letters, 299(1–2), 69–79. https://doi.org/10.1016/j.epsl.2010.08.013 Fei, H., Faul, U., & Katsura, T. (2021). The grain growth kinetics of bridgmanite at the topmost lower mantle. Earth and Planetary Science Letters, 561, 116820. https://doi.org/10.1016/J.EPSL.2021.116820 Flament, N., Bodur, Ö. F., Williams, S. E., & Merdith, A. S. (2022). Assembly of the basal mantle structure beneath Africa. 603(March). https://doi.org/10.1038/s41586-022-04538-y Ford, S. R., Garnero, E. J., & Mcnamara, A. K. (2006). A strong lateral shear velocity gradient 134 and anisotropy heterogeneity in the lowermost mantle beneath the southern Pacific. J. Geophys. Res, 111(3), 3306. https://doi.org/10.1029/2004JB003574 French, S. W., & Romanowicz, B. A. (2014). Whole-mantle radially anisotropic shear velocity structure from spectral-element waveform tomography. Geophysical Journal International, 199(3), 1303–1327. https://doi.org/10.1093/gji/ggu334 Garnero, E. J., & McNamara, A. K. (2008). Structure and Dynamics of Earth’s Lower Mantle. Science, 320(5876), 626–628. https://doi.org/10.1126/science.1148028 Grand, S. P. (2002). Mantle shear–wave tomography and the fate of subducted slabs. Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 360(1800), 2475–2491. https://doi.org/10.1098/rsta.2002.1077 He, Y., & Wen, L. (2009). Structural features and shear‐velocity structure of the “Pacific Anomaly” J. Geophys. Res, 114(2), 2309. https://doi.org/10.1029/2008JB005814 He, Y., & Wen, L. (2012). Geographic boundary of the “Pacific Anomaly” and its geometry and transitional structure in the north. Journal of Geophysical Research: Solid Earth, 117(9). https://doi.org/10.1029/2012JB009436 Heyn, B. H., Conrad, C. P., & Trønnes, R. G. (2018). Stabilizing Effect of Compositional Viscosity Contrasts on Thermochemical Piles. Geophysical Research Letters, 45(15), 7523– 7532. https://doi.org/10.1029/2018GL078799 Heyn, B. H., Conrad, C. P., & Trønnes, R. G. (2020). How Thermochemical Piles Can (Periodically) Generate Plumes at Their Edges. Journal of Geophysical Research: Solid Earth, 125(6), e2019JB018726. https://doi.org/10.1029/2019JB018726 Hofmann, A. W. (1997). Mantle geochemistry: the message from oceanic volcanism. Nature 1997 385:6613, 385(6613), 219–229. https://doi.org/10.1038/385219a0 Hosseini, K., Sigloch, K., Tsekhmistrenko, M., Zaheri, A., Nissen-Meyer, T., & Igel, H. (2020). Global mantle structure from multifrequency tomography using P, PP and P-diffracted waves. Geophysical Journal International, 220(1), 96–141. https://doi.org/10.1093/gji/ggz394 Houser, C., Masters, G., Shearer, P., & Laske, G. (2008). Shear and compressional velocity models of the mantle from cluster analysis of long-period waveforms. Geophysical Journal International, 174(1), 195–212. https://doi.org/10.1111/J.1365-246X.2008.03763.X/3/174- 1-195-FIG016.JPEG Ishii, M., & Tromp, J. (1999). Normal-mode and free-air gravity constraints on lateral variations in velocity and density of earth’s mantle. Science, 285(5431), 1231–1236. https://doi.org/10.1126/SCIENCE.285.5431.1231/SUPPL_FILE/1040466S9_THUMB.GIF 135 Ishii, M., & Tromp, J. (2001). Even-degree lateral variations in the Earth’s mantle constrained by free oscillations and the free-air gravity anomaly. Geophysical Journal International, 145(1), 77–96. https://doi.org/10.1111/j.1365-246X.2001.00385.x Ishii, M., & Tromp, J. (2004). Constraining large-scale mantle heterogeneity using mantle and inner-core sensitive normal modes. Physics of the Earth and Planetary Interiors, 146(1–2), 113–124. https://doi.org/10.1016/j.pepi.2003.06.012 Karato, S. I. (1984). Grain-size distribution and rheology of the upper mantle. Tectonophysics, 104(1–2), 155–176. https://doi.org/10.1016/0040-1951(84)90108-2 Karato, S. I., & Wu, P. (1993). Rheology of the upper mantle: A synthesis. Science, 260(5109), 771–778. https://doi.org/10.1126/science.260.5109.771 Karato, S. I., Zhang, S., & Wenk, H. R. (1995). Superplasticity in Earth’s Lower Mantle: Evidence from Seismic Anisotropy and Rock Physics. Science, 270(5235), 458–461. https://doi.org/10.1126/SCIENCE.270.5235.458 Koelemeijer, P., Deuss, A., & Ritsema, J. (2017). Density structure of Earth’s lowermost mantle from Stoneley mode splitting observations. Nature Communications, 8. https://doi.org/10.1038/ncomms15241 Koelemeijer, P., Ritsema, J., Deuss, A., & van Heijst, H. J. (2016). SP12RTS: A degree-12 model of shear- and compressional-wave velocity for Earth’s mantle. Geophysical Journal International, 204(2), 1024–1039. https://doi.org/10.1093/GJI/GGV481 Labrosse, S., Hernlund, J. W., & Coltice, N. (2007). A crystallizing dense magma ocean at the base of the Earth’s mantle. Nature, 450(7171), 866–869. https://doi.org/10.1038/NATURE06355 Langemeyer, S. M., Lowman, J. P., & Tackley, P. J. (2020). The dynamics and impact of compositionally originating provinces in a mantle convection model featuring rheologically obtained plates. Geophysical Journal International, 220(3), 1700–1716. https://doi.org/10.1093/GJI/GGZ497 Langemeyer, S. M., Lowman, J. P., & Tackley, P. J. (2022). Contrasts in 2-D and 3-D system behaviour in the modelling of compositionally originating LLSVPs and a mantle featuring dynamically obtained plates. Geophysical Journal International, 230(3), 1751–1774. https://doi.org/10.1093/GJI/GGAC143 Lau, H. C. P., Mitrovica, J. X., Davis, J. L., Tromp, J., Yang, H. Y., & Al-Attar, D. (2017). Tidal tomography constrains Earth’s deep-mantle buoyancy. Nature, 551(7680), 321–326. https://doi.org/10.1038/nature24452 Lee, C. T. A., Luffi, P., Höink, T., Li, J., Dasgupta, R., & Hernlund, J. (2010). Upside-down differentiation and generation of a primordial lower mantle. Nature, 463(7283), 930–933. https://doi.org/10.1038/NATURE08824 136 Lekic, V., Cottaar, S., Dziewonski, A., & Romanowicz, B. (2012). Cluster analysis of global lower mantle tomography: A new class of structure and implications for chemical heterogeneity. Earth and Planetary Science Letters, 357–358, 68–77. https://doi.org/10.1016/J.EPSL.2012.09.014 Li, M., McNamara, A. K., & Garnero, E. J. (2014). Chemical complexity of hotspots caused by cycling oceanic crust through mantle reservoirs. Nature Geoscience, 7(5), 366–370. https://doi.org/10.1038/ngeo2120 Li, M., & Zhong, S. (2017). The source location of mantle plumes from 3D spherical models of mantle convection. Earth and Planetary Science Letters, 478, 47–57. https://doi.org/10.1016/j.epsl.2017.08.033 Li, M., Zhong, S., & Olson, P. (2018). Linking lowermost mantle structure, core-mantle boundary heat flux and mantle plume formation. Physics of the Earth and Planetary Interiors, 277, 10–29. https://doi.org/10.1016/j.pepi.2018.01.010 Li, Y., Deschamps, F., & Tackley, P. J. (2014). The stability and structure of primordial reservoirs in the lower mantle: Insights from models of thermochemical convection in three- dimensional spherical geometry. Geophysical Journal International, 199(2), 914–930. https://doi.org/10.1093/GJI/GGU295 Li, Y., Deschamps, F., Yang, J., Chen, L., Zhao, L., & Tackley, P. J. (2019). Effects of the Compositional Viscosity Ratio on the Long-Term Evolution of Thermochemical Reservoirs in the Deep Mantle. Geophysical Research Letters, 46(16), 9591–9601. https://doi.org/10.1029/2019GL083668 Li, Y., Vilella, K., Deschamps, F., Zhao, L., & J. Tackley, P. (2018). Effects of Iron Spin Transition on the Structure and Stability of Large Primordial Reservoirs in Earth’s Lower Mantle. Geophysical Research Letters, 45(12), 5918–5928. https://doi.org/10.1029/2018GL078125 Masters, G., Laske, G., Bolton, H., & Dziewonski, A. (2000). The relative behavior of shear velocity, bulk sound speed, and compressional velocity in the mantle: Implications for chemical and thermal structure. In Geophysical Monograph Series (Vol. 117, pp. 63–87). Blackwell Publishing Ltd. https://doi.org/10.1029/GM117p0063 McNamara, A. K. (2019). A review of large low shear velocity provinces and ultra low velocity zones. Tectonophysics, 760, 199–220. https://doi.org/10.1016/j.tecto.2018.04.015 McNamara, A. K., Garnero, E. J., & Rost, S. (2010). Tracking deep mantle reservoirs with ultra- low velocity zones. Earth and Planetary Science Letters, 299(1–2), 1–9. https://doi.org/10.1016/j.epsl.2010.07.042 Mcnamara, A. K., Van Keken, P. E., Karato, I., Mcnamara, C. :, Van Keken, P. E., & Karato, S.- I. (2003). Development of finite strain in the convecting lower mantle and its implications for seismic anisotropy. J. Geophys. Res, 108(B5), 2230. https://doi.org/10.1029/2002JB001970 137 McNamara, A. K., & Zhong, S. (2004). Thermochemical structures within a spherical mantle: Superplumes or piles? Journal of Geophysical Research: Solid Earth, 109(B7). McNamara, A. K., & Zhong, S. (2005). Thermochemical structures beneath Africa and the Pacific Ocean. Nature, 437(7062), 1136–1139. https://doi.org/10.1038/NATURE04066 Meade, C., Silver, P. G., & Kaneshima, S. (1995). Laboratory and seismological observations of lower mantle isotropy. Geophysical Research Letters, 22(10), 1293–1296. https://doi.org/10.1029/95GL01091 Megnin, C., & Romanowicz, B. (2000). The three-dimensional shear velocity structure of the mantle from the inversion of body, surface and higher-mode waveforms. Geophysical Journal International, 143(3), 709–728. https://doi.org/10.1046/J.1365- 246X.2000.00298.X/2/M_143-3-709-IEQ052.JPEG Nakagawa, T., Tackley, P. J., Deschamps, F., & Connolly, J. A. D. (2010). The influence of MORB and harzburgite composition on thermo-chemical mantle convection in a 3-D spherical shell with self-consistently calculated mineral physics. Earth and Planetary Science Letters, 296(3–4), 403–412. https://doi.org/10.1016/j.epsl.2010.05.026 Ni, S., & Helmberger, D. V. (2003a). Further constraints on the African superplume structure. Physics of the Earth and Planetary Interiors, 140(1–3), 243–251. https://doi.org/10.1016/j.pepi.2003.07.011 Ni, S., & Helmberger, D. V. (2003b). Ridge-like lower mantle structure beneath South Africa. Journal of Geophysical Research: Solid Earth, 108(B2). https://doi.org/10.1029/2001jb001545 Niu, F., & Perez, A. M. (2004). Seismic anisotropy in the lower mantle: A comparison of waveform splitting of SKS and SKKS. Geophysical Research Letters, 31(24), 1–4. https://doi.org/10.1029/2004GL021196 Niu, Y., & O’Hara, M. J. (2003). Origin of ocean island basalts: A new perspective from petrology, geochemistry, and mineral physics considerations. Journal of Geophysical Research: Solid Earth, 108(B4), 2209. https://doi.org/10.1029/2002JB002048 Olson, P., Schubert, G., & Anderson, C. (1987). Plume formation in the D″-layer and the roughness of the core–mantle boundary. Nature 1987 327:6121, 327(6121), 409–413. https://doi.org/10.1038/327409a0 Porcelli, D., & Elliott, T. (2008). The evolution of He Isotopes in the convecting mantle and the preservation of high 3He/4He ratios. Earth and Planetary Science Letters, 269(1–2), 175– 185. https://doi.org/10.1016/j.epsl.2008.02.002 Ranalli, G., & Fischer, B. (1984). Diffusion creep, dislocation creep, and mantle rheology. Physics of the Earth and Planetary Interiors, 34(1–2), 77–84. https://doi.org/10.1016/0031- 9201(84)90086-4 138 Richard, P., Shimizu, N., & Allègre, C. J. (1976). 143Nd/146Nd, a natural tracer: an application to oceanic basalts. Earth and Planetary Science Letters, 31(2), 269–278. https://doi.org/10.1016/0012-821X(76)90219-3 Ritsema, J., Deuss, A., Van Heijst, H. J., & Woodhouse, J. H. (2011). S40RTS: A degree-40 shear-velocity model for the mantle from new Rayleigh wave dispersion, teleseismic traveltime and normal-mode splitting function measurements. Geophysical Journal International, 184(3), 1223–1236. https://doi.org/10.1111/j.1365-246X.2010.04884.x Ritsema, J., & van Heijst, H. J. (2002). Constraints on the correlation of P- and S-wave velocity heterogeneity in the mantle from P, PP, PPP and PKPab traveltimes. Geophysical Journal International, 149(2), 482–489. https://doi.org/10.1046/J.1365- 246X.2002.01631.X/2/M_149-2-482-EQ057.JPEG Samuel, H., & Farnetani, C. G. (2003). Thermochemical convection and helium concentrations in mantle plumes. Earth and Planetary Science Letters, 207(1–4), 39–56. https://doi.org/10.1016/S0012-821X(02)01125-1 Schierjott, J., Rozel, A., & Tackley, P. (2020). On the self-regulating effect of grain size evolution in mantle convection models: Application to thermochemical piles. Solid Earth, 11(3), 959–982. https://doi.org/10.5194/SE-11-959-2020 Schuberth, B. S. A., Bunge, H.-P., & Ritsema, J. (2009). Tomographic filtering of high- resolution mantle circulation models: Can seismic heterogeneity be explained by temperature alone? Geochemistry, Geophysics, Geosystems, 10(5), n/a-n/a. https://doi.org/10.1029/2009GC002401 Simmons, N. A., Forte, A. M., Boschi, L., & Grand, S. P. (2010). GyPSuM: A joint tomographic model of mantle density and seismic wave speeds. Journal of Geophysical Research: Solid Earth, 115(B12), 12310. https://doi.org/10.1029/2010JB007631 Solomatov, V. S. (2001). Grain size-dependent viscosity convection and the thermal evolution of the Earth. Earth and Planetary Science Letters, 191(3–4), 203–212. https://doi.org/10.1016/S0012-821X(01)00426-5 Solomatov, V. S., El-Khozondar, R., & Tikare, V. (2002). Grain size in the lower mantle: constraints from numerical modeling of grain growth in two-phase systems. Physics of the Earth and Planetary Interiors, 129(3–4), 265–282. https://doi.org/10.1016/S0031- 9201(01)00295-3 Solomatov, V. S., & Reese, C. C. (2008). Grain size variations in the Earth’s mantle and the evolution of primordial chemical heterogeneities. Journal of Geophysical Research: Solid Earth, 113(B7), 7408. https://doi.org/10.1029/2007JB005319 Stadler, G., Gurnis, M., Burstedde, C., Wilcox, L. C., Alisic, L., & Ghattas, O. (2010). The dynamics of plate tectonics and mantle flow: From local to global scales. Science, 329(5995), 1033–1038. https://doi.org/10.1126/SCIENCE.1191223/SUPPL_FILE/STADLER-SOM.PDF 139 Su, W. J., & Dziewonski, A. M. (1997). Simultaneous inversion for 3-D variations in shear and bulk velocity in the mantle. Physics of the Earth and Planetary Interiors, 100(1–4), 135– 156. https://doi.org/10.1016/S0031-9201(96)03236-0 Sun, D., Tan, E., Helmberger, D., & Gurnis, M. (2007). Seismological support for the metastable superplume model, sharp features, and phase changes within the lower mantle. Proceedings of the National Academy of Sciences of the United States of America, 104(22), 9151–9155. https://doi.org/10.1073/PNAS.0608160104 Sun, X., Song, X., Zheng, S., & Helmberger, D. V. (2007). Evidence for a chemical-thermal structure at base of mantle from sharp lateral P-wave variations beneath Central America. Proceedings of the National Academy of Sciences of the United States of America, 104(1), 26–30. https://doi.org/10.1073/PNAS.0609143103 Tackley, P. J. (1998). Self-consistent generation of tectonic plates in three-dimensional mantle convection. Earth and Planetary Science Letters, 157(1–2), 9–22. https://doi.org/10.1016/S0012-821X(98)00029-6 Tackley, P. J. (2000). Mantle convection and plate tectonics: Toward an integrated physical and chemical theory. In Science (Vol. 288, Issue 5473, pp. 2002–2007). American Association for the Advancement of Science. https://doi.org/10.1126/science.288.5473.2002 Tackley, P. J. (2002). Strong heterogeneity caused by deep mantle layering. Geochemistry, Geophysics, Geosystems, 3(4), 1–22. https://doi.org/10.1029/2001GC000167 Tackley, P. J. (2012). Dynamics and evolution of the deep mantle resulting from thermal, chemical, phase and melting effects. Earth-Science Reviews, 110(1–4), 1–25. https://doi.org/10.1016/J.EARSCIREV.2011.10.001 Tackley, P. J., & King, S. D. (2003). Testing the tracer ratio method for modeling active compositional fields in mantle convection simulations. Geochemistry, Geophysics, Geosystems, 4(4). https://doi.org/10.1029/2001GC000214 Tan, E., Leng, W., Zhong, S., Gurnis, M., Tan, E., Leng, W., Zhong, S., & Gurnis, M. (2011). On the location of plumes and lateral movement of thermochemical structures with high bulk modulus in the 3-D compressible mantle. Geochemistry, Geophysics, Geosystems, 12(7). https://doi.org/10.1029/2011GC003665 Thorne, M. S., Garnero, E. J., & Grand, S. P. (2004). Geographic correlation between hot spots and deep mantle lateral shear-wave velocity gradients. Physics of the Earth and Planetary Interiors, 146(1–2), 47–63. https://doi.org/10.1016/J.PEPI.2003.09.026 To, A., Romanowicz, B., Capdeville, Y., & Takeuchi, N. (2005). 3D effects of sharp boundaries at the borders of the African and Pacific Superplumes: Observation and modeling. Earth and Planetary Science Letters, 233(1–2), 137–153. https://doi.org/10.1016/j.epsl.2005.01.037 Torsvik, T. H. (2019). Earth history: A journey in time and space from base to top. 140 Tectonophysics, 760, 297–313. https://doi.org/10.1016/J.TECTO.2018.09.009 Torsvik, T. H., Burke, K., Steinberger, B., Webb, S. J., & Ashwal, L. D. (2010). Diamonds sampled by plumes from the core–mantle boundary. Nature 2010 466:7304, 466(7304), 352–355. https://doi.org/10.1038/nature09216 Torsvik, T. H., Smethurst, M. A., Burke, K., & Steinberger, B. (2006). Large igneous provinces generated from the margins of the large low-velocity provinces in the deep mantle. Geophysical Journal International, 167(3), 1447–1460. https://doi.org/10.1111/j.1365- 246X.2006.03158.x Torsvik, T. H., Smethurst, M. A., Burke, K., & Steinberger, B. (2008). Long term stability in deep mantle structure: Evidence from the ~ 300 Ma Skagerrak-Centered Large Igneous Province (the SCLIP). Earth and Planetary Science Letters, 267(3–4), 444–452. https://doi.org/10.1016/J.EPSL.2007.12.004 Torsvik, T. H., Van Der Voo, R., Doubrovine, P. V., Burke, K., Steinberger, B., Ashwal, L. D., Trønnes, R. G., Webb, S. J., & Bull, A. L. (2014). Deep mantle structure as a reference frame for movements in and on the Earth. Proceedings of the National Academy of Sciences of the United States of America, 111(24), 8735–8740. https://doi.org/10.1073/PNAS.1318135111 Trampert, J., Deschamps, F., Resovsky, J., & Yuen, D. (2004). Probabilistic tomography maps chemical heterogeneities throughout the lower mantle. Science, 306(5697), 853–856. https://doi.org/10.1126/SCIENCE.1101996/SUPPL_FILE/TRAMPERT.SOM.PDF Tucker, J. M., & Mukhopadhyay, S. (2014). Evidence for multiple magma ocean outgassing and atmospheric loss episodes from mantle noble gases. Earth and Planetary Science Letters, 393, 254–265. https://doi.org/10.1016/j.epsl.2014.02.050 Wang, Y., & Wen, L. (2007). Geometry and P and S velocity structure of the “African Anomaly.” Journal of Geophysical Research: Solid Earth, 112(5). https://doi.org/10.1029/2006JB004483 Wasserburg, G. J., & DePaolo, D. J. (1979). Models of earth structure inferred from neodymium and strontium isotopic abundances. Proceedings of the National Academy of Sciences, 76(8), 3594–3598. https://doi.org/10.1073/PNAS.76.8.3594 Wen, L., Silver, P., James, D., & Kuehnel, R. (2001). Seismic evidence for a thermo-chemical boundary at the base of the Earth’s mantle. Earth and Planetary Science Letters, 189(3–4), 141–153. https://doi.org/10.1016/S0012-821X(01)00365-X White, W. M. (2015). Probing the earth’s deep interior through geochemistry. Geochemical Perspectives, 4(2), 95–251. https://doi.org/10.7185/GEOCHEMPERSP.4.2 Williams, C. D., Li, M., McNamara, A. K., Garnero, E. J., & Van Soest, M. C. (2015). Episodic entrainment of deep primordial mantle material into ocean island basalts. Nature Communications, 6. https://doi.org/10.1038/NCOMMS9937 141 Yuan, Q., & Li, M. (2022). Instability of the African large low-shear-wave-velocity province due to its low intrinsic density. Nature Geoscience 2022 15:4, 15(4), 334–339. https://doi.org/10.1038/s41561-022-00908-3 Zhang, N., Zhong, S., Leng, W., & Li, Z. X. (2010). A model for the evolution of the Earth’s mantle structure since the Early Paleozoic. Journal of Geophysical Research: Solid Earth, 115(6). https://doi.org/10.1029/2009JB006896 Zhao, C., Garnero, E. J., McNamara, A. K., Schmerr, N., & Carlson, R. W. (2015). Seismic evidence for a chemically distinct thermochemical reservoir in Earth’s deep mantle beneath Hawaii. Earth and Planetary Science Letters, 426, 143–153. https://doi.org/10.1016/J.EPSL.2015.06.012 Zhong, S. (2006). Constraints on thermochemical convection of the mantle from plume heat flux, plume excess temperature, and upper mantle temperature. Journal of Geophysical Research: Solid Earth, 111(4), 4409. https://doi.org/10.1029/2005JB003972 142 APPENDIX Table 5.3 Effects of grain size on diffusion creep viscosity* Non-dimensional Diffusion Non-dimensional grain creep viscosity of dense size of dense material2 material1 1 500 1,000 2,000 5,000 1.00 12.01 15.85 20.91 30.17 * The table shows how the intrinsic diffusion creep viscosity of dense material (thermochemical pile and ULVZ material) is related to the grain size of it. The calculated values are based on the power-law relationship between grain-size and viscosity in diffusion creep as stated in McNamara et al., 2003: 𝜂𝑑𝑖𝑓𝑓 = 𝜇 𝐴′ 𝑑𝑖𝑓𝑓 ( 𝑑 𝑏 ) 𝑚 exp( 𝑔𝑑𝑖𝑓𝑓𝑇𝑚 𝑇𝑑𝑖𝑚 ), where 𝑇𝑚is the dimensional melting temperature, 𝑇𝑑𝑖𝑚 is the dimensional temperature, d is the grain size, 𝑔𝑑𝑖𝑓𝑓 is an empirical parameter related to activation free energy, 𝜇 is the reference rigidity, b is the reference Burgers vector, and 𝐴′ 𝑑𝑖𝑓𝑓 is a prefactor. m is a constant and is taken as 2.5 according to Karato and Wu, 1993. Note here that = 𝜂𝑑 𝜂𝑚 ( 𝑑𝑑 )2.5 , where 𝜂𝑑 and 𝜂𝑚 are the intrinsic viscosity of the dense material and the 𝑑𝑚 background mantle, and 𝑑𝑑 and 𝑑𝑚 are the grain size of the dense material and the background mantle. 1. Non-dimensional diffusion creep viscosity of dense material, taking the diffusion creep viscosity of the background mantle as 1. 2. Non-dimensional grain size, taking the grain size of the background mantle as 1.00. Note: Although grain size affects only the diffusion creep viscosity, because thermochemical piles are low stress regions and are dominated by diffusion creep14, it is a reasonable approximation to multiply the viscosity formula by a composition- dependent prefactor 𝜂𝑐(𝐶) in the model. 143 Figure 5.10 The logarithm of non-dimensional initial effective viscosity fields used for the cases in experiment sets C and D in this study. The intrinsic viscosity contrasts between thermochemical piles and the background mantle are none in case C1 and D1 (a), 100x in case C2 and D2 (b), 1,000x in case C3 and D3 (c), 2,000x in case C4 and D4 (d), and 5,000x in case C5 and D5 (e), respectively. Rayleigh number applied in these sets is 9.34 x 107. 144 Figure 5.11 Initial temperature and compositional conditions for Cases 21-24 (a) and the logarithm of non-dimensional initial effective viscosity fields used for Case 21 (b), Case 22 (c), Case 23 (d), and Case 24 (e). The intrinsic viscosity contrasts between thermochemical piles and the background mantle are none for (b) and (d), and 2,000x in (c) and (e). The Rayleigh number applied in these cases is 9.34 x 107. 145 Figure 5.12 Snapshots in a time sequence for Cases 21-22 with 1x and 2,000x of pile intrinsic viscosity, respectively. Time t is scaled to geologic time after incorporating hot anomalies and is listed on top of the snapshots. The figure on the left shows the initial thermal and compositional condition of the two cases. Each panel shows snapshots for the same case, while each column compares snapshots for the 2 cases at the same time. Color scales here have the same meaning as that in Figure 5.4. Rayleigh number applied in these sets is 9.34 x 107. 146 Figure 5.13 Snapshots in a time sequence for Cases 23 and 24 with 1x and 2,000x of pile intrinsic viscosity, respectively. Time t is scaled to geologic time after the incorporation of hot anomalies and is listed on top of the snapshots. The figure on the left shows the initial thermal and compositional condition of the two cases. Each panel shows snapshots for the same case, while each column compares snapshots for the 2 cases at the same time. Color scales here have the same meaning as that in Figure 4. The Rayleigh number applied in these sets is 9.34 x 107. 147 CHAPTER 6: CONCLUSIONS AND PERSPECTIVES This dissertation advanced the understanding of the nature, physical properties, dynamics, and thermochemical structures of the Earth’s deep mantle by performing geodynamic modeling. Chapter 1 introduced the background and importance of mantle convection, especially its lower boundary layer ---- the Earth’s lowermost mantle. It described limitations in the observations towards the Earth’s lowermost mantle. It demonstrated the importance and necessity of multi- disciplinary studies or perspectives for understanding the structures, origin, and dynamics of the Earth’s lower mantle. It emphasized the role geodynamic modeling plays in advancing the fundamental understanding of the Earth’s lowermost mantle. Chapter 2 illustrated the governing equations for mantle convection (equations of conservation of mass, momentum, and energy) and the numerical methods to solve these equations. It also discussed some model setups, such as model domain selections and formulation of thermochemical piles for investigations in Chapters 3-5. In Chapter 3, we performed 3D partial spherical calculations to investigate if the lowermost mantle flow patterns can reflect slab trench length and the viscosity contrast between slabs and the background mantle (slab strength). We found that strong and wide slabs induced unidirectional lowermost mantle flow patterns, while weak and narrow slabs induced more azimuthally divergent flow patterns. This difference will likely cause a characteristically different directionality of the crystal-preferred orientation that seismic anisotropy observations can detect. With more seismic anisotropy observations in the future, the results could provide more constraints on slab strength and trench length. Chapter 4 presented that the cross-sectional shapes of ULVZs are controlled by the viscous coupling between the interface of ULVZs and their surrounding materials. We performed high- 148 resolution 2.5 D spherical calculations for the investigation under the assumption that ULVZs are caused by ultra-dense material accumulations (UDMAs) and the LLSVPs are caused by thermochemical piles. We find that differential viscous coupling at the edge of thermochemical piles leads to differential stress acting on the two sides of UDMAs and causes asymmetrical UDMA shapes; uniform viscous coupling within the pile interior leads to uniform stress acting on the two sides of UDMAs and generates symmetrical UDMA shapes. For UDMAs at pile edges, if UDMAs and pile are intrinsically more viscous than the background mantle, the cross-sectional UDMA shapes are more symmetrical or exhibit the opposite asymmetry compared to the ones with similar intrinsic viscosity as the background mantle. Thus, in the future, if regional high-resolution seismic studies can resolve the cross-sectional shapes of ULVZs, we may be able to use the cross- sectional shapes of ULVZs to infer the rheological properties of the LLSVPs and ULVZs. Chapter 5 examined if and how the increased intrinsic viscosity of thermochemical piles affects the lateral mobility of thermochemical piles along the CMB and the stability of cross- sectional shapes of piles. Assuming the LLSVPs originate from dense thermochemical piles, this chapter aims to discuss this mechanism that can cause the LLSVPs to be stationary for up to a few hundred million years, under debate in the recent decade in the geoscience community. We performed 2.5 D thermochemical calculations to investigate how piles with varied intrinsic viscosities respond to changing upwelling and downwelling flow induced by thermal anomalies introduced into the models. We find that the increased intrinsic viscosity of piles will not cause them to be stationary along the CMB but will stabilize their cross-sectional morphology. Therefore, if the LLSVPs are indeed fixed, another mechanism must be found to fix them. The results of this dissertation open a few interesting questions to explore in the future. For example, 149 1) In Chapter 3, I explored slab deformation in the lower mantle in isochemical 3D models. It would be interesting to experiment with how slabs deform in the lower mantle and the fate of slabs in the lower mantle with the presence of piles. Furthermore, how will the interactions between slab and piles affect the lowermost mantle flow patterns near pile edges? Will this explain the seismic anisotropy features near the edges of LLSVPs? 2) In Chapter 4, we focused on UDMAs within or at edges of thermochemical piles, but it’s also interesting to think about ULVZs outside of LLSPVs. In the ULVZ detection map in a recent study (Yu & Garnero, 2018), it was shown that some ULVZs are outside of the LLSVP margins. What if those instinct small patches are actually pieces of continuous larger patches, just not resolved due to seismic data coverage? In our models, we sometimes see UDMAs within tiny piles. Can the ULVZs outside of the LLSVPs, if actually more continuous and larger than they appeared now, actually be within or be tiny piles distributed outside of the LLSVPs? 3) In the results of Chapter 5, we observed that plumes can consistently form outside of the piles with increased intrinsic viscosity. In conventional geodynamic models, we always see plumes forming on top of the piles because piles swiftly accommodate themselves to upwellings. However, if we consider oceanic island basalts (OIBs) to originate from deep mantle plumes rooted in the lowermost mantle, it would not be ideal to use conventional geodynamic models to explain the OIBs outside the LLSVP margins. The preliminary findings regarding plume source locations in this chapter lead to the question: would the increased intrinsic viscosity of piles be possible to explain 150 the location and geochemical signature of the OIBs outside of the present-day LLSVP margins? 4) In the results of Chapter 5, we observed that the entrainment of pile materials into the upwellings reduces as the pile intrinsic viscosity increases. However, we did not quantify this effect because it’s beyond the scope of this chapter. This preliminary result opens the question: Would this explain the variety of geochemical signatures of OIBs, i.e., some entrained more primordial materials as it’s formed while others did not? 5) Based on the conclusion of Chapter 5, the cross-sectional shape of the LLSVPs can be affected by their intrinsic viscosity. From seismic tomography models, the Pacific LLSVP has a smaller height than that of the Africa LLSVP (e.g., Y. He & Wen, 2009; Wang & Wen, 2007). Can this large-scale feature be explained with different intrinsic viscosity between them? What would it mean if the two LLSVPs have different intrinsic viscosity? Can they have different compositions or different evolution histories? 6) One hypothesis for the origin of thermochemical piles not focused in this dissertation is that they are caused by the accumulation of oceanic crust subducted to the CMB (e.g., M. Li et al., 2014). A recent study (Dannberg et al., 2017) incorporated grain size dependent rheology into isochemical models and found that slabs underwent lateral periodic buckling, unlike observed in models with constant grain size. In the model, slabs buckled so much that they can accumulate in the lower mantle for a time frame of a few hundred million years. It might be worthwhile to see if dense slabs can 151 accumulate more easily into thermochemical piles in grain size dependent thermochemical models. 7) To further (6) and combine methodology in Chapters 3 and 5, it is natural to think: What if the LLSVPs are caused by a combination of the two mechanisms, with oceanic crust subduction being a recharge for the pile accumulations? What geochemical signature will it lead to in the OIBs in this situation? How can we reconcile surface plate motions with the existence of the two LLSVPs? 152 REFERENCES Dannberg, J., Eilon, Z., Faul, U., Gassmöller, R., Moulik, P., & Myhill, R. (2017). The importance of grain size to mantle dynamics and seismological observations. Geochemistry, Geophysics, Geosystems, 18(8), 3034–3061. https://doi.org/10.1002/2017GC006944 He, Y., & Wen, L. (2009). Structural features and shear‐velocity structure of the “Pacific Anomaly” J. Geophys. Res, 114(2), 2309. https://doi.org/10.1029/2008JB005814 Li, M., McNamara, A. K., & Garnero, E. J. (2014). Chemical complexity of hotspots caused by cycling oceanic crust through mantle reservoirs. Nature Geoscience, 7(5), 366–370. https://doi.org/10.1038/ngeo2120 Wang, Y., & Wen, L. (2007). Geometry and P and S velocity structure of the “African Anomaly.” Journal of Geophysical Research: Solid Earth, 112(5). https://doi.org/10.1029/2006JB004483 Yu, S., & Garnero, E. J. (2018). Ultralow Velocity Zone Locations: A Global Assessment. Wiley Online Library, 19(2), 396–414. https://doi.org/10.1002/2017GC007281 153