University of Southern California Dissertations and Theses

Deep structure of Earth's tectonic regions

(USC Thesis Other)

Deep structure of Earth's tectonic regions

PDF

Download

Open document

Flip pages

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content 1 DEEP STRUCTURE OF EARTH’S TECTONIC REGIONS by Elizabeth M. Paulson A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY GEOLOGICAL SCIENCES AUGUST 2016 2 Table of Contents Table of Contents ............................................................................................................................ 2 Acknowledgements ......................................................................................................................... 6 Preface............................................................................................................................................. 7 1 Convergence Depths of Tectonic Regions from an Ensemble of Global Tomographic Models............................................................................................................................................. 8 1.1 Abstract ............................................................................................................................ 8 1.2 Introduction ...................................................................................................................... 9 1.3 Model Selection.............................................................................................................. 10 1.4 Projection Method .......................................................................................................... 14 1.5 Profile Variance Analysis............................................................................................... 16 1.6 Tectonic Correlations ..................................................................................................... 21 1.6.1 Significant regional features ................................................................................... 22 1.6.2 Variability among tomographic models.................................................................. 24 1.7 Estimation of Convergence Depths ................................................................................ 29 1.7.1 Apparent Convergence Depths ............................................................................... 30 1.7.2 Vertical-Smearing Bias ........................................................................................... 34 1.7.3 Smoothing Kernels.................................................................................................. 38 1.7.4 Bias-Corrected Lower Bounds on Convergence Depths ........................................ 43 1.7.5 Ensemble Lower Bounds on Convergence Depths ................................................. 46 3 1.8 Other Sources of Bias in Convergence-Depth Estimates ............................................... 49 1.8.1 Form of the Bias Correction ................................................................................... 50 1.8.2 Inter-Regional Correlations .................................................................................... 52 1.8.3 Tomographic Inversion Filters ................................................................................ 54 1.8.4 Tectonic Regionalization ........................................................................................ 56 1.9 Discussion ...................................................................................................................... 56 1.9.1 Intra-Oceanic Convergence Depths ........................................................................ 57 1.9.2 Continent-Ocean Convergence Depths ................................................................... 61 1.9.3 Correlated Mantle Flow: an Alternative Hypothesis? ............................................ 64 1.10 Conclusions .................................................................................................................... 66 2 Thermal Aging of the Oceanic Asthenosphere ..................................................................... 69 2.1 Abstract .......................................................................................................................... 69 2.2 Introduction .................................................................................................................... 70 2.3 Model Selection.............................................................................................................. 73 2.4 Projection Method .......................................................................................................... 76 2.5 Profile Variance Analysis............................................................................................... 78 2.6 Tectonic Correlations ..................................................................................................... 79 2.7 Estimation of convergence depths.................................................................................. 80 2.8 Travel Time Variance Analysis...................................................................................... 85 2.9 Tectonic TT Correlations ............................................................................................... 89 4 2.10 Travel Time by Age Profiles .......................................................................................... 91 2.11 Best-fitting HSC Models ................................................................................................ 95 2.12 Discussion .................................................................................................................... 100 2.13 Conclusions .................................................................................................................. 105 3 Testing the Isopycnic Hypothesis Using Garnet Lherzolite Xenoliths ............................... 107 3.1 Abstract ........................................................................................................................ 107 3.2 Introduction .................................................................................................................. 108 3.3 Data Set ........................................................................................................................ 112 3.4 Pressure and Temperature ............................................................................................ 114 3.5 Density ......................................................................................................................... 117 3.6 Ol Mg# as a Proxy for Density..................................................................................... 117 3.7 Estimated Geotherms ................................................................................................... 118 3.8 Oceanic Adiabat ........................................................................................................... 120 3.9 Test of Isopycnicity ...................................................................................................... 120 3.10 Theil Slope Estimation Method.................................................................................... 122 3.11 Discussion .................................................................................................................... 127 4 Appendices .......................................................................................................................... 130 4.1 Appendix A: Calculation of Profile Variance .............................................................. 130 4.2 Appendix B: Regional Profiles .................................................................................... 132 4.3 Appendix C: Bias and Variance of Regional Profiles .................................................. 138 5 4.4 Appendix D: A Bias Model for Vertical Smearing ..................................................... 141 4.5 Appendix E: Regional Radial Correlation Functions ................................................... 142 4.6 Appendix F: Bayesian Interpretation of Convergence Depth Probabilities ................. 144 4.7 Appendix G: Density vs. Depth for Thermobarometers .............................................. 147 List of Figures ............................................................................................................................. 151 References ................................................................................................................................... 159 6 Acknowledgements I would like to express my gratitude to my advisor, Professor Thomas H. Jordan for guiding me in my study and research at USC, for his unfailing patience and support, his assistance in improving my research, presentation and writing skills, for the opportunities he has given me to be a part of the larger scientific community, as well as his guidance through the qualification exam, dissertation writing, and dissertation defense process. Secondly, I would like to thank my dissertation and qualification committee members, Dr. Charles Sammis, Dr. Patrick Lynett, Dr. Thorsten Becker, Dr. Meghan Miller, and Dr. Lloyd Armstrong for the generous contribution of their time and expertise. I would also like to thank Cindy Waite for going above and beyond to help me handle registration issues, the degree process, and paperwork, John McRaney for his help with fellowship and payroll issues, John Yu for helping with computing resources, and Barbara Grubb for help with teaching assistantships and lab setups. Finally, I would like to thank the USC Provost’s PhD Fellowship as well as the Earth Sciences department for their generous support. 7 Preface The objective of this thesis is to study the structure of the lithosphere and underlying asthenosphere as well as the correlation of this structure with its overlying tectonic regions. We use several approaches to examining this structure. The first part involves examining the average aspherical velocity structure beneath the Earth’s tectonic regions and determining convergence depths below which the velocity profiles of these regions become indistinguishable from each other. The second part involves an in-depth examination of the Earth’s oceanic regions and their underlying thermal structure. We examine the convergence depths of the oceanic regions for a finer tectonic ocean age regionalization (OAR1) in order to determine related convergence depths, as well as average travel time models for each region in order to examine the thermal structure of that region while minimizing vertical smearing bias which is introduced when examining a region based simply on the underlying velocity structure. The third part uses garnet lherzolite xenoliths to statistically test whether we can reject the null hypothesis that the xenolith’s density changes with depth in accordance with the isopycnic hypothesis, which has implications regarding the structure of cratonic continental lithosphere and asthenosphere. All of these approaches allow us to examine the LAB hypothesis as it applies to both oceanic and continental regions, as well as comparing models of the thermal structure of the oceanic asthenosphere. 8 1 Convergence Depths of Tectonic Regions from an Ensemble of Global Tomographic Models 1.1 Abstract The aspherical variations in shear wave velocities (v S ) from 27 whole mantle tomographic models were projected onto the GTR1 global tectonic regionalization. A set of convergence depths was estimated, below which the v S profiles for the GTR1 regions become statistically indistinguishable from each other. These convergence depths were corrected for vertical smearing bias. At 90% confidence, the estimated lower bound on the convergence depth between A (young- < 25 Ma ) and B (intermediate- 25-100 Ma ) age oceans was found to be >152 km, while the lower bound on the convergence depth between C (old >100 Ma) and B (intermediate) age oceans was calculated to be >178 km. Regions P (platforms) & S (shields) were combined to define stable continental areas and regions B & C were combined to define mature oceans. The lower bound on the convergence depths between these regions was found to be >340 km. The oceanic convergence depths are much deeper than the G discontinuity, implying that more than 100 km of oceanic asthenosphere is advected by plate flow and cools according to a square- root of age relationship to >200 my. The stable continents, in addition to showing a deep convergence depth, have high radial correlation, which suggests that the convergence depth between stable continental areas and mature oceans defines the base of the kinematic plate. This is contrary to the view of the lithosphere-asthenosphere boundary as the kinematic base of the tectonic plates, and more consistent with a thick tectosphere decoupled from mantle flow by a low-viscosity layer immediately above the 410-km discontinuity. This study is an expansion on 9 an earlier work (Jordan and Paulson, 2013) , with a larger model list and updated calculated values, figures, and discussions. 1.2 Introduction The lithosphere-asthenosphere hypothesis (LAB) states that the plates (lithosphere) are decoupled from deeper mantle flow by a weak zone of lateral shearing (asthenosphere). Under the LAB hypothesis, the continental keels cannot be too deep, or they would impede the flow of the asthenospheric mantle below them. However, this assumption was challenged by early studies showing high velocity continental keels extending to depths of greater than 300 km under cratons (e.g.,(Jordan, 1975); (Woodhouse and Dziewoński, 1984); (Lerner-Lam and Jordan, 1987)). Recent regional studies have both supported (e.g.,(Grand, 1994); (Gaherty and Jordan, 1995); (Fouch et al., 2004); (Wittlinger and Farra, 2007)), and questioned ((Gung et al., 2003); (Fishwick et al., 2005); (Yuan and Romanowicz, 2010)) the deep keel assumption. The latter studies support a maximum depth for the continental keels of 200-250 km. The purpose of this study is to use global tomographic models to analyze, in a statistically rigorous way, the correlation between surface tectonic structure and the mantle beneath it. Global tomographic models are generally independent of geologic observations, so the correlation between them and known surficial tectonic structure should give reliable information about the relationship between surface tectonics and the underlying mantle structure. As there are many tomographic models, this study aims to evaluate an ensemble of models, analyze the similarities and differences which are observed among models, and search for globally consistent features related to Earth’s tectonic regions. 10 This study determines the convergence depths at which perturbations in S-wave velocity (dv S ) between continental and oceanic mantle regions are no longer distinguishable from each other. It also determines convergence depths between oceanic age regions. The oceanic convergence depths are thought to indicate the depth of the thermal boundary layer below oceanic age regions (Jordan and Paulson, 2013), while the stable continent- mature ocean convergence depth is thought to be the “decoupling zone” between Earths kinematic plates in the upper mantle and the lower mantle beneath continents. The results indicate that the LAB hypothesis should be amended. 1.3 Model Selection We considered only global tomographic models of whole-mantle shear-velocity (v S ) structure. To ensure adequate vertical resolution in the upper mantle, we focused on models that included travel times of multiply-reflected S waves with shallow turning points or equivalent higher-mode dispersion data. We did not consider models derived solely from fundamental-mode surface waves, which do not resolve small-scale structure in the critical 200-400 km depth interval. We also excluded models derived using explicit tectonic regionalizations (e.g. Khan et al., 2011) or strong mineralogical hypotheses (Cammarano and Romanowicz, 2007; Simmons et al., 2010), as well as those with variable lateral resolution designed for regional studies (Nettles and Dziewoński, 2008). Among the selected models (Table 1) are representatives from four centers of sustained tomographic research: the University of Texas (Grand et al., 1997; Grand, 2002; Simmons et al., 2007; Simmons et al., 2009), Harvard University (Su and Dziewoński, 1997; Ekström and Dziewoński, 1998; Gu et al., 2001; Ishii and Tromp, 2001; Kustowski et al., 2008), the Scripps Institution of Oceanography (Masters et al., 1996; Masters et al., 2000; Montelli et al., 2006; 11 Houser et al., 2008), and the University of California at Berkeley (Megnin and Romanowicz, 2000; Panning and Romanowicz, 2006; Lekić et al., 2010; Lekić and Romanowicz, 2011). We also included three models by Ritsema et al., ( 2000; 2004; 2010) and two models by Takeuchi (2007; 2012) and a single model by Auer et al., ( 2014). The model sequences from each group reflect increases in data coverage as well as evolving approaches to data measurement, model parameterization, and the type of wave theory employed to relate the data functionals to elastic structure. The models in Table 1 exhibit several levels of complexity in the representation of upper- mantle anisotropy. Some authors parameterized the elastic structure solely in terms of SH velocities and relied exclusively, or primarily, on measurements of SH-polarized waves: all six Texas models, SAW24B16 of Mégnin and Romanowicz (2000), and both Takeuchi models (2007; 2012). Others inverted data sensitive to both SH and SV velocities assuming isotropic perturbations to a radially anisotropic structure: all five Scripps models, all three Ritsema models, and MK12WM13 (Su and Dziewoński, 1997), SPRD6 (Ishii and Tromp, 2001), and S362D1 (Gu et al., 2001). A third set of authors inverted for lateral heterogeneity in the radial anisotropy: S20A (Ekström and Dziewoński, 1998), S362ANI and S362WMANI (Kustowski et al., 2008), SAW642AN (Panning and Romanowicz, 2006), SAW642ANb (Lekić et al., 2010), SEMum (Lekić and Romanowicz, 2011), and SAVANI (Auer et al., 2014). The latter set of tomographic studies suggest that there are interesting tectonic correlations in the global distribution of upper-mantle anisotropy (see also Montagner and Tanimoto, 1991; Nettles and Dziewoński, 2008). However, the global tomographic models still show significant inconsistencies in anisotropic structure, leading Kustowski et al. ( 2008) to conclude that “the determination of radial anisotropy in the mantle is still an ongoing experiment and only few 12 anisotropic anomalies may be considered robust.” With this guidance (see also Ferreira et al., 2010), we reduced the perturbations for the seven radially anisotropic 3D models to their approximate Voigt averages, v S = ( v SH + 2 v SV )/3. For all models, we subtracted out (if necessary) the spherical average of the isotropic shear-velocity perturbation at each depth z to ensure that the variation was aspherical: 1 4𝜋 ∫ v 𝑆 𝑆 1 ( 𝑧 , Ω) 𝑑 Ω = 0 (1-1) Here  = ( , ) is geographic position on the unit sphere S 1 ,  is colatitude, and d = sin  d d is an element of solid angle. Figure 1-1 Global tomographic models analyzed in this study, listed by research group. Maps are Mollweide equal- area projections showing the relative aspherical variation in isotropic or Voigt-averaged shear-wave speed at 150 km depth; positive variations are in blue and negative variations in red. See Table 1 for model descriptions and references. 13 Table 1-1 Tomographic Models Selected for Correlation Analysis Group Model Name Inner Scale S Type References Horizontal Vertical Texas GRAND 275 km 75-150 km isotropic (Grand et al., 1997) TXBW 275 km 75-150 km isotropic (Grand, 2002) TX2007 250 km 75-150 km isotropic (Simmons et al., 2007) TX2008 250 km 75-240 km isotropic (Simmons et al., 2009) TX2008J 250 km 75-240 km isotropic (Simmons et al., 2009) GyPSuM 275 km 75-240 km isotropic (Simmons et al., 2010) Harvard MK12WM13S 16º 220 km isotropic (Su and Dziewoński, 1997) S20A 10º 100-400 km voigt (Ekström and Dziewoński, 1998) SPRD6 29º 220 km isotropic (Ishii and Tromp, 2001) S362D1 11º 100-400 km isotropic (Gu et al., 2001) S362WMANI 11º 60-400 km voigt (Kustowski et al., 2008) S362ANI 11º 60-400 km voigt (Kustowski et al., 2008) Berkeley SAW24B16 8º 100-300 km isotropic (Megnin and Romanowicz, 2000) SAW642AN 8º 100-300 km voigt (Panning and Romanowicz, 2006) SAW642ANb 8º 100-300 km voigt (Lekić et al., 2010) SEMum 8º 100-300 km voigt (Lekić and Romanowicz, 2011) Scripps S16B30 12º 100 km isotropic (Masters et al., 1996) SB4L18 4º 100-200 km isotropic (Masters et al., 2000) SB10L18 10º 100-200 km isotropic (Masters et al., 2000) PRI-S05 300-800 km adaptive grid isotropic (Montelli et al., 2006) HMSL-S06 4º 100-200 km isotropic (Houser et al., 2008) Ritsema S20RTS 10º 60-300 km isotropic (Ritsema and Van Heijst, 2000) S20RTSb 10º 60-300 km isotropic (Ritsema et al., 2004) S40RTS 5º 60-300 km isotropic (Ritsema et al., 2010) Takeuchi SH18CE 11º 60-500 km isotropic (Takeuchi, 2007) SH18CEx 11º 60-500 km isotropic (Takeuchi, 2012) Auer SAVANI 1.25° -5° variable voigt (Auer et al., 2014) The model diversity is evident in Figure 1-1, which plots the aspherical variations v S at z = 150 km. The maps show a range of amplitude and lateral roughness, reflecting differences in the data sets, model parameterizations, and dampings used in the inversions. In particular, the lateral roughness is limited by the model’s horizontal inner scale  0 , the smallest angular dimension of the model parameterization. For models represented by a spherical-harmonic expansion to angular degree l max , the inner scale (in radians) is  0 =√ 4𝜋 ( 𝑙 max + 1) ⁄ . The values of  0 , listed 14 in Table 1, range from less than 3º for the geographic grids used in Texas models to almost 30º for SPRD6, for which l max = 6. 1.4 Projection Method The tomographic models were projected onto a 5º  5º representation of the GTR1 global tectonic map (Jordan, 1981b), which divided Earth’s surface into three oceanic and three continental domains (Figure 1-2). Oceanic crust was partitioned according to the square-root of its magmatic age, as inferred from marine magnetic anomalies: region A (0-5 Ma 1/2 ), region B (5-10 Ma 1/2 ), and region C (10-15 Ma 1/2 ). Although continents can also be regionalized according to tectonic age (e.g., Artemieva, 2006), the complex deformation and magmatic history of the continents often makes this age ambiguous. In GTR1, the continents were instead divided according to their generalized tectonic behavior during the Phanerozoic eon. This interval of nearly 600 million years comprised a complete Wilson cycle of drifting, supercontinent amalgamation, breakup, and redrifting, and thus constitutes an extended period for assessing continental stability. The Phanerozoic orogenic zones and magmatic belts (region Q) were separated from stable, but slightly subsided and thus thinly sedimented, Phanerozoic platforms (region P), which in turn were distinguished from the stable and less subsided Precambrian shields and platforms (region S). 15 Figure 1-2 The GTR1 global tectonic regionalization on a 5º  5º geographic grid (Jordan, 1981). Area of each region is percentage of the sphere. GTR1 was synthesized in the late 1970s, and much more is now known about the history of the continents, especially along the passive continental margins where information about crustal tectonics has been substantially improved (e.g., Kearey et al., 2009). However, GTR1 is still used in seismological studies (e.g., Becker et al., 2008; Nettles and Dziewoński, 2008; Rychert and Shearer, 2009; Dalton and Faul, 2010; Rychert and Shearer, 2011), and its antiquity assures its independence of current knowledge of mantle structure. Tests described later in this paper Color Symbol Description Area (%) A Young ocean 13 B Middle-age ocean 35 C Old ocean 13 Q Phanerozoic orogenic and magmatic zones 22 P Phanerozoic platforms 10 S Precambrian shields and platforms 7 16 (§7.4) demonstrate that our conclusions are robust with respect to alternative tectonic regionalizations. The projection of a geographic function F( ) onto GTR1 yields six regional averages: 𝑓 𝑥 = 1 𝐴 𝑋 ∫ 𝐹 ( Ω) 𝑑 Ω 𝑋 , X  {A, B, C, Q, P, S}, (1-2) where 𝐴 𝑋 = ∫ 𝑑 Ω 𝑋 is the area of the region X. The fractional areas A X /4 , listed in Figure 1-2, sum to unity. Because we wanted to preserve the lateral structure of the tomographic models, we limited our choices of F( ) to radial functionals; i.e., functionals of v S (z, ) at fixed geographic position . In this paper, we will be primarily concerned with the relative aspherical perturbation 𝛿 v S ( 𝑧 , Ω) /v S ( 𝑧 ) , where v S (z) is the spherically averaged shear velocity. Other examples of radial functionals are the average aspherical gradient and the aspherical variation in the vertical travel time within some layer of the mantle. 1.5 Profile Variance Analysis For any radial functional F( ) with regional averages {f X }, the statistical significance of the inter-regional variations can be assessed from the intra-regional variances, which we will now define. We suppose that the observed variations in each region sample a Gaussian stochastic process  X with a mean value  X =〈𝐹 ( Ω) 〉,   X, which is independent of , and we define the inter-regional covariance function 𝑉 𝑋𝑌 ( Ω, Ω′ )= 〈( 𝐹 ( Ω)− 𝛷 𝑋 ) ( 𝐹 ( Ω′ )− 𝛷 𝑌 ) 〉, Ω ∈ 𝑋 , Ω ′ ∈ 𝑌 . (1-3) Then, the regional average f X is a sample from a Gaussian random variable with mean  X and variance 17 𝜎 𝑋 2 ≡ 〈( 𝑓 𝑋 − 𝜙 𝑋 ) 2 〉 = 1 𝐴 𝑋 2 ∫ 𝑋 ∫ 𝑋 〈( 𝐹 ( Ω)− 𝜙 𝑋 ) ( 𝐹 ( Ω ′ )− 𝜙 𝑌 ) 〉 𝑑 Ω𝑑 Ω′ = 1 𝐴 𝑋 2 ∫ 𝑋 ∫ 𝑋 𝑉 𝑋𝑋 ( Ω, Ω ′ ) 𝑑 Ω𝑑 Ω′ (1-4) Equation (4) provides the standard error in the mean (s.e.m.) needed to test the significance of the regional variations in the radial functional F( ). For example, if  X and  Y are Gaussian and uncorrelated (V XY =0), the difference in two regional averages f X – f Y will have a standard error of √𝜎 𝑋 2 + 𝜎 𝑌 2 , and the null hypothesis H 0 :  X =  Y can be rejected in favor of the alternative H 1 :  X >  Y at the (1 – ) confidence level if f X – f Y > ( ) ,√𝜎 𝑋 2 + 𝜎 𝑌 2 (1-5) where ( ) delimits the tail of area  for a normal distribution. At the 95% confidence level,  = 0.05 and  = 1.64. This one-sided significance test is appropriate when considering, say, the maximum depth of high-v S cratonic keels. A closely related concept is the statistical p-value, defined as the probability of obtaining a regional difference at least as large as the observed value f X – f Y , assuming that H 0 is true. The p-value for a Gaussian distribution can be expressed in terms of the error function: .𝑝 𝑋 −𝑌 = 𝑃 ( 𝛷 𝑋 − 𝛷 𝑌 > 𝑓 𝑋 − 𝑓 𝑌 |𝐻 0 ) = 1 2 [1 − 𝑒𝑟𝑓 ( 𝑓 𝑋 −𝑓 𝑌 √2( 𝜎 𝑋 2 +𝜎 𝑌 2 ) )] (1-6) In cases where the ordering is not determined by the alternative hypothesis, such as when examining the significance of regional variations in the lower mantle, we evaluate the null hypothesis H 0 :  X =  Y against the alternative H 1 :  X   Y using a two-sided test: 18 | f X – f Y | > ( /2) √𝜎 𝑋 2 + 𝜎 𝑌 2 . (1-7) Then, the 95% confidence interval corresponds to  = 1.96. To estimate the intra-regional variances, we adopt the Fisher-von Mises covariance function, which is commonly used in the statistical analysis of data on a sphere (e.g., Watson, 1982), 𝑉 𝑋𝑋 = 𝘀 𝑋 2 e 𝜅 𝑋 ( 𝑐𝑜𝑠𝛹 −1) (1-8) This covariance function is stationary and isotropic, depending on  and ' only through their angular separation, cos 𝛹 = Ω · Ω′= cos𝜃 cos𝜃 ′ + sin θ sin θ’ cos( φ -φ’) . (1-9) It has two parameters, the single-point variance 𝘀 𝑋 2 ≡ 〈( 𝐹 ( Ω)− 𝜙 𝑋 ) 2 〉 and a dispersion parameter  X , which we assume are constant throughout the region X. When  << 1, the cosine can be approximated as a quadratic, cos ≈ 1 −  2 2 ⁄ , and the covariance function reduces to a Gaussian with an angular half-width 𝜆 𝑋 = 𝜅 𝑋 −1/2 , 𝑉 𝑋𝑋 ( )= 𝘀 𝑋 2 e −  2 2𝜆 𝑋 2 ⁄ (1-10) When the region is the whole sphere, S 1 , the integral of the Fischer-von Mises covariance (8) reduces to 𝜎 2 = 𝘀 2 ( 1 − e −2𝜅 ) /2𝜅 . Thus,  =  when  = 0 (completely correlated process), whereas   0 as    (completely independent process). Algorithms for estimating the correlation angle  X (z) and s.e.m.  X (z) from discretized values of a radial functional are described in Appendix A. Generally speaking, the inter-regional variations in  X (z) are fairly small and not significantly different from the mean correlation angle obtained by averaging over all tectonic regions. is the horizontal analog of the radial l(z) l(z) 19 correlation length (z), the half-width of a model’s radial correlation function (Puster and Jordan, 1997). Radial correlation functions have been widely applied to assess tomographic structures and their implications for mantle convection (Jordan et al., 1993; Puster and Jordan, 1994; Puster and Jordan, 1997; Gu et al., 2001; Becker and Boschi, 2002; Panning et al., 2010). In this paper, we use the radial correlation lengths for individual tectonic regions,  X (z), to constrain the inversion bias in the regional profiles. For most of the tomographic models considered here, the depth variations of and are similar, as would be expected for heterogeneities of nearly constant aspect ratio. Minimum values of tend to occur in the upper part of the lower mantle, where the aspherical variations in shear velocities are weakest, or near the surface, where they are strongest. Most models show an increase in in the lowermost mantle, and some, such as S362D1 and S362ANI, give high values of in the transition zone, where the aspect ratio of the heterogeneity, / , becomes quite large, owing to model parameterizations that are decorrelated at the 660 km discontinuity (see Figure 1-9). r l r l l l l r 20 Figure 1-3 Minimum value of the horizontal correlation half-width 𝜆̅ min plotted against horizontal inner scale  0 for the 27 tomographic models in Table 1. Scaling relation 𝜆̅ min /  0 = ½ (dashed line) is set by the lower limit of model resolution. The minimum half-widths 𝜆 ̅ min are plotted against the model inner scales  0 in Figure 1-3. The calculated values are bounded by the scaling relation 𝜆 ̅ min =  0 /2, which is set by the lower limit of model resolution. The distance of a model above this line is indicative of the smoothing imposed in the inversion to suppress short-wavelength aspherical variations. Except for models with the largest inner scales (SPRD6 and MK12WM13), the horizontal correlation angles are small compared to the characteristic dimensions of the regions, which implies that the cross- correlation between two different regions is weak; i.e., 21 1 𝐴 𝑋 𝐴 𝑌 ∫ 𝑋 ∫ 𝑌 𝑉 𝑋𝑌 ( ΩΩ′ ) 𝑑 Ω𝑑 Ω ′ ≪ 𝜎 𝑋 2 (X ≠ Y). (1-11) Consequently, inequalities (5) and (7) should provide good statistical tests for the inter-regional differences in the aspherical averages. 1.6 Tectonic Correlations Identical correlation analyses were performed on the 27 tomographic models in Table 1-1. In Appendix B, we display the whole-mantle regional profiles of the relative aspherical variations v S /v S for all of the models. The profiles for each model are banded with their 1  errors, computed by the method described in Appendix A. We began our ensemble analysis by selecting a representative model, or “exemplar,” from each of the five multi-model groups as well as Auer’s SAVANI model. Regional profiles of the exemplars are plotted in Figure 1-4. 22 Figure 1-4 Profiles of relative aspherical variations v S /v S obtained by projecting the group exemplars onto the six tectonic regions of GTR1 (X = A, B, C, Q, P, S). The profiles are plotted as a function of depth from 50 km to the core-mantle boundary. Profiles for each region (dark lines) were calculated from equation (A2), and the 1  uncertainties (shaded bands) were calculated from equation (A3). Color code is the same as in Figure 1-2. Similar plots for the other 16 tomographic models considered here are given in Appendix B. 1.6.1 Significant regional features We use “significant regional features” to describe attributes of the regional profiles that can be statistically justified by the variance analysis for each and every exemplar. In other words, a regional feature must be common to all five of the exemplar models for it to be a significant feature. We applied the two-sided test (7) with  = 0.05 to judge whether two regional averages from the same exemplar are significantly different. Among the regional features that stand out in Figure 1-4 are the following: 23 1. The aspherical variations of shear velocity in the uppermost mantle correlate strongly with GTR1. They are most negative beneath young oceans (region A) and are most positive beneath stable continents (regions P and S). The latter are significantly higher than the shear velocities of region Q. Overall, the inter-regional variations observed in the uppermost mantle are substantially greater than the uncertainties of the regional averages deduced from the intra-regional variations. 2. The exception to the last statement is the similarity between Precambrian shields and platforms (S) and Phanerozoic platforms (P). Above 200 km, the region S averages are slightly higher than those of region P, but the differences are not significant at the 95% confidence level. In fact, there are no statistically significant differences between the P and S profiles anywhere in the mantle for any of the exemplars. 3. Shear velocities in the uppermost mantle beneath ocean basins increase systematically with crustal age in all exemplars. Each model exhibits an “apparent oceanic convergence depth”, below which the three oceanic regions (A, B, C) show no statistically significant variations. The apparent convergence depths of regions A and B are not significantly different from those of regions B and C. 4. There are statistically significant variations between continents and oceans at depths greater than the oceanic convergence depth. Each model is characterized by an “apparent continental convergence depth” at which the profiles for the stable-continent regions (P and S) merge with the oceanic profiles (A, B, and C) and below which these tectonic regions show no statistically significant differences. Among the exemplars, the apparent continental convergence depths exceed the apparent oceanic convergence depths by more than 100 km. 24 5. Throughout most of the lower mantle, the tectonic profiles generally hover within a few tenths of a percent of the zero line. The fluctuations are not statistically significant, consistent with an aspherical distribution of heterogeneity that is completely uncorrelated with GTR1. 6. The dispersion among the tectonic profiles broadens at the base of the mantle, showing a weak, but in some cases statistically significant, correlation with the surface distribution of continents. In the D" layer just above the core-mantle boundary (CMB), the average v S in oceanic regions tends to be lower than in continental regions. 1.6.2 Variability among tomographic models These significant regional features are shared by all six models in the Texas group (Figure 1-5 and B-1). The Texas models were derived from the travel times of SH-polarized phases in the interval between the direct S wave and the Love wave, including sS, ScS n , sScS n , and the multiply-reflected S waves. SKS and SKKS times were added to enhance lower-mantle resolution, but the regional profiles are SH-dominated throughout the most of the mantle. The waveform measurement techniques, model parameterizations (all isotropic), and inversion methods have been fairly stable during the Texas iterations, so that any structural changes mainly reflect improvements in the global sampling of SH travel times. The dataset has increased from 25,000 travel times, inverted by Grand (2002) for TXBW, to more than 46,000 travel times, inverted by Simmons et al. (2009) for TX2008, which we chose as the Texas exemplar. 25 Figure 1-5 GTR1 projections of the relative aspherical variations v S /v S for the Texas models, banded with their 1  uncertainties. Same layout as Figure B-1, but truncated at a depth of 700 km to highlight upper mantle structure. Dashed line at 300 km approximates the oceanic. The regional features are robust with respect to the improvements in data coverage. The more recent Texas models show slightly amplified fluctuations in the lower-mantle profiles relative to Grand et al.’s (1997) original model and small variations in amplitude and extent of tectonically correlated D" structure, but these minor features are not statistically significant at the 95% confidence level. The four models are consistent with apparent oceanic convergence depths near 300 km and apparent continental convergence depths around 400-500 km (Figure 1-5). The profiles for regions P and S are statistically indistinguishable at the 95% confidence level. 26 Region Q is significantly slower than regions P and S, but the apparent convergence depths for all three continental regions are similar. Five of the Harvard models are compatible with the regional features of §5.1 (Figure B-2). This agreement provides additional tests of robustness, because the Harvard datasets are larger and more diverse than the Texas datasets, incorporating measurements of P-SV waves, including Rayleigh surface waves, as well as SH waves, and the model parameterizations span a considerable range of inner scales. The mid-mantle profiles for the lower-resolution models SPRD6 (l max = 6) and MK12WM13 (l max = 12) are more variable than the more recent S362 models, but the profile uncertainties are correspondingly larger, primarily because the correlation angles increase with inner scale (Figure 1-3). Hence, the inter-regional differences in the lower mantle can be explained as fluctuations due to the intra-regional variance. The regional profiles of S362ANI, in which the aspherical variations were allowed to be radially anisotropic in the uppermost mantle, have features similar to those of its isotropic predecessor, S362D1 (Figure B-2), which indicates robustness with respect to the parameterization of upper-mantle anisotropy. S362ANI, chosen as the Harvard exemplar, was derived from more than 138,000 shear-wave travel times, 290,000 fundamental-mode phase anomalies, both Rayleigh and Love, and 47,000 long-period waveforms (Kustowski et al., 2008). The regional features in this mature model are compatible with those in less well-constrained and more simply parameterized models (e.g., MK12WM13 of Su and Dziewoński, 1997). This and other comparisons among groups of models suggest that increases in the data coverage are unlikely to alter the large-scale regional features considered here. One of the Harvard models, S20A, shows statistically significant fluctuations in the transition zone and upper part of the lower mantle that are incompatible with the other models (Figure B- 27 2). S20A was an early attempt to invert for aspherical variations in radial anisotropy on a global scale (Ekström and Dziewoński, 1998). The harmonic character of the fluctuations, especially the large excursions in the profiles for regions P and S between 300 km and 1300 km, suggests some kind of inversion instability, perhaps associated with under-constrained anisotropic perturbations. In any case, no such fluctuations are observed in the similarly parameterized and better constrained anisotropic model S362ANI. We have therefore labeled S20A as “anomalous” in Table 2. Four of the five Scripps models generally conform with the regional features of §5.1, except in the lower mantle, where the stable-continent profiles show systematically higher velocities than the oceanic profiles (Figure B-3). To evaluate the significance of these deep regional differences, we combined regions P and S into a “stable continent” profile (abbreviated PS) and combined regions B and C into a “mature ocean” profile (abbreviated BC); we then tested the difference between these profiles using inequality (5). The PS-BC differences in the lower mantle were not statistically significant at the 95% confidence level for S16B30 or SB4L18 (the Scripps exemplar), but they were for HMSL-S06, at depths 900 < z < 1400 km, and for SB10L18, at depths 700 < z < 1700 km. The PRI-S05 model lacks surface wave information and has poor resolution in the uppermost 400 km of the Earth’s mantle, and is considered anomalous (A) for the purposes of this study. The datasets used by the Scripps group included SS travel times but only those with lower-mantle turning points and not those of the higher-order multiple- S reflections or equivalent higher modes. Owing to this data deficiency, the Scripps inversions may be less capable than those of Texas and Harvard of separating upper and lower mantle heterogeneity. 28 A similar lower-mantle anomaly characterizes the anisotropic Berkeley model SAW642AN, which shows significantly positive PS-BC differences over the depth range of 900-2000 km (Figure B-4). The PS-BC differences in this interval are much suppressed in a revised version of the model, SAW642ANb (Lekić et al., 2010), which we chose as the Berkeley exemplar, and are absent (i.e., not statistically significant) in the other three Berkeley models. The latest Berkeley model, SEMum, is also anomalous, though at shallower depths, including A-C differences that swing between positive and negative values from the top of the mid-mantle transition zone to 1000 km. SEMum was derived by the inversion of low-frequency seismic waveforms and group velocity dispersion maps using the spectral element method to solve the forward problem (Lekić and Romanowicz, 2011). However, mantle structure was modified only above a depth of 1000 km; deeper structure was fixed to that of SAW24B16 (Megnin and Romanowicz, 2000). The SEMum profiles are identical with the starting model below 1000 km, as expected, but inconsistent above this depth (Figure B-4). The regional features of SAW24B16 are compatible with the exemplar set in §5.1, whereas those of SEMum are anomalous. The regional profiles of the three Ritsema models, S20RTS, S20RTSb and S40RTS, are very similar (Figure B-5), even though the latter was constrained by a much larger data set and parameterized with l max = 40 rather than just 20. Both models have regional features compatible with the §5.1 list, except that the PS-BC differential profile for S40RTS is significantly positive from the uppermost mantle to almost 1300 km. The apparent convergence depth is less than 500 km for S20RTS, which we selected as the exemplar. This comparison primarily reflects the much finer inner scale of S40RTS, which lowers its intra-regional correlation angles by a factor of two, decreasing the intra-regional variance in the lower-mantle. At issue is the degree to which regional variations in the upper mantle may have been “vertically smeared” into the lower mantle 29 owing to inadequate resolution of vertical structure. We use a model of vertical smearing to obtain bias-corrected estimates of convergence depths, as described in the next section. The SH18CE and SH18CEX models, which were based on SH waveform inversion using a Born approximation (Takeuchi, 2007; Takeuchi, 2012), are compatible with the group exemplars, although the amplitudes of their regional profiles are substantially lower and their convergence points are somewhat deeper (Figure 1-4). The ensemble comparison suggests that the upper-mantle anomalies in SH18CE and SH18CEx have been underestimated and vertically smeared to greater depths. Auer’s recent SAVANI model which relies on ~10,000,000 surface wave phase delays from fundamental modes up to the sixth overtone between 25 and 350s, as well as ~500,000 cross-correlation travel times of major body wave phases is compatible with the other models and is considered an exemplar (E). Table 1-2 separates the tomographic models into the three classes based on the ensemble of correlation analyses: the six exemplars (E) from which we derived the significant regional features of §5.1, fifteen models compatible with these regional features (C), and six models judged from spurious lower-mantle correlations to be anomalous (A). 1.7 Estimation of Convergence Depths All of the tomographic models in Table 1-1 show strong regional variations in f = v S /v S that generally agree with plate-tectonic expectations for the boundary-layer structure of the uppermost mantle (Figures B-1 to B-5). The regional profiles for most of the models converge monotonically to the spherically averaged Earth above the 660-km discontinuity, and the deeper fluctuations are consistent with lower-mantle heterogeneities dominated by convective structures uncorrelated with surface tectonics. If this behavior applies to the real Earth, then any two 30 tectonic regions X and Y will share a convergence depth somewhere in the upper mantle, below which the differences between the regional averages, f X-Y  f X – f Y , can be explained by a combination of modeling errors from the inversion process, which include the vertical-smearing bias and statistical fluctuations in the projections. The convergence of the regional profiles may be asymptotic rather than abrupt, so the problem of estimating convergence depths is asymmetric. If the shear-velocity structure of oceanic regions were governed by simple conductive cooling of the upper mantle, for example, the true convergence depths for young and intermediate-age ocean floor, , would not be well-defined, owing to the asymptotic convergence of the geotherms. However, we can use the depth ranges where the profile differences are significantly different from zero to place lower bounds on the convergence depths. 1.7.1 Apparent Convergence Depths We focus on the estimation of three convergence depths: , between young and intermediate-age ocean basins; , between old and intermediate-age ocean basins; and , between stable continent and mature ocean basin. The latter is estimated by lumping P and S into a single stable continent region (PS) and B and C into a mature-ocean region (BC), consistent with the ensemble analysis, which indicates that there are no significant differences between these regional pairs below the uppermost mantle. The A-B and C-B differential profiles for the six exemplar models are shown in Figure 1-6, and the PS-BC profiles in Figure 1-7. ˆ z X-Y ˆ z A-B ˆ z A-B ˆ z C-B ˆ z PS-BC 31 Figure 1-6 Differential profiles for oceanic region A relative to B (yellow curves) and oceanic region C relative to B (red curves) for the group exemplars. The curves are banded with their two-sided 95% confidence intervals. The model-specific A-B and C-B convergence depths, listed for all models in Table 1-2, are the minimum depths where 95% confidence bands intersect zero (orange line). 32 Figure 1-7 Differential shear-velocity profiles for the stable-continent region PS relative to the mature-ocean region BC (blue curves) for the group exemplars. The curves are banded with their two-sided 95% confidence intervals. The model-specific PS-BC convergence depths, listed for all models in Table 2, are the minimum depths where the 95% confidence bands intersect zero (red line). Above the apparent convergence depths, f A-B is strictly negative, whereas f C-B is strictly positive, as expected for an oceanic thermal boundary layer that cools with age. f PS-BC is also strictly positive, consistent with a cold, isopycnic tectosphere beneath the stable continents. A one-sided statistical test can thus be used to sharpen our terminology: we define the “apparent convergence depth at confidence level (1 – )” to be the minimum depth at which the null hypothesis H 0 : 𝑓 𝑋 −𝑌 = 0 can be rejected in favor of the alternative H 1 : 𝑓 𝑋 −𝑌 > 0 using inequality (5). The apparent convergence depths at the 97.5% confidence level are listed for the tomographic models in Table 2. 33 Table 1-2 Properties of Tomographic Models Group Model Name Class* Apparent Convergence Depths (km)** PS- BC A-B C-B Texas GRAND C 830 634 241 TXBW C 472 287 250 TX2007 C 466 286 270 TX2008 E 466 286 270 TX2008J C 430 286 243 GyPSuM C 448 282 246 Harvard MK12WM13 C 468 228 246 S20A A 482 209 232 SPRD6 C 462 202 223 S362D1 C 485 196 203 S362WMANI C 398 195 259 S362ANI E 428 192 253 Scripps S16B30 C 541 264 332 SB4L18 E 405 217 212 SB10L18 A 324 215 207 PRI-S05 A 1163 937 883 HMSL-S06 A 389 222 209 Berkeley SAW24B16 C 498 306 295 SAW642AN A 329 236 216 SAW642ANb E 374 244 271 SEMum A 851 192 250 Ritsema R20RTS E 475 179 187 S20RTSb C 485 190 198 R40RTS C 1273 220 234 Takeuchi SH18CE C 859 389 361 SH18CEx C 850 400 405 Auer SAVANI E 369 190 215 * Exemplar (E), compatible (C), and anomalous (A). ** Estimates based on a one-sided test ( α = 0.025) of profile differences uncorrected for bias. 34 1.7.2 Vertical-Smearing Bias The apparent depth criterion accounts for the statistical fluctuations in the observed values of 𝑓 𝑋 −𝑌 , but it completely ignores the bias caused by the vertical smearing of the inversion filter. To correct for vertical smearing, we assume the errors in the differential profiles 𝑓 𝑋 −𝑌 sample a Gaussian process with a depth-dependent bias, 𝛽 𝑋 −𝑌 ( 𝑧 )≡ 〈𝛷 𝑋 −𝑌 ( 𝑧 )− 𝑓 ̂ 𝑋 −𝑌 ( 𝑧 ) 〉 = 𝜙 𝑋 −𝑌 ( 𝑧 )− 𝑓 ̂ 𝑋 −𝑌 ( 𝑧 ) , (1-12) and a variance 𝜎 𝑋 −𝑌 ( 𝑧 )= √𝜎 𝑋 2 ( 𝑧 )+ 𝜎 𝑌 2 ( 𝑧 ) . Here 𝑓 ̂ 𝑋 −𝑌 is the true differential profile which, by definition, is zero below the actual convergence depth 𝑧 ̂ 𝑋 −𝑌 . Because the tectonic regions are large compared to correlation angles  X , the inversion process can be approximated as a linear low-pass filter of the differential profile (Appendix C), 𝜙 𝑋 −𝑌 ( 𝑧 )≈ ∫𝑔 𝑋 −𝑌 ( 𝑧 , 𝑧 ′ ) 𝑓 ̂ 𝑋 −𝑌 ( 𝑧 ′ ) 𝑑𝑧 ′ . (1-13) This deterministic equation states that the expected value of the observed differential regional profile, 𝜙 𝑋 −𝑌 , is the true differential profile 𝑓 ̂ 𝑋 −𝑌 filtered by a differential smoothing kernel g X-Y . In Appendix C, we show that g X-Y is a “bi-regional” average of the inversion filter (see equations C3 and C10a). To estimate the vertical smearing bias, we make three assumptions: (i) The upper part of the smoothing kernel {g X-Y (z, z'): z' < z} can be approximated as a Gaussian of half-width (z). (ii) g X-Y (z, z') is unimodular in z and therefore preserves the integral of 𝑓 ̂ 𝑋 −𝑌 ; i.e.,  ̂ 𝑋 −𝑌 ≡ ∫𝑓 ̂ 𝑋 −𝑌 ( 𝑧 ) 𝑑𝑧 = ∫∫𝑔 𝑋 −𝑌 ( 𝑧 , 𝑧 ′ ) 𝑓 ̂ 𝑋 −𝑌 ( 𝑧 ′ )𝑑𝑧𝑑 𝑧 ′ =∫𝛷 𝑋 −𝑌 ( 𝑧 ) 𝑑𝑧 . (1-14) (iii)The differential profile 𝑓 ̂ 𝑋 −𝑌 increases (or decreases) linearly with height above the convergence depth 𝑧 ̂ 𝑋 −𝑌 . 35 The first two statements are concordant with the basic principles for constructing inverse filters (Backus and Gilbert, 1968). In a well-conditioned inversion, the smoothing kernel 𝑔 𝑋 −𝑌 ( 𝑧 , 𝑧 ′ ) will be localized in a Gaussian-like peak, as expressed in (i). Tomographic data sets typically provide good integral constraints on regional upper-mantle structure, such as the vertical shear- wave travel times; hence, the kernel should also be nearly unimodular in z and preserve 𝛹 ̂ 𝑋 −𝑌 , as expressed in (ii). Figure 1-8 schematically illustrates the vertical smearing of one-sided, linear differential profile by a unimodular Gaussian filter. The results are qualitatively consistent with the shapes of the observed differential profiles in the vicinity of the apparent convergence depths (Figure 1-6 and Figure 1-7), supporting approximation (iii). 36 Figure 1-8 Schema for the bias calculation using the linear model of equation (17). Left panel: The bi-regional inversion kernel 𝑔 𝑋 −𝑌 centered at the actual convergence depth 𝑧 ̂ 𝑋 −𝑌 is characterized as a Gaussian with an upper half-width 𝛾 ̂ 𝑋 −𝑌 . Right panel: Vertical smearing of the actual structural difference 𝑓 ̂ 𝑋 −𝑌 (black line) results in an observed difference 𝑓 𝑋 −𝑌 (red line) with bias 𝛽 ̂ 𝑋 −𝑌 ≈ 0.4 𝑓 ̂ 𝑋 −𝑌 ( 𝑧 ̂ 𝑋 −𝑌 − 𝛾 ̂ ) . The bias model in equation (18) assumes that the inversion process conserves the model integral; i.e., that the areas under the red and black curves are the same. The apparent convergence depth at the (1 – ) confidence level is shown as the red dot. In Appendix D, we show that, under these assumptions, the bias evaluated at the convergence depth is 𝛽 ̂ 𝑋 −𝑌 = √ 2 𝜋  ̂ 𝑋 −𝑌 𝛾 ̂ ( 𝑧 ̂ 𝑋 −𝑌 − 𝑧 0 ) 2 , (1-15) 37 where z 0 is the minimum depth of integration and 𝛾 ̂ is the upper half-width of the kernel 𝑔 𝑋 −𝑌 ( 𝑧 𝑋 −𝑌 , 𝑧 ′ ) . For each depth z, we consider the null hypothesis H 0 ': z = 𝑧 𝑋 −𝑌 . The p-value for this hypothesis is the bias-corrected version of equation (8): 𝑝 𝑋 −𝑌 ( 𝑧 )= 1 2 [1 − erf ( 𝑓 𝑋 −𝑌 ( 𝑧 )− 𝑏 𝑋 −𝑌 ( 𝑧 ) √2𝜎 𝑋 −𝑌 ) ], (1-16) where b X-Y(z) is the bias conditional on H 0 '; i.e., obtained by substituting z for 𝑧 ̂ 𝑋 −𝑌 in equation (1-15. By construction, p X-Y (z) is the probability that a sample of the inversion process would equal or exceed the observed value f X-Y (z) if z were equal to 𝑧 ̂ 𝑋 −𝑌 . For the exemplars, as well as most of the other tomographic models, f X-Y is a monotonic function of z above the apparent convergence depth, and its magnitude decreases with depth more rapidly than either the conditional bias b X-Y or standard error 𝜎 𝑋 −𝑌 , implying that p X-Y (z) increases with depth. Monotonicity ensures that, if H 0 ' can be rejected in favor of the alternative, H 1 ': z < 𝑧 ̂ 𝑋 −𝑌 , then the augmented null hypothesis H 0 ": z ≥ 𝑧 ̂ 𝑋 −𝑌 can also be rejected in favor of its complement H 1 " = H 1 '. Therefore, 1 – p X-Y (z) measures our confidence that the trial depth z is less than the actual convergence depth 𝑧 ̂ 𝑋 −𝑌 ; i.e., that z is a lower bound on 𝑧 ̂ 𝑋 −𝑌 . These relationships are summarized in the one-sided convergence depth test: If p X-Y(z) is smaller than the critical value , then we reject the null hypothesis, H 0" : z ≥ 𝑧 ̂ 𝑋 −𝑌 , in favor of the alternative, H 1" : z < .𝑧 ̂ 𝑋 −𝑌 (1-17) To calculate the bias-corrected p-values, we use the tomographic estimate of (1-14),  𝑋 −𝑌 = ∫ 𝑓 𝑋 −𝑌 ( 𝑧 ) 𝑑𝑧 , 𝑧 1 𝑧 0 (1-18) 38 which yields the conditional bias formula, (1-19) We adopted z 0 = 50 km as the minimum depth of integration and chose the upper limit z 1 to be the apparent convergence depth at a relatively low confidence level (50-60%) in order to capture the full data constraint (Table 3), but  X-Y is insensitive to this cutoff. Table 1-3. Parameters of Bias Models Exemplar  (km)   0 *  c PS-BC A-B C-B TX2008 8.3 -3.6 2.7 31 0.12 S362ANI 8.5 -2.8 3.0 29 0.08 SB4L18 8.2 -2.6 2.5 32 0.07 SAW642ANb 7.6 -2.9 2.9 21 0.08 S20RTS 7.3 -2.3 2.3 20 0.13 SAVANI 8.0 -3.0 2.8 9 0.17 * Smoothing half-width at z 0 = 50 km 1.7.3 Smoothing Kernels Although the smoothing kernels for the regional differences g X-Y (z, z') can be directly calculated during the inversions, they were not available for the tomographic models analyzed here. Nevertheless, we can recover bounds on the smoothing half-widths (z) from the radial correlation properties of the models themselves (Appendix E). Because the global data coverage is fairly uniform, (z) should not vary much from region to region. Under the Gaussian approximation, the half-width of the regional radial correlation function is , where is the vertical scale lengths of the intra-regional heterogeneity (equation E9). In regions where the latter is small, (z) can therefore be approximated by . b X -Y (z) = 2 π Y X -Y g(z) (z-z 0 ) 2 r X 2 (z)= 2g 2 +h X 2 h X r X (z) / 2 39 We illustrate the method using the Harvard model S362ANI. The radial correlation functions for region BC are shown in Figure 1-9 and those for region PS in Figure 1-10. The correlation half-widths in the uppermost mantle are much larger for region PS than for region BC, indicating that structural coherence, rather than inversion smoothing, dominates beneath the stable continents. The strong radial correlation in the PS region, observed for all exemplars (Figure 1-10), is important evidence that the stable continents are vertically coherent structures comprising thick tectosphere. For the BC oceanic region, on the other hand, inversion smoothing appears to dominate: the radial correlation is weak, and the correlation lengths decrease towards the surface, as they should owing to the increasing resolving power of the surface-wave dispersion data at shallower depths. We therefore consider the approximation, 𝛾 ≈ 𝜌 BC ( 𝑧 ) /√ 2 , using as an estimator the distance from the diagonal to the = 0.78 contour in Figure 1-9a. e -1/ 2 40 Figure 1-9 Regional radial correlation functions for mature ocean basins (region BC), computed from relative aspherical variations v S /v S for the exemplar models. The radial correlation lengths in the uppermost mantle of region PS are significantly larger than those of region BC, indicating that the v S /v S structure beneath the stable continents is more vertically correlated. The strong decorrelation at 660 km for S362ANI is an artifact of the Harvard inversion scheme, which allows a discontinuity in the model perturbation at this depth (Kustowski et al., 2008). 41 Figure 1-10 Regional radial correlation functions for stable continent (region PS), computed from relative aspherical variations v S /v S for the exemplar models. The radial correlation lengths in the uppermost mantle of region PS are significantly larger than those of region BC (see Figure 1-9), indicating that the v S /v S structure beneath the stable continents is more vertically correlated. The relatively high radial correlation coefficient (~0.6) between PS structure in the uppermost mantle and PS structure near the core-mantle boundary is a consistent feature of the model ensemble. In Figure 1-11, we compare these values with the half-widths of the cubic B-splines used by Kustowski et al. (2008) to represent S362ANI structure. The resolving-power calculations of Gu et al. (2003), Panning and Romanowicz (2006), and Kustowski et al. (2008) indicate that the coefficients of the B-splines for the isotropic v S perturbations should be individually well resolved when averaged over the very large BC region (48% of Earth’s surface); hence, the 42 spline half-widths should be good approximations to the inversion filter half-widths. Although the radial correlation function shows statistical and perhaps structural fluctuations (e.g. in the mid-mantle transition zone), the overall agreement between the two estimates in Figure 1-11 indicates that the smoothing half-widths can be recovered from the radial correlation function. Figure 1-11 Smoothing half-widths (z) of the relative aspherical variations v S /v S for S362ANI. Three estimates are shown:  BC / 2 (solid curve), half-widths of the B-splines used by Kustowski et al. (2008) to represent S362ANI upper-mantle structure (points), and a linear fit to the latter with  0 = 33 km, b = 0.08 (thick dashed line). The thin dashed lines are the ±20% uncertainty bands on the linear model. The region-BC radial correlation functions for the other exemplars were fit to a linear model of the smoothing half-widths, (z) =  0 + c (z – z 0 ). (1-20) 43 The results for SB4L18, SAW642ANb, and S20RTS are similar in slope and offset to those of S362ANI (Table 3). The apparent smoothing half-width is largest for TX2008: 55 km at z = 250 km, compared with 45 km for S362ANI. TX2008 is the only exemplar derived by an inversion lacking surface-wave data, so its lower resolving power at shallow depths is not surprising. SAVANI shows a smoothing half width at 250 km similar to those of the other models (43 km at z = 250 km), however, the slope of the line is steeper than the other exemplars, SAVANI utilizes variable finely parameterized upper mantle and a depth–variable regularization scheme, choosing to down weigh multiply bouncing body wave phases. At shallow depths from 0-300 km, SAVANI vertical roughness was given 2x the weight of the horizontal damping terms, while below 300 km this weighting term increases linearly from 2 to 5. 1.7.4 Bias-Corrected Lower Bounds on Convergence Depths For the exemplars, the PS-BC conditional bias models decrease from 1-1.5% at 200 km to 0.2-0.3% at 450 km (Figure 1-12). The decrease comes from the integral constraint (18), which implies a shallower slope in the true differential profile, and therefore less bias, if 𝑧 ̂ PS−BC is deeper. For all exemplars, this effect is generally stronger than the bias increase due to the increased smoothing with depth implied by the positive values of the slope c in (1-20). 44 Figure 1-12 Conditional bias models for the PS-BC differential profiles, computed from equation (1-19) using the parameter values from Table 3. Dashed lines are ±20% uncertainty bands on the S362ANI bias model. Applying the bias corrections yields the 𝑝 𝑋 −𝑌 ( 𝑧 ) curves of Figure 1-13 and Table 1-4. The observed differential profiles usually decrease more rapidly with depth than the conditional bias, so that the p-values increase monotonically with depth. In accordance with the one-sided convergence depth test (17), 1 − 𝑝 𝑋 −𝑌 ( 𝑧 ) measures our confidence in asserting that a trial z is less than the X-Y convergence depth. In the case of S362ANI, 𝑝 PS−BC equals 0.1 at 358 km and 0.5 at 440 km (Table 4c). Therefore, we interpret 440 km to be a lower bound on 𝑧 ̂ PS−BC at the 50% confidence level (the median estimate), and 358 km to be a lower bound at the 90% 45 confidence Figure 1-13 Convergence depths of the differential profiles for A-B (left), C-B (middle), and PS-BC (right) as a function of p-value, computed from equation (17) using the bias parameters of Table 3. Dashed black line is the ensemble mean (Bayesian model average), which combines the aleatory and epistemic uncertainties. We interpret the p-value as the probability that a test depth z is greater than 𝑧 ̂ 𝑋 −𝑌 . The cumulative distributions 𝑝 𝑋 −𝑌 ( 𝑧 ) in Figure 1-13 can differ substantially from a Gaussian shape. At the high end, the p-values for some of the models flatten out as the differential profiles asymptotically approach zero below the apparent convergence depth. For differential profiles with strong overshoot (e.g., the A-B profile for SAW642ANb), the flattening is absent or occurs only at high p-values. At the low end, the p-value distributions have fatter tails than Gaussian, owing to the change in conditional bias with z (Figure 1-12). For some of the 46 A-B and C-B profiles, the growth in the conditional bias at shallow depth is sufficiently large that the p-value reaches a minimum. In these cases, we extrapolated the p-value to shallower depths as a constant equal to this minimum. Table 1-4 Convergence Depths Table 1-4a. B-A Convergence Depths (km) Exemplar p = 0.5 0.4 0.3 0.2 0.1 TX2008 277 271 264 255 243 S362ANI 183 179 173 164 143 SB4L18 248 237 224 204 - SAW642ANb 263 257 250 243 233 S20RTS 173 166 155 145 - SAVANI 196 194 191 188 184 Table 1-4b. C-B Convergence Depths (km) Exemplar p = 0.5 0.4 0.3 0.2 0.1 TX2008 287 277 266 251 233 S362ANI 282 270 254 242 223 SB4L18 249 240 228 213 - SAW642ANb 311 301 292 281 258 S20RTS 196 188 173 131 - SAVANI 227 223 218 212 202 Table 1-4c. PS-BC Convergence Depths (km) Exemplar p = 0.5 0.4 0.3 0.2 0.1 TX2008 440 434 427 419 408 S362ANI 440 419 397 382 358 SB4L18 422 411 397 380 350 SAW642ANb 385 378 370 362 350 S20RTS 466 457 446 433 412 SAVANI 342 338 334 329 322 1.7.5 Ensemble Lower Bounds on Convergence Depths The confidence levels for an individual exemplar in Table 1-4 are conditional on the tomographic model being “right.” In statistical jargon, they describe only the aleatory 47 (stochastic) variability intrinsic to regional averaging of model heterogeneities, which we have treated as random noise. They do not account for the epistemic (lack-of-knowledge) uncertainty associated with errors in model construction, parameterization, and inversion. The variation of the p X-Y curves among models in Figure 1-13 reflects such epistemic errors. Some contributions to epistemic uncertainty can be assessed through sensitivity analysis. For example, a 20% increase in the S362ANI bias model lowers the confidence level of ≥ 358 km from 90% to 78%, and a 20% decrease raises it to 96%. At the 90% confidence level, a 20% increase in bias drops the lower bound by 18 km, from 358 km to 340 km. In comparison, if the aleatory distribution were Gaussian, one standard deviation would be the depth difference between the p-values of 0.34 and 0.5, which is 35 km. Therefore, the 20% bias uncertainty is smaller than the aleatory uncertainty and also smaller than the variability among the exemplars. Sensitivity tests for other types of bias are described in §7. Thus far, our procedures for statistical inference, such as the use of the convergence-depth test (17), have conformed to the Neyman-Pearson paradigm of hypothesis testing and decision making (Cox and Hinkley, 1974). However, frequentist procedures of this sort are inadequate for assessing the epistemic error implied by the ensemble variation (e.g., Goodman, 1999), so we turn to a simple form of Bayesian model averaging (Hoeting et al., 1999). We interpret the p- value for the model M n as the probability that a test depth z is greater than or equal to 𝑧 ̂ 𝑋 −𝑌 conditional on M n being correct, which we express as 𝑝 𝑋 −𝑌 ( 𝑧 |𝑀 𝑛 ) . A Bayesian analysis of the specific conditions required for this interpretation, which include equally likely priors on H 0 " and H 1 ", is given in Appendix F. We suppose the n th model has a probability  n of correctly representing the real Earth. If the model ensemble {M n : n = 1, … , N} were a complete and mutually exclusive set—i.e., if one ˆ z PS-BC 48 and only one of the models were a correct representation of the real Earth—then the sum of  n over the set must be unity, and the unconditional ensemble probability of z ≥ is just the Bayesian model average (BMA), 𝑝 ̃ 𝑋 −𝑌 ( 𝑧 )= ∑ 𝜋 𝑛 𝑝 𝑋 −𝑌 ( 𝑧 |𝑀 𝑛 ) 𝑁 𝑛 −1 (1-21) The BMA distributions are plotted for the exemplar ensemble in Figure 1-13 assuming all five models are equally probable representations (  n = 1/6). Table 5 lists the corresponding lower bounds on the convergence depths for confidence levels greater than 50%. Table 1-5 Ensemble Bounds on Convergence Depths Confidence* A-B C-B PS-BC (%) (km) (km) (km) 99 - - 304 95 - - 330 90 152 188 340 80 179 218 356 70 189 233 379 60 196 246 398 50 217 268 421 * Confidence that the listed depth is less than the true convergence depth Practitioners of probabilistic seismic hazard analysis (PSHA) will recognize that (21) is analogous to averaging alternative hazard curves to get a mean hazard curve. In PSHA, these alternative probability distributions, which collectively represent the epistemic uncertainties, are usually organized as branches of a logic tree and assigned weights by expert opinion (e.g., Budnitz et al., 1997; Field et al., 2009). No finite set of models is likely to be complete, of course; instead, it is assumed that the chosen set, in our case the tomographic exemplars, ˆ z X-Y 49 adequately samples the range of real-world possibilities (Hoeting et al., 1999; Bommer and Scherbaum, 2008). The results of Table 5 rest on this assumption. Other assignments of the model probabilities {  n } are plausible. In particular, S362ANI was derived from the most informative set of observational constraints on upper mantle structure and would therefore be the model to which most expert opinion (at least ours) would assign the highest weight. For C-B and PS-BC, the S362ANI p-values are very close to the ensemble mean, so increasing its weighting would have little effect on the unconditional probabilities other than narrowing the distribution about this mean and thus increasing the lower bounds. This robustness is reassuring. On the other hand, because the A-B curves for the exemplars cluster into two groups, the ensemble probability distribution is distinctly bimodal (Figure 1-13a). The S362ANI bounds on A-B convergence depth fall in the shallower group, so that increasing its weight would reduce the ensemble lower bounds. For the PS-BC SAVANI model, vertical smoothing was done in a stepwise method, pre-setting the vertical roughness terms between 0 and 300 km at twice the weight of the horizontal damping terms, and at depths larger than 300 km linearly increasing this factor from 2 to a value of 5 (Auer et al., 2014). Since our estimates of vertical bias are based on an estimate from a radial correlation function from ~ 50 to 400 km, this may cause an overestimation of the vertical smearing bias of ~5 km. Another possible source of error introduced by this stepwise method of vertical smoothing is possible over damping of perturbations below 300 km, which may bias our convergence depth results to shallower values. Due to these concerns, an argument can be made assign a lower weight to SAVANI. 1.8 Other Sources of Bias in Convergence-Depth Estimates We have derived lower bounds on convergence depths from an ensemble of the global tomographic models by an objective procedure that corrects for vertical-smearing bias and 50 accounts for both aleatory and epistemic uncertainties in the regionalized shear-velocity profiles. In this section, we consider biases that could arise from “ontological” uncertainties; i.e., mistakes in the assumptions used to quantify the other two types of error. Robustness with respect to one of these, model weighting, has already been discussed. 1.8.1 Form of the Bias Correction The conditional bias (1-19) assumes that the true differential shear-velocity profile decreases linearly down to its convergence depth, corresponding to an exponent  = 1 in the general functional form of Appendix D. We tested this bias-correction procedure by applying it to Faul and Jackson's, (2005) theoretical models of boundary-layer structure (Figure 1-14). These authors obtained an internally consistent description of shear-velocity variations with frequency, temperature, and grain size by fitting experimental shear-modulus and attenuation data to a common mineralogical model, and they presented a set of v S profiles derived from geotherms for continents and oceans. Using their numerical code, we calculated area-weighted averages of v S (z) for the oceanic regions A, B, and C. To represent the combined P (stable platform) and S (shield) region, we averaged their models for the stable platform of northern Australia and the Yilgarn shield of southwestern Australia. These profiles and their associated geotherms are compared with the estimated convergence depths in Figure 1-16. The PS average of the Faul and Jackson, (2005) continental geotherms does not converge with the oceanic adiabat in the upper mantle (dashed blue line in Figure 1-16). To force convergence, we added a temperature correction that varied linearly from zero at the surface to 112ºC at 400 km. This choice sets the PS-BC convergence depth exactly (though arbitrarily) at 400 km (solid blue line in Figure 1-16). 51 Figure 1-14 Test of the bias-correction procedure using the differential shear velocity profiles from the Faul- Jackson theoretical models (solid lines). These theoretical profiles were smoothed using the (z) values for TX2008 (dashed lines), and the bias-correction procedure was applied to recover the convergence depths (solid points). We computed A-B, C-B, and PS-BC differential profiles from the Faul-Jackson models, smoothed them using the (z) values for TX2008, and calculated the convergence depths by applying the bias-correction scheme to the smoothed profiles. For a p-value of 0.5 (the median estimator), the bias correction based on a linear model (  = 1) yields a PS-BC depth of 382 km, which corresponds to profile differences of +0.17%. This 18-km negative bias is due to the concave curvature of the theoretical profile, which is better approximated by   1.2. If we use 52 this value in the general bias-correction model (equation D1), we recover the correct convergence depth of 400 km. Faul and Jackson (2005) assumed half-space cooling in constructing their oceanic geotherms, so their differential profiles converge only asymptotically. The intra-oceanic convergence depths are not well defined for this type of model, but we can obtain lower bounds using the one-sided test (17). For p = 0.5, a bias correction with  = 1 gives convergence depths of 157 km for A-B and 231 km for C-B, corresponding to offsets of about –0.16% and +0.20% in the differential profiles (Figure 1-14). Using these offsets to measure the sizes of the bias corrections, we find they are close to the values of the smoothing half-width , which ranges from about 50 km for A- B to about 80 km for PS-BC. We retain the linear bias model (  = 1) with the understanding that, if the true profiles are concave like the Faul-Jackson models, the estimates will be biased to slightly shallower depths by an over-correction for vertical smearing. 1.8.2 Inter-Regional Correlations Our analysis assumes that the cross-correlation in the aspherical variations between two different tectonic regions is weak, as expressed in inequality (1-11). If this is correct, then the “horizontal smearing” of the tomographic inversions can probably be neglected when considering large-scale regional averages. We can check this assumption by computing the inter- regional cross-correlation coefficient, 𝑅 𝑋𝑌 ( 𝑧 , 𝑧 ) in the notation of Appendix E. In the depth range of 200-400 km, S362ANI gives averages of 𝑅 AB  0.41 and 𝑅 CB  0.30 for the oceanic regions. These modest cross-correlations do not affect the depth estimates very much; e.g., they reduce the differential s.e.m. 𝜎 A−B by a factor of only √1 − 𝑅 AB  0.8. They allow for some lateral smearing across juxtaposed oceanic regions, but the smearing of low-velocity structure from 53 region A to region B should be roughly compensated by the smearing of high-velocity structure from C to B. The correlation between regions A and C is less than 0.1, so any bias in the intra- oceanic convergence depths should be fairly small. The intra-continental correlations are larger, as expected for the smaller continental regions. In the S362ANI upper mantle, 𝑅 PS  0.60, R QP  0.48, and 𝑅 QS  0.26. The lateral smearing allowed by the relatively high correlation between regions P and S may help to explain the similarity of their averaged v S profiles in essentially all of global tomographic models (Appendix B). Regional studies (e.g., James et al., 2001; Nettles and Dziewoński, 2008), as well as regionalizations of global data sets (e.g., Jordan, 1981), indicate that the continental shields have higher average shear velocities than the surrounding continental platforms. This difference appears to be smoothed out of the global inversions. The cross-correlation between regional combinations PS and BC is negligible (~ 0.05), but what about the correlation between regions P and Q? In the uppermost mantle, the region-Q profiles show significantly lower shear velocities than regions P and S, as expected for regions of thin lithosphere. The apparent convergence depth of region Q with the convecting oceanic mantle is not shallow, however; in fact for many of the models, it is comparable to that of stable continent (Appendix B). Because a variety of tectonic processes contribute to mantle heterogeneities beneath orogenic zones and magmatic belts, including the subduction of high- velocity oceanic plates, we will defer a full analysis of the region-Q structure. For our present purposes, it is sufficient to note that lateral smearing between region Q and the combined PS region is not likely to introduce much bias to the estimated PS-BC convergence depth. 54 1.8.3 Tomographic Inversion Filters From the angular and radial correlation functions of the individual structural models, we derive an approximation to the tomographic inversion filter implied by the model’s data set, which we then use to make the bias corrections. The main approximations are: (i) the expected value of an observed differential profile is the true differential profile smoothed only over depth (equation 13); (ii) the smoothing can be approximated by a unimodular, Gaussian kernel with a half-width ; and (iii) the depth dependence of  can be estimated from the model’s radial correlation function computed for the combined BC region. Sufficient conditions for the strictly regional relationship of equation (13), described in Appendix C, are closely related to the decorrelation condition of equation (11) and the issues of lateral smearing. However, for tomographic data set with good global coverage, the character of the inversion filter should not be a strong function of geographic position, in which case the inter-regional correlation bias must be small. According to inference principles first articulated by Backus and Gilbert, (1968), the ideal inversion filter from a good global data set would conform with (ii), which allows the essential non-uniqueness to be expressed as a localized smoothing of the actual structure. In practice, tomographic inversions often employ some form of a Gaussian-Bayesian inversion (Tarantola, 2005), rather than a Backus-Gilbert construction. The resulting inversion filters can be complex with large sidebands (under-damped), or they can be sub-unimodular (over-damped). Usually, they are both. We rely on heavy averaging to suppress the sidebands. As shown in Appendix C, the projection process averages the complete tomographic inverse filter 𝐺 ( 𝑧 , Ω; 𝑧 ′ , Ω ′ ) over both angular positions,  and ', to obtain the depth-dependent kernels 𝑔 𝑋𝑋 ( 𝑧 , 𝑧 ′ ) for each region. 55 The kernels for two regions must then be averaged to obtain the effective inversion kernel for their differential profile, 𝑔 𝑋 −𝑌 ( 𝑧 , 𝑧 ′ ) . For example, the PS-BC kernel is the bi-regional average of the inversion kernel over two-thirds of Earth’s surface. It is plausible to expect this averaging to yield kernels that are well-localized in depth, though this assertion should be verified by direct calculation. Evidence for sub-unimodularity comes from Table 3, in which the tomographic values of the integrated differentials show relative variations of 10-20%. These variations are positively correlated, suggesting that the differences among the exemplars are in part due to an “under- fitting bias.” This type of bias is common to tomographic inversions in which the maximization of the data fit is traded off against the minimization of a model norm. Owing to under-fitting bias, the smoothing kernels will not be unimodular as assumed in (14); the integral ∫𝑔 𝑋 −𝑌 ( 𝑧 , 𝑧 ′ ) 𝑑𝑧 will instead be less than unity. However, so long as this integral is a constant over the depth range of interest (here the upper mantle), it can be factored out of the bias analysis in Appendix D, and the p-value in (17) will not itself be biased. Our use of regionalized radial correlation functions to estimate (z) is novel. If our estimates are low, then the bias corrections will be too small, and the convergence depths will be biased deep. The estimation is based on an upper bound (E10) proportional to the radial correlation length (z). By assuming  equals this bound, we bias it to be high, because radial correlations associated any true Earth structure in the BC region will increase . Figure 1-11compares this bound with an independent measure of smoothing half-width derived from the S362ANI parameterization. The comparison suggests that an error of ±20% may be typical of the  estimates. Figure 1-13 shows the variations in convergence depths implied by this uncertainty, which is smaller than the dispersion of the ensemble. Y X -Y 56 1.8.4 Tectonic Regionalization Our results are based on a 5º  5º digitization of the GTR1 map. As an alternative, we constructed a 1º  1º tectonic regionalization from Müller et al.'s, (2008) magmatic ages for oceans and Artemieva's, (2006) tectonic ages for continents. To represent the A, B, and C regions, we included oceanic crust in the same age ranges as GTR1: 0-25 Ma, 25-100 Ma, and 100-205 Ma, respectively. As a proxy for the combined P and S regions (stable continent), we included all continental area with tectonic ages greater than 600 Ma. We obtained shear-velocity profiles by projecting the tomographic exemplars onto this high- resolution regionalization and applied the same procedures to estimate convergence depths. Averaged over the six exemplars, the apparent convergence depths at the 97.5% confidence interval are about 3 km shallower for A-B, 14 km shallower for C-B, and 18 km deeper for PS- BC. This comparison indicates that the convergence-depth biases from errors in GTR1 are small compared to the other estimation uncertainties. Replacing GTR1 with an alternative regionalization based on the Muller et al. (2008) and Artemieva (2006) crustal ages would tend to favor a greater separation between the intra-oceanic and continent-ocean convergence depths; in particular, it would yield a deeper value for the PS-BC convergence depth. 1.9 Discussion We have placed lower bounds on inter-regional convergence depths by detecting where the signals of tectonically correlated boundary-layer structure, corrected for vertical-smearing bias, rise above the background noise of uncorrelated lateral heterogeneity. Results have obtained by an objective analysis without regard for their implications for plate tectonics and mantle flow, which we will now consider. Table 5 lists the ensemble lower bounds 𝑧 ̃ 𝑋 −𝑌 at confidence 57 levels1-𝑝 ̃ 𝑋 −𝑌 > 50%, where the ensemble probability 𝑝̃ 𝑋 −𝑌 is the Bayesian model average (21). By definition, the median ensemble estimator of the X-Y convergence depth is 𝑧 ̃ 𝑋 −𝑌 (50%). 1.9.1 Intra-Oceanic Convergence Depths Confinement of the age-dependent component of intra-oceanic heterogeneity to the uppermost mantle cannot be rejected. The median estimates are𝑧 ̃ A−B (50%) = 236 km and 𝑧 ̃ C−B (50%) = 284 km, and the lower bounds at 90% confidence are 146 km and 178 km, respectively. The statistical analysis does not establish lower bounds on the intra-oceanic convergence depths at confidence levels above 90% because of the inconsistencies among the tomographic models (Figure 1-13) and the vertical-smearing bias, which increases as z decreases (equation 19). Models of Pacific Ocean upper mantle derived from high-resolution regional waveform inversions are consistent with global tomography (Figure 1-15). The uncertainties in the regional models have not been quantified, but the apparent convergence depths indicate a plausible range of 150-230 km, subject to some unknown vertical smearing bias. This bias is likely to be small, because the multiple-S data used in the regional inversions are especially sensitive to changes in the v S gradients at these depths (Gaherty et al., 1996; Xu and Wiens, 1997; Katzman et al., 1998; Gu, 2005; Tan and Helmberger, 2007). 58 Figure 1-15 Comparison of shear-velocity profiles from Pacific waveform inversion studies (lines) with the lower bounds on the A-B and C-B convergence depths determined from the exemplar ensemble of this study (box). East Pacific Rise (EPR) model is by Gu et al. (2005), PA5 model is by Gaherty et al. (1996), and Old Pacific and N. Fiji Basin models are by Xu and Wiens (1997). The Pacific convergence depths from these regional models (stippled) are consistent with the tomographic lower bounds, which are labeled with their confidence levels. 59 Figure 1-16 Comparison of the convergence-depth probabilities from the tomographic ensemble (left side) with the Faul and Jackson (2005) models of upper mantle temperature (middle) and shear velocity (right). The PS average of the original Faul-Jackson model (dashed blue lines) does not converge with the oceanic averages in the upper mantle, so a linear correction was added to force convergence at 400 km (solid blue lines). Colored dots on the shear-velocity profiles show the convergence depths recovered from the Faul-Jackson models in Figure 1-14. The probable depth range for the G discontinuity, 55-70 km, is from Schmerr, (2012)’s SS precursor data. The medians recovered from the Faul-Jackson half-space cooling model are 157 km and 231 km, respectively, also consistent with the tomographic bounds. The ensemble probabilities for C-B are shifted deeper than those for A-B, but this distribution difference is insignificant by any plausible statistical test, owing in part to the large discrepancies in the 𝑧 ̃ A−B values among the exemplars (Figure 1-13). More generally, the GTR1 age categories are too coarse and the tomographic convergence depths are too uncertain to discriminate between alternative models of oceanic thermal evolution, such as the plate model against the half-space cooling model (Sclater et al., 1981; Stein and Stein, 1992; Korenaga and Korenaga, 2008a). It can only be said that the A-B -B= C-B PS-BC 60 tomographic constraints on the intra-oceanic convergence depths are consistent with half-space cooling and less consistent with hypotheses implying smaller convergence depths. Stronger inferences can be derived from the hypothesis that the Gutenberg (G) discontinuity is the oceanic LAB. The G discontinuity marks a sharp drop in v S that separates the seismic lid from the shear-wave low-velocity zone (Revenaugh and Jordan, 1991) and may indicate a compositional, and therefore rheological, boundary generated by the partial melting of peridotites at spreading centers (Hirth and Kohlstedt, 1996; Gaherty et al., 1999; Karato, 2012). The G discontinuity beneath the Pacific and its marginal basins appears to be confined to a relatively narrow depth range, with most reflected-wave observations falling between 50 and 100 km (Revenaugh and Jordan, 1991; Bagley and Revenaugh, 2008; Kawakatsu et al., 2009; Rychert and Shearer, 2011; Schmerr, 2012). The lack of substantial thickening of the lid with thermal age is consistent with a chemical control of the G discontinuity. For example, Schmerr’s (2012) SS precursor data give mean G depths of 57 ± 3 km for region A, 64 ± 3 km for region B, and 64 ± 4 km for region C. These averages are much shallower than either 𝑧 ̃ A−B (90%) or 𝑧 ̃ C−B (90%) . According to the seismological models, the mantle more than 100 km below the G discontinuity participates in the conductive cooling that forms the kinematic plates. If the G discontinuity coincides with the base of the mechanical boundary layer (lithosphere), as the sharp drop in v S suggests (Hirth and Kohlstedt, 1996; Karato, 2012), then the weakest part of the asthenosphere must ride with the lithosphere in an augmented kinematical entity we will call oceanic tectosphere (Jordan, 1981a, Figure 1). The global tomographic models, as well as other structural studies (e.g., Katzman et al., 1998; Montagner, 2002; Ritzwoller et al., 2004; Cao et al., 2011), provide convincing evidence 61 for thermal anomalies in the oceanic upper mantle between the G and 410-km discontinuities. Some form of small-scale convection is likely to be active beneath old oceanic lithosphere (Huang, 2003; van Hunen, 2003; Korenaga and Jordan, 2003; Korenaga and Jordan, 2004; Landuyt and Ierley, 2012). However, the large-scale flow evidently advects the small-scale heterogeneity in the uppermost mantle along with the plates, allowing the upper part of the asthenosphere to continue cooling with lithospheric age. The dominance of this large-scale horizontal flow may be related to the high stresses associated with its channelization in a thin (~100 km) asthenosphere (Höink et al., 2012), as well as the possible focusing of the subtectospheric strain in a low-viscosity channel immediately above the 410-km discontinuity (discussed below). 1.9.2 Continent-Ocean Convergence Depths Firm lower bounds can be placed on the convergence depth between stable continents and mature oceans. We can say with 95% confidence that 𝑧 ̂ PS−BC is greater than 330 km (Table 5). The high value of this bound is notable because 𝑧 ̂ PS−BC is an average over all Precambrian crustal terrains, platforms as well as shields, and thus samples a broader distribution of continental structure than just the Archean cratons, where continental keels are thought to be the thickest. Surface-wave studies over such large regions have often been interpreted as favoring substantially shallower convergence depths, usually quoted in range of 200-250 km (Okal, 1977; Fishwick et al., 2005; Van der Lee and Frederiksen, 2005; Priestley et al., 2006). This underestimation appears to result from two factors: (i) the poor resolution of fundamental-mode data for structure below 250 km, and (ii) a high-v S bias in the reference structure used as the baseline for calculating structural anomalies. In regional studies, the latter is often set by an 62 average heavily weighted towards the continents. Because all continental regions, even Phanerozoic orogenic zones and magmatic belts (region Q), have higher average shear velocities than mature oceanic upper mantle (regions B and C) at depths greater than 250 km (Figure 1-5), cross-sections that show stable continental mantle in blue and unstable continental mantle in red will underestimate 𝑧 ̂ PS−BC . Regional models of stable continents and mature ocean basins that incorporate higher modes and multiple-S phases are generally consistent with the global tomographic constraints on 𝑧 ̂ PS−BC (e.g., Grand and Helmberger, 1984; Grand, 1994; Gaherty and Jordan, 1995). Lerner-Lam and Jordan, (1987) sought to place a firm lower bound on this convergence depth by using an inversion procedure they called “squeezing”, which forced structural differences to be small below some variable test depth z*. From an analysis of higher-mode dispersion along Eurasian and Pacific paths, they found no models with z* < 220 km that could fit the data adequately, and they concluded that convergence depth for the stable continents and mature ocean is not likely to be much less than 400 km. Squeezing experiments have not been published with the tomographic data sets and model parameterizations considered in this paper, although we hope our results will motivate such studies. A useful measure of the lithosphere’s strength is its effective elastic thickness T e (Burov and Diament, 1995; Watts, 2001). Audet and Bürgmann, (2011) mapped T e of continental lithosphere on a global scale by comparing the spectral coherence between topography and gravity anomalies and the flexural response of an equivalent elastic plate to loading. Averaging their maps over GTR1 yields T e = 36 ± 4 km for region Q, 91 ± 8 km for region P, and 78 ± 8 km for region S, in good agreement with previous summaries (Watts, 2001). Realistic modeling of lithospheric stresses indicates that these data require an LAB depth of less than 200 km beneath 63 stable continents (Burov and Diament, 1995; Watts, 2001). The stable continental lithosphere, if defined by its mechanical response—i.e., consistent with Barrell's (1914)Barrell’s original usage—is clearly much thinner than the PS-BC convergence depth obtained from mantle tomography. The regional radial correlation functions in Figure 1-10 indicate that v S heterogeneities in the upper mantle beneath stable continents are vertically coherent in the mantle above 400 km. The most straightforward interpretation identifies the PS-BC convergence depth as the thickness of the stable continental tectosphere, in which case we can conclude that the LAB is an internal boundary within the continental tectosphere, not its base (Jordan, 1975; Jordan, 1981a). Under this interpretation, the “mid-lithosphere” discontinuities imaged beneath stable continents by reflected-wave and receiver-function methods (Artemieva et al., 2006; Wittlinger and Farra, 2007; Rychert and Shearer, 2009; Fischer et al., 2010), including the Lehmann discontinuity in anisotropy (Revenaugh and Jordan, 1991; Gaherty and Jordan, 1995), could be features associated with the stable-continent LAB. The median estimate 𝑧 ̃ PS−BC (50%) = 435 km is suspiciously close to the 410-km discontinuity, which suggests that the base of the continental tectosphere is more likely to be dynamically controlled by processes associated with this phase boundary. One interesting possibility is that the decoupling of the continental tectosphere from mantle convection is regulated by a low-viscosity zone just above the 410-discontinuity, marked by a shear-wave low- velocity zone at these depths (Revenaugh and Sipkin, 1994; Vinnik and Farra, 2002; Vinnik, 2003; Song et al., 2004; Courtier and Revenaugh, 2007; Vinnik and Farra, 2007; Jasbinsek et al., 2010). The transition-zone water-pump of Bercovici and Karato, (2003) provides a plausible mechanism for lubricating the olivine-rich mantle at these depths on a global scale (Leahy and 64 Bercovici, 2007; Tauzin et al., 2010). Although little has been written about the role of this low- velocity layer in decoupling tectosphere from deep convection, its low viscosity could help to stabilize plate-driven flow in both the continental and oceanic upper mantle. 1.9.3 Correlated Mantle Flow: an Alternative Hypothesis? Gung et al., (2003) and Yuan and Romanowicz, (2010) have suggested that the cratonic tectosphere is decoupled from mantle flow by an LAB at a depth of 200-250 km, but, because the asthenospheric layer below the cratonic LAB is characterized by a strong flow-induced anisotropy, the cratons appear to be thicker than 250 km. This model retains the LAB hypothesis by explaining the high shear velocities below the uppermost mantle as a correlation between cratonic structure and mantle flow. According to this type of correlated-flow model, the layer of subcratonic flow has three properties: (a) radial anisotropy is characterized by v SH > v SV , (b) azimuthal anisotropy aligned with absolute plate motions, and (c) a configuration that makes the cratons appear to have anomalously high shear velocities at depths greater than 250 km. The evidence for any of these properties is less than compelling. Property (a) was supported by Panning and Romanowicz, (2006) using model SAW642AN; however, as noted by Kustowski et al., (2008) and Ferreira et al., (2010), the global tomographic models show substantial inconsistencies in radial anisotropy below the uppermost mantle. In particular, more recent inversions published by the Berkeley group (e.g., SEMum) found no significant cratonic signature in the radial anisotropy at these depths (Lekić and Romanowicz, 2011). Yuan and Romanowicz, (2010) used the joint inversion of SKS splitting observations with surface-wave dispersion from North America to argue for property (b). However, this interpretation depends on the assumption that variations in upper-mantle structure are rather smooth. For example, small-scale anisotropic heterogeneity in which the vertical correlation 65 length is smaller than the horizontal correlation length preferentially scatters SV waves. Such scattering decreases the magnitude of SKS splitting (see Fig. 10 of Saltzer et al., (2000)), while it increases the polarization anisotropy of the surface waves (Gaherty et al., 1996). The cratonic mantle imaged using receiver functions from dense arrays (Wittlinger and Farra, 2007; Bostock et al., 2010) shows just this type of heterogeneity. Jordan and Gaherty, (1995) and Saltzer et al., (2000) invoked small-scale, anisotropic heterogeneity to explain the apparent discrepancy between the low vertical splitting times from SKS waves and the high horizontal (polarization) splitting times from Love and Rayleigh surface waves, which is evident in many cratons and the basis for the Yuan and Romanowicz, (2010) conclusions. Moreover, property (b) is not consistent with observations of azimuthal anisotropy beneath a number of cratons (Conrad and Lithgow-Bertelloni, 2006), including the Baltic shield (Pedersen et al., 2006) and the Kaapvaal-Zimbabwe craton complex of Southern Africa (Silver et al., 2001). Rather than correlating with plate-motion directions, the azimuthal orientation of the anisotropy appears to be more aligned with ancient crustal tectonic structures, suggesting that the anisotropy beneath stable continents is frozen in the tectosphere, not caused by present-day mantle flow. Property (c) is the most problematic aspect of the Gung et al., (2003) hypothesis or, for that matter, any other correlated-flow model, because the asthenospheric flow must be arranged so that it induces anomalously high values of the Voigt-averaged shear velocities right beneath the cratons. Whereas the wave velocities in particular directions (e.g., aligned with plate motions) or with particular polarizations (e.g., v SH ) can be increased by flow-induced anisotropy, the isotropic shear velocities obtained by Voigt averaging cannot, at least to first order; they primarily reflect temperature and, to a lesser degree, composition (Jordan, 1979; Lee, 2003). 66 Therefore, it is hard to explain the high velocities observed in the depth range 250-400 km beneath the stable continents without appealing to lower temperatures. In this depth range, the global tomography models with the shallowest GTR1 heterogeneity display an average PS-BC difference in isotropic v S of about +1% (Figure 1-7). For a temperature derivative of dv S /dT = 3.3  10 –4 km/s/ºC (Lee, 2003), this corresponds to an average temperature difference of about 150ºC. In the Gung et al., (2003) model, there is no obvious mechanism for maintaining a temperature anomaly of this magnitude. 1.10 Conclusions A diverse ensemble of global tomographic models show a consistent set of correlations with large-scale tectonic regions. As described in §5, the tectonic correlations with isotropic shear velocities are very robust with respect to data selection and model parameterization, including the treatment of shear-wave anisotropy, and they are unlikely to change as new data are inverted and more detailed parameterizations are explored. This paper has demonstrated how lower bounds on the convergence depths for large-scale tectonic regions can be objectively estimated at fixed confidence levels from global tomographic models (§6). Because the convergence of two regional profiles can be asymptotic rather than abrupt, the estimation of convergence depths is an asymmetric estimation problem. We have provided a solution to this problem in the form of a one-sided convergence-depth test (17). We have shown how the covariance analysis of individual tomographic models can yield bias-corrected convergence depths and their aleatory uncertainties, and we have conducted a comprehensive ensemble analysis to estimate epistemic uncertainties in the convergence depths (§1-6). We have checked the bias-correction procedure by applying it to the Faul-Jackson theoretical models of upper-mantle v S structure (§1-7). We have also considered ontological 67 uncertainties associated with inter-regional correlations, tomographic inversion filters, and tectonic regionalizations, concluding that most would introduce a negative bias; i.e., increase the convergence depths at a fixed confidence level (§1-7). Therefore, we can assert with 90% confidence that age-correlated variations in oceanic shear velocities persist to depths greater than 170 km and with even higher confidence (95%) that the high shear velocities correlated with stable continents persist below 350 km (Table 5). The inferences regarding tectonic convergence depths have been vigorously debated for several decades. In particular, while it is now generally recognized that the cratons are underlain by thick, basalt-depleted keels in near-isopycnic balance (Forte and Perry, 2000; Carlson et al., 2005; Simmons et al., 2009; Lee et al., 2011), the maximum thickness of the keels has remained controversial since the isopycnic theory of the continental tectosphere was first proposed in the mid-1970s (Jordan, 1975; Jordan, 1978b). Several recent reviews have combined tomographic models with other information to conclude that cratonic keels achieve maximum depths of 200- 250 km; i.e., they are confined to the uppermost mantle (Eaton et al., 2009; Fischer et al., 2010; Khan et al., 2011; Lee et al., 2011). This conclusion does not agree with the results from global tomography presented here. Our meta-analysis of the published tomographic models confirms the existence of coherent continental structures below the uppermost mantle that are inconsistent the LAB hypothesis. In particular, the depth where the average shear-velocity profile beneath stable continents converges with the average profile for sub-oceanic mantle is very likely to exceed 350 km. Our results, which pertain to Voigt-averaged shear velocities, are inconsistent with the hypothesis that this continental deep structure is explained by anisotropy generated by correlated mantle flow (Yuan and Romanowicz, 2010). 68 We infer that the tectosphere beneath stable continents extends below acceptable definitions of the LAB; i.e., the kinematically coherent cratonic keels comprise both lithosphere and asthenosphere. We conclude that the plate decoupling zone beneath stable continents is closer to the top of the mid-mantle transition zone, perhaps involving a low-viscosity layer immediately above the 410-km discontinuity (Bercovici and Karato, 2003; Tauzin et al., 2010). The case against the LAB hypothesis is strengthened by lower bounds on the intra-oceanic convergence depths, which indicate that the mantle more than 100 km below the G discontinuity, which we take to be the oceanic LAB, participates in the conductive cooling that forms the kinematic plates. Oceanic tectosphere, like its continental counterpart, comprises more than lithosphere. The ensemble analysis has allowed us to identify tomographic models with anomalous structure, as well as other significant regional features (§5). The average convergence depth for unstable regions of the continents (region Q) is greater than the intra-oceanic convergence depths (Figure 1-5), although the interpretation of this observation is problematic due to correlations with mantle flow (e.g., subducted oceanic lithosphere). Another interesting feature is the peculiar correlation between continental and D" structure. The shear-velocity profiles near the CMB systematically higher beneath stable continents than beneath the oceanic regions (Figure 1-4); moreover, the intra-continental heterogeneities in D" also correlate with those in upper mantle, with correlation coefficients exceeding 0.6 for most of the exemplars (Figure 1-9). Although these intra-continental correlations may be spurious, they deserves further scrutiny in the context of hypothesized plume-mediated interactions between D" and continents (Burke et al., 2012). 69 2 Thermal Aging of the Oceanic Asthenosphere 2.1 Abstract The lithosphere-asthenosphere boundary (LAB) is traditionally defined as the boundary between Earth’s kinematic plates in the upper mantle, and a weak zone of lateral shearing, the asthenosphere, above the lower mantle. This study treats this definition as a testable hypothesis. The purpose of this study is to use global tomographic models to analyze, in a statistically rigorous way, the correlation between surface tectonic structure and the mantle beneath it. As there are many tomographic models, this study aims to evaluate an ensemble of models, analyze the similarities and differences which are observed among models, and search for globally consistent features related to Earth’s oceanic regions. This study determines the convergence depths at which perturbations in S-wave velocity (v S ) between oceanic mantle age regions could no longer be discerned. The oceanic convergence depths (>170 km) are thought to indicate the depth of the thermal boundary layer below oceanic age regions. However this is not consistent with the standard plate cooling model, some seismic studies, and it is deeper than the observed G-discontinuity. Small-scale convection may still create heterogeneities in the oceanic upper mantle; however, the large-scale flow evidently advects these small-scale heterogeneities along with the plates, allowing the upper part of the asthenosphere to continue cooling with lithospheric age. The dominance of this large-scale horizontal flow may be related to the high stresses associated with its channelization in a thin (~100 km) asthenosphere, as well as the possible focusing of the subtectospheric (sub-thermal boundary layer) strain in a low-viscosity channel immediately above the 410-km discontinuity. The coupling of the lithosphere and asthenosphere for times scales on the order of the age of the oceans is also supported by seismic observations of apparent subduction of 100+ 50 km thick asthenospheric mantle at convergent 70 margins. The observed thermal aging of oceanic asthenosphere is inconsistent with the LAB hypothesis, which states that lithospheric plates are decoupled from deeper mantle flow by a shear zone in the upper part of the asthenosphere. Our results show (at 90% confidence) that age-correlated variations in oceanic shear velocities persist to depths greater than 170 km. We also conclude that the thermal structure of the oceanic asthenosphere above the oceanic convergence depths decays on average linearly with the square-root of age to >200 my, indicating that the upper part of the oceanic asthenosphere participates in the cooling that forms the kinematic plates (tectosphere). The ocean basins show convergence depths much greater than the oceanic LAB. 2.2 Introduction Part 1 of this dissertation concerned the analysis of the correlation between the large-scale tectonic structure of Earth, as delineated by GTR1, and Earth’s aspherical seismic S-velocity structure. It was found that the Earth’s aspherical variations of shear velocity in the uppermost mantle correlate strongly with GTR1. In particular, for Earth’s oceanic regions, the aspherical shear velocity (v S ) variations were most negative beneath young oceans (region A) and increased systematically with crustal age to least negative beneath old oceans (region C). The models exhibited an “apparent oceanic convergence depth”, below which depth profiles of the three oceanic regions (A, B, C) showed no statistically significant variations. These convergence depths, when corrected for vertical smearing bias, persisted to depths >170 km. This finding is only partially explained by seafloor topography, and is inconsistent with some models of the cooling structure beneath Earth’s ocean basins as well as the traditional definition of the LAB hypothesis. 71 Seafloor topography, which reflects the thermal state of the oceanic lithosphere, becomes shallower than predicted by simple half-space cooling after ~70 million years (Parsons and Sclater, 1977; Stein and Stein, 1992). This trend is consistent even with the removal of areas affected by hotspots and oceanic plateaus, although some seafloor appears unaffected by this flattening (Korenaga, 2015). There is also an observed heat flow decrease with age of the oceanic floor (Sclater and Francheteau, 1970). Our observation of deep (>170 km) convergence depths for oceanic age regions is inconsistent with the observed flattening of seafloor topography, which indicates the presence of half-space cooling to depths much greater than that the lithosphere-asthenosphere boundary (LAB). Several models of have been proposed for the cooling structure beneath the ocean basins. One of the models which is consistent with our results is the half-space cooling model (Turcotte and Oxburgh, 1967; Parker and Oldenburg, 1973) for oceanic thermal structure, which delimits plate thickness by an isotherm, which allows the plate thickness to grow indefinitely with age. Other models which are inconsistent with our results are the standard plate cooling models (Houseman et al., 1981, Parsons and McKenzie, 1978, Stein and Stein, 1992) , which sets a temperature at a fixed depth as the bottom of the plate, allowing for conductive cooling of the lithosphere above this fixed depth. The plate cooling models attempt to explain observed oceanic seafloor topography, which shows subsidence consistent with a √ 𝑡 dependence up to approximately 70 million years, with a flattening of seafloor topography past this age. This flattening is explained as a result of small-scale convection at the bottom of the lithospheric plate destabilizing the half-space cooling profile of the plate bottom, allowing flattening of the overlying seafloor bathymetry past these ages (Korenaga, 2015). The fixed temperature depth is set between 95km (Stein and Stein, 1992) and 125 km (Parsons and McKenzie), depending on 72 the details of the thermal model used. Our results which show significant differences between ocean age regions above 170 km are not consistent with the shallower plate depths proposed by these models. Another model of the cooling structure of the oceanic basins is the CHABLIS model of Doin and Fleitout, 1996, which models subplate convection as a constant heat flow boundary condition. This model allows for continued subsidence beyond 70 million years, but results in plate thicknesses on the order of 80 km, (Doin and Fleitout, 1996). While the plate model and the CHABLIS model can explain the flattening of seafloor bathymetry observed beyond 70-100 million years, they are inconsistent with our convergence depth observations, which indicate thermal cooling is taking place in the asthenosphere well below the depths of the thickest plates allowed by either the plate cooling models or the CHABLIS model. In addition to inconsistency with the oceanic models for thermal structure, our results also are inconsistent with the traditional LAB hypothesis, which states that the Earth’s lithospheric plates are decoupled from deeper mantle flow by a shear zone in the upper asthenosphere. These depths are well below the 40-75 km depth of the seismic G discontinuity (Revenaugh and Jordan, 1991; Beghein et al., 2009; Schmerr, 2012) , or various recent estimates of LAB depth, whether derived from receiver functions (55-80 km) (Kawakatsu et al., 2009; Olugboji et al., 2016), P-S converted phases (70-81 km) (Rychert and Shearer, 2009), magnetotelluric data (45- 70 km) (Naif et al., 2013), seismic ambient noise and teleseismic wave forms (40-80 km) (Takeo et al., 2016) or surface wave dispersion data (20-130 km) (Burgos et al., 2014). This study attempts to describe the correlation of Earth’s surface tectonics and seismic structure for oceanic regions in finer detail than in Part 1. First, for a suite of models, we utilize a finer age regionalization than that used by GTR1 in order to investigate the oceanic v S convergence depths for those regions. Then, in order to minimize the effects of vertical 73 smearing present in the tomographic inversion process, we analyze v S travel times between 50- 300 km for a suite of models by both age and ocean basin. We then compare our results for a suite of tomographic models to those of a √𝑡 dependent model, which is consistent with the convergence depths estimated in Part 1 of this study. We also evaluate these results to see if we observe the flattening signal (Korenaga, 2015) seen in seafloor topography and whether we can see evidence of a Pacific reheating event as described by Ritzwoller et al., (2004). 2.3 Model Selection Table 2-1 Tomographic models selected for oceanic correlation analysis Inner Scale Apparent Convergence Depths Group Modelname Horizontal Vertical S- type λ OAR –All/ Pacific/ Indian/ Atlantic rcf BC (50- 300 km) A-B C-B References Harvard S362ANI 11° 60-400 km Voigt 11.5/ 13.6/ 9.0/ 9.0 34 192 252 (Kustowski et al., 2008) Berkeley SEMum2 8° 100-300 km Voigt 7.3/ 7.8/ 6.4/ 6.0 23 244 225 (French et al., 2013) Ritsema S40RTS 5° 60-300 km isotropic 10.3/ 11.7/ 7.9/ 7.2 51 220 233 (Ritsema et al., 2010) Auer SAVANI 1.25° -5° variable layers Voigt 10.0/ 12.0/ 8.2/ 7.8 31 189 215 (Auer et al., 2014) Debayle DR2012 400 km 50-500 km Sv 9.5/ 10.0/ 8.3/ 6.3 16 172 198 (Debayle and Ricard, 2012) Schaeffer SL2013sv 250km- 290 km 7-198 km Sv 4.9/ 4.9/ 4.6/ 4.3 38 175 188 (Schaeffer and Lebedev, 2013) The models in Table 2-1 differ in their representation of upper mantle anisotropy. The Schaeffer model SL2013sv parameterized the elastic structure solely in terms of SV velocities 74 and relied on measurements of surface and S, and multiple S (up to S 7 ) wave forms from ~750,000 broad-band seismograms from 11-450 s. The Debayle model parameterized the elastic structure solely in terms of SV velocities and relied solely on measurements of approximately 375,000 Rayleigh waveforms at fundamental and higher modes. The Ritsema model S40RTS inverted data sensitive to both SH and SV velocities, including ~20,000,000 Rayleigh wave phase velocity measurements (sensitive to the upper part of the mantle) as well as 500,000 S velocity travel time measurements (primarily sensitive to the lower mantle) and assumed isotropic perturbations to a radially anisotropic structure. Another set of authors inverted for lateral heterogeneity in the radial anisotropy: the Berkeley model SEMum2 which used full- waveform inversion, with synthetic seismogram computation using the Spectral Element Method, relies on ~100,000 long period seismic waveforms (fundamental mode Love and Rayleigh, long period overtones, and long period (>60 s) body waves) and ~100,000 group velocity dispersion maps between 25s and 150 s, and images global radially anisotropic shear velocity (v S ) structure at upper mantle and transition zone depths; the Harvard model S362ANI uses ray theory which relies on ~250,000 surface wave phase velocities (35-150 s) to constrain structure above ~300 km, ~50,000 long-period waveforms (including body waves (>50 s), and mantle wave (>125 s), including ~1,000 mantle waves (>200 s)) to constrain structure of the transition zone and the uppermost lower mantle, and ~140,000 body wave travel times (~20 s) to constrain lower mantle and the topography of the transition zone discontinuities, they inverted for isotropic and anisotropic variations in shear wave velocity, constraining anisotropy to vanish at 410 km; and Auer’s SAVANI model which relies on ~10,000,000 surface wave phase delays from fundamental modes up to the sixth overtone between 25 and 350s, as well as ~500,000 cross-correlation travel times of major body wave phases (S, SS, SKS, ScS, sScS, ScS2, ScS3, 75 ScS4, (SSm), (SSSm), (SSSSm)), it inverted for radially anisotropic velocity perturbations and made ray-theoretical corrections for azimuthal anisotropy before inversion. We limited the perturbations for the radially anisotropic 3-D models to their approximate Voigt averages, 𝛿 v S = 𝛿 v SH +2𝛿 v SV 3 . For all models we subtracted out (if necessary) the spherical average of the isotropic shear velocity perturbation at each depth z to ensure that the variation was aspherical: ∫ 𝛿 v S ( 𝑧 , Ω) 𝑑 Ω = 0 𝑠 1 (2-1) Here Ω = ( 𝜃 , 𝜑 ) is the geographic position on the unit sphere S 1, θ is the colatitude, and dΩ=sin θdθdφ is an element of solid angle. Figure 2-1 plots the aspherical variations of δv s at z=150 km for the ensemble of models. The maps show a range of amplitude and lateral roughness, reflecting differences in the data sets, model parameterizations, and dampings used in the inversions. 76 Figure 2-1 Selected tomographic models at 150 km depth 2.4 Projection Method In Part 1, the tomographic models were projected onto a 5º  5º representation of the GTR1 global tectonic map (Jordan, 1981b), which divided Earth’s surface into three oceanic and three continental domains (Figure 1-2). In order to estimate large-scale convergence depths, oceanic crust was partitioned according to the square-root of its magmatic age, as inferred from marine magnetic anomalies: region A (0-5 Ma 1/2 ), region B (5-10 Ma 1/2 ), and region C (10- 15 Ma 1/2 ). For this study, we created a regionalization with higher temporal and spatial resolution, Ocean Age Regionalization 1 (OAR1). OAR1 was derived from a 1x1° map of oceanic crustal ages from (Müller et al., 2008)’s paper divided into bins delineated by equal 77 increments of the square root of crustal age. The bins were designated as follows : A1, youngest oceans, (0-6.25 my or 0-2.5 my 1/2 ), A2, young oceans, (6.25-25 my or 2.5-5 my 1/2 ), B1, younger intermediate oceans (25-56.25 my or 5-7.5 my 1/2 ), B2, intermediate oceans (56.25-100 my or 7.5-10 my 1/2 ), C1, younger old oceans (100-156.25 my or 10-12.5 my 1/2 ), and C2, old oceans (>156.25 my or >12.5 my 1/2 ). We removed marginal basins/seas from the model in order to avoid the effects of subducting slabs and back-arc basins on our average S-velocities, as we are interested purely in the thermal cooling of oceanic crust with age as it moves away from the ridge. In order to estimate the effects of the presence of marginal basins on our results, we created a mask of back arc basins and inland seas and calculated travel times for the oceans including marginal basins, without marginal basins, and for the marginal basins alone. The removal of the marginal basins from the travel time calculations did not significantly affect the average S-velocities of the OAR1 regions. The six regional averages were calculated using 𝑓 𝑋 = 1 𝐴 𝑋 ∫ 𝐹 ( Ω) 𝑑 Ω 𝑋 , X  {A1, A2, B1, B2, C1, C2}, (2-2) where F( ) is the radial functional and the realization of a Gaussian random field with random process fluctuations, the regional average f X is a sample from a Gaussian random variable and 𝐴 𝑋 = ∫ 𝑑 Ω 𝑋 is the area of the region X. 78 Figure 2-2- Regionalization in 1˚x1˚ and 2.5√My age bins based on oceanic ages from (Müller et al., 2008) with marginal basins excluded 2.5 Profile Variance Analysis In order to perform the profile variance analysis for our selection of models, we followed the procedure detailed in Part 1, which calculates the correlation angle  X (z) and s.e.m.  X (z) from discretized values for the model ensemble. We then created differential velocity profiles for OAR1 regions (Figure 2-4) for each model. We also calculated the radial correlation length (z), the half-width of a model’s radial correlation function (Puster and Jordan, 1997). In this paper, we use the radial correlation lengths for individual tectonic regions,  X (z), to constrain the inversion bias in the regional profiles. r 79 Figure 2-3- Shear-velocity profiles and 1-sigma error bars obtained by projecting the models onto OAR1. The shear velocities increase as the square root of crustal age with apparent convergence depths greater than 200 km. 2.6 Tectonic Correlations Identical correlation analyses were performed on the 6 tomographic models in Table 2-1. We display the regional profiles above 700 km of the relative aspherical variations v S /v S for all of the models using the OAR1 in Figure 2-3. The profiles for each model are banded with their 1 errors, computed by the method described in Part 1. For OAR1, we find that the non-bias corrected shear velocities increase as the square-root of crustal age, and the age-correlated variations persist to >170 km (Figure 2-3), more than 100 km below the estimated depth of the G-discontinuity, (~70 km (Revenaugh and Jordan, 1991; Schmerr, 2012), and nearly 50 km below the plate bottom for the thickest plate cooling model (~125 km, (Parsons and Sclater, 80 1977) ). These are similar results to those obtained in Part 1, albeit with a different set of tomographic models. Regional features strongly correlate with oceanic age regions in the OAR1 projections. The youngest oceans (region A1) have the slowest aspherical shear-velocity structure, while the oldest oceans (region C2) have the fastest shear-velocity structure. For OAR1, most of the models show some overlap of the 1-σ uncertainty bands between the two oldest age regions (C1-C2) and the two youngest age regions (A1-A2) above the estimated convergence depths, This is likely due to correlation between these age regions caused by the relatively small surface area of these youngest and oldest regions and their proximity to each other, for example, in the SAVANI model regions A1 and A2 have correlation coefficients of 0.51 and regions C1 and C2 have correlation coefficients of 0.32, these correlations are similar to those obtained for the rest of the models in this study. Some of the models (SEMum2, SAVANI, DR2012, SL2013sv) show significant regional fluctuations below the estimated convergence depths (Figure 2-3 and Figure 2-4), these may be due to lack of resolution below 300 km for the surface wave dominated models (DR2012 and SL2013sv, and SAVANI), the change of vertical resolution imposed below 300 km in the SAVANI model (Auer et al., 2014) or artefacts of the model inversions. 2.7 Estimation of convergence depths All of the models in Figure 2-3 show oceanic regional profiles which converge above 400 km. We followed the method of Part 1, in calculating the apparent convergence depths for the models with the OAR1 projection. Due to strong correlations between adjacent regions, we combined regions A1 & A2, B1 & B2, and C1 & C2 when calculating these convergence depths which lowered the inter-regional correlations (Table 2-4). We estimated two convergence depths: 𝑧 ̂ 𝐴 1𝐴 2−𝐵 1𝐵 2 , between young and intermediate-age ocean basins and 𝑧 ̂ 𝐶 1𝐶 2−𝐵 1𝐵 2 , between 81 old and intermediate-age ocean basins. The A-B and C-B differential profiles for the six models are shown in Figure 2-4 and the convergence depths are listed in Table 2-1. Figure 2-4- Differential shear-velocity profiles and 1-sigma error bars obtained by projecting the models onto OAR1’s combined A, B, and C regions. The shear velocities increase as the square root of crustal age with apparent convergence depths greater than 200 km. To correct for vertical smearing bias caused by the tomographic models’ inversion filter, we follow the method of Jordan and Paulson, (2013) described in Part 1, where the inversion process is approximated as a linear low-pass filter of the differential profile : 82 𝜙 𝑋 −𝑌 ( 𝑧 )≈ ∫ 𝑔 𝑋 −𝑌 ( 𝑧 , 𝑧 ′ ) 𝑓 ̂ 𝑋 −𝑌 ( 𝑧 ′ ) 𝑑𝑧 ′ (2-3) where the expected value of the observed differential regional profile, 𝛷 𝑋 −𝑌 , is the true differential profile 𝑓 ̂ 𝑋 −𝑌 filtered by a differential smoothing kernel 𝑔 𝑋 −𝑌 . To calculate bias- corrected p-values, we used 𝛹 𝑋 −𝑌 = ∫ 𝑓 𝑋 −𝑌 ( 𝑧 ) 𝑑𝑧 𝑧 1 𝑧 0 (2-4) which yields the conditional bias formula, (2-5) where z 0 is the minimum depth of integration and is the upper half-width of the kernel . We adopted z 0 = 50 km as the minimum depth of integration and chose the upper limit z 1 to be the apparent convergence depth at a relatively low confidence level (50-60%). The smoothing kernel g X-Y was approximated using the bounds on the smoothing half- (z) calculated from the regional radial correlation properties between 50 and 400 km depth of the models where (z) is approximated by where ρ x 2 (z) is the half-width of the regional radial correlation function (Figure 2-5). The BC-regional radial correlation functions for the models were fit to a linear model of the smoothing half-widths, (z) =  0 + c (z – z 0 ). using the bias parameters in Table 2-2. b X -Y (z) = 2 π Y X -Y g(z) (z-z 0 ) 2 ˆ g g X-Y ( ˆ z X-Y ,z') r X (z) / 2 83 Figure 2-5 Regional radial correlation functions for mature ocean basins (region BC), computed from relative aspherical variations vS/vS for the model ensemble. 0.78 contours are outlined in black. Table 2-2 –Model bias parameters- calculated for depths between 50-300 km Model Ψ ( k m )  0 *  c A-B C-B  S40RTS -2.1 2.1 20 0.19 S362ANI -3.6 2.7 29 .08 SEMum2 -3 2.7 11 0.14 SAVANI -3.0 2.8 15 0.11 DR2012 -2.5 2.3 21 0.00 SL2013sv -3 2.3 12 0.18 84 Applying the bias correction as described in Part 1 to the oceanic OAR1 regions yields the p X-Y (z) curves of Figure 2-6. The profiles show a consistent deepening of convergence depths with age. In order to estimate the ensemble lower bounds on convergence depths, we used an estimate of the Bayesian model average of all of the models. 𝑝 ̃ 𝑋 −𝑌 ( 𝑧 )= ∑ 𝜋 𝑛 𝑝 𝑋 −𝑌 ( 𝑧 |𝑀 𝑛 ) 𝑁 𝑛 −1 (2-6) assuming all six models are equally probable (  n = 1/6). These are plotted on Figure 2-6 as black lines. These ensemble convergence depths also show a consistent increase in depth with the age of the region. 85 Figure 2-6- Convergence depths of the differential profiles for A-B (left), C-B (right) as a function of p-value, computed from using the method described in Part 1 and bias parameters of Table 2-2. The black line is the ensemble mean (Bayesian model average), which combines the aleatory and epistemic uncertainties. We interpret the p-value as the probability that a test depth z is greater then 𝑧 ̂ 𝑋 −𝑌 . Most of the models in Figure 2-6 show a consistently shaped probability curve, with 𝑧̂ 𝑋 −𝑌 increasing with increasing probability. However, SEMum2’s A-B’s probability curve truncates at ~0.6. The reason for this is apparent in Figure 2-4 where the A-B never equals zero and runs nearly parallel with the zero line below ~230 km depth. For this model, there is no convergence depth at which we can guarantee irreducible probability 2.8 Travel Time Variance Analysis In order to study the relationship between oceanic age regions, while negating the effects of vertical smearing bias and minimizing the effects of correlation between tectonic regions, we 86 chose to analyze vertical travel times 𝑇 = ∫1./𝑣𝑠 ( 𝑧 ) 𝑑𝑧 between 50 and 300 km depth. The use of vertical travel times to integrate the velocity structure between these depths minimizes the effects of vertical smearing bias as this depth range allows us to include all structure from slightly below the estimated convergence depth to 50 km depth as well as reducing the correlation between age ranges (an example for SAVANI is shown in Table 2-3). This range is thought to encompass the oceanic asthenosphere, while excluding most of the oceanic lithospheric structure. A 1˚x1˚ map of vertical travel times was created from each velocity model in this study. We then calculated average travel times and standard deviations for OAR1 regions using the method outlined in the Part 1 of this paper. Figure 2-7 shows the aspherical travel times for the models as well as the Bayesian ensemble average, by OAR1 region, along with the 1σ error bar for the ensemble average. The travel times show a clear increase in travel times with a decrease in age, consistent with a cooling asthenosphere. The black error bars show the aleatory uncertainty in the Bayesian model average (average of the individual aleatory uncertainties for the models) while the scatter of the models estimates the epistemic uncertainty (variability between models). The aleatory 1σ Ensemble error bars encompass most of the epistemic variability of the models. 87 Table 2-3 Correlation between aspherical variations in TT’s for OAR regions from 50-300 km depth for SAVANI A1 A2 B1 B2 C1 C2 Q P S A1 1.00 0.51 0.19 0.09 0.02 0.01 0.06 0.04 0.04 A2 1.00 0.51 0.18 0.05 0.02 0.13 0.04 0.02 B1 1.00 0.49 0.11 0.02 0.19 0.05 0.06 B2 1.00 0.36 0.03 0.17 0.06 0.09 C1 1.00 0.32 0.18 0.07 0.11 C2 1.00 0.11 0.06 0.03 Q 1.00 0.42 0.21 P 1.00 0.62 S 1.00 Table 2-4 Correlation between aspherical variations in TT’s for combined OAR regions and GTR1 continental regions from 50-300 km depth for SAVANI A B C Q P S A 1.00 0.38 0.05 0.10 0.05 0.04 B 1.00 0.23 0.25 0.07 0.09 C 1.00 0.21 0.09 0.11 Q 1.00 0.42 0.21 P 1.00 0.62 S 1.00 88 Figure 2-7 - OAR1 region ensemble and individual aspherical travel times with 1σ error bars 89 2.9 Tectonic TT Correlations Figure 2-8 - Regionalization of individual ocean basins In order to study differences between ocean basins, the travel time models were further divided into the three major ocean basins (Pacific, Indian, and Atlantic), as seen in Figure 2-8 with the regional averages calculated as in equation (2-2) . Figure 2-9 shows the aspherical travel times for the models as well as the Bayesian Ensemble average, by ocean basin, along with the 1σ error bar for the Ensemble average. The mean travel times for each basin are significantly different from each other, with the Atlantic basin showing the shortest average travel time and the Pacific showing the longest. Assuming the velocity temperature relationship from Lee, 90 (2003) of 3.3 x 10 -4 km/s/˚C, we get estimated average temperature differences of 36 ˚C between the Indian and Atlantic oceans, 91˚C between the Pacific and Indian Oceans, and 128˚ between the Pacific and Atlantic oceans, for depths between 50-250 km. These differing travel time profiles by basin, are consistent with the range of upper mantle potential temperatures estimated for varying oceanic ridge areas of 150˚C - 200˚C (Dalton et al., 2009) . Our numbers are slightly lower than the range of estimates from Dalton et al., (2009). This is likely due to our averaging of velocities from 50-300 for all regions, while (Dalton et al., 2009) limited their depth of interest between 150 and 200 km. The travel times by ocean basin ensemble’s aleatory variability encompasses the majority of the epistemic variability observed in the scatter shown by the model ensemble. 91 Figure 2-9- Ensemble and individual travel times by ocean basin with 1 σ error bar 2.10 Travel Time by Age Profiles In order to further analyze oceanic structure we plotted the tomographic model and ensemble travel times by both OAR age delineation and ocean basin in Figure 2-10 below. A consistent decrease in travel time with depth is observed, consistent with a cooling asthenosphere, as well as distinct differences in the vertical travel times between ocean basins. The travel time curves in Figure 2-10 show a flattening toward shorter travel times in the lowest age ranges associated with the truncation of the upper limit of integration at 50 km. A similar 92 effect is seen at the oldest ages, where cold lithosphere likely deeper than 50 km is included in the travel time calculations, decreasing the travel times at these ages, causing a downward kink in the travel times. The exception to this at the oldest ages is the Indian Ocean, which shows an upward kink, this is likely due to the extremely small surface area of Indian Ocean in this age range, allows the angular smearing from the previous age bin to dominate the v S values for this age bin (see Figure 2-12). Figure 2-10- Travel time curves by OAR age bins and ocean basin for ensemble and individual tomographic models. 93 In order to look at finer age increments of the travel time profiles, we plotted the tomographic model and Ensemble travel times in by √5 million year age bins, as well as ocean basins in Figure 2-11. Similar cooling profiles and flattening at youngest and oldest ages are observed as seen in the OAR age delineations. However, smaller-scale features are also observed, which will be discussed later in this paper. In order to determine if our travel time profiles follow a half-space cooling model, while also accounting for the effects the tomographic smoothing process for the basins, and basinal averaging for the PacIndAtl model due to changes in relative weighting between individual basins (see Figure 2-12), we chose to compare the tomographic results to an smoothed half-space cooling model as described below. 94 Figure 2-11- Travel time curves by 5 million year age bins and ocean basin for ensemble and individual tomographic models 95 Figure 2-12- Relative surface area contributions for OAR and √𝟓 my age bins by ocean basin. Bins are centered at mean area-weighted age location within each bin. 2.11 Best-fitting HSC Models In order to compare an ideal half space cooling travel time model to our tomographic travel time models, we created a “raw’ vertical S-wave travel time model for each tomographic model, assuming a strict square-root of age relationship for S-wave travel times for the ocean basins. We calibrated the slope of the initial square-root of age half-space cooling (HSC) model to that of Faul and Jackson, 2005, in order to limit the slope of the initial model to a physically reasonable, experimentally determined linear model. We then replaced the non-oceanic regions 96 with travel times equivalent to those for the individual tomographic models. A smoothed model was then obtained by filtering the raw model using a half-width derived from the model’s angular correlation function as an estimate of the model smoothing as shown in Figure 2-13 below. Figure 2-13 Method used to create the “ideal” travel time model. A strict square root of age ocean velocity model is created, then the continental velocities are derived from the tomographic model, and the resulting raw model is smoothed to create a smoothed ideal travel time model. We compared our tomographic model travel time results to this idealized travel time model. Figure 2-14 illustrates the resulting travel time curves for each ocean basin for the mean of the ensemble models. The travel times from the filtered model fit the tomographic models well and show that the S travel time above the estimated convergence depths decays nearly 97 linearly with the square-root of crustal age out to 225 million years. The individual ocean basins vary in quality of fit, but the best-fitting smoothed half-space cooling model fits within the 1-σ error bars for the mean tomographic model’s aleatory variability. Figure 2-14- Best fitting ensemble smoothed tomographic models (red solid line) compared to initial tomographic model (blue solid line) by ocean basin. The HSC input model pre-smoothing is shown as a red dashed line. The √𝟓 my tomographic age bins are shown as blue dots. 98 Figure 2-15 95% error ellipses on the slope and intercept of the best-fitting half-space cooling model for the individual models in the ensemble. Figure 2-15 illustrates the 95% error ellipses for the best-fitting slopes and intercepts used for the smoothed half-space cooling models. The slopes were pre-set to 0.20 in order to fit the slope of a travel time curve derived from Faul and Jackson, 2005’s velocity model. It is clear from the highly elliptical shape of the error ellipses that a wide range of initial slopes result in a reasonable fit to the tomographic models, necessitating the calibration of slopes to a physically reasonable, non-seismically determined value. One commonality between models is the consistent change in intercept value of the initial HSC model with ocean basin. In every case, the Atlantic has the lowest intercept, followed by the Indian, and then the Pacific Ocean with the highest intercept. These intercept differences result in travel time profiles for the oceans which are shortest for the Atlantic Ocean, and longest 99 for the Pacific Ocean. Assuming the primary driver of the travel time differences is adiabatic temperature rather than mineralogical differences (Dalton et al., 2014), this would result in a colder Atlantic Ocean, intermediate Indian Ocean, and hotter Pacific Ocean. The error ellipses overlap between the Atlantic and Indian oceans for most of the models, indicating their differences are not significant at the 95% confidence level, however, the Pacific Ocean is significantly different from the Indian and Atlantic Oceans for most models. This difference necessitates a change in our calculation of the combined PacIndAtl basin’s best fitting HSC values. Previously, we had treated the combined basins as one large single basin, fitting a single straight-line half space cooling model to it. However, given the baseline differences in intercepts for each ocean basin, as well as the large surface area differences between ocean basins, it is more accurate to treat each ocean basin as a separate entity within an aggregate Pacific-Indian- Atlantic basin. Figure 2-16 illustrates this improvement for the SAVANI model. Going forward, we will be referring to these improved fits for the combined Pacific-Indian-Atlantic region. 100 Figure 2-16 Initial fit of HSC travel time curve, assuming a purely linear starting HSC model, compared to the travel time curve for the SAVANI model’s combined PacIndAtl region, assuming a half-space cooling starting model which varies for each ocean basin. This results in a slightly non-linear starting HSC model, and better fits of our smoothed model to the tomographic model. The results for the remainder of the models are plotted in Appendix G. 2.12 Discussion For the large scale OAR1 age bins, we achieve a reasonable fit for a half-space cooling model to our tomographic model values. However, when we look at the finer scale √5 million year age bins there appear to be several small-scale departures from the larger scale averages (Figure 2-11Figure 2-14). In order to ascertain whether these “bumps” in the tomographic models correspond to real features in the Earth’s upper mantle, we discretized our best fitting smoothed half-space cooling model into √5 million year age bins. Figure 2-17 shows the mean 101 ensemble tomographic model aspherical travel times plotted in √5 million year age bins, compared to the smoothed best fitting half-space cooling model. It can be seen that the “bumps” at ages greater than ~140 my are recreated by the smoothing process in the half-space cooling model, although somewhat dampened in intensity. For the individual ocean basins (Pacific, Indian, Atlantic), this anomaly is thought to be caused by smearing of continental velocities into oceanic regions by the tomographic inversion process as well as our starting depth of 50 km, which at older ages (> ~120 my), causes our travel time models to integrate some lithospheric structure. The fact that we can recreate these bumps in our half-space cooling model shows us that they are not caused by actual asthenospheric structure, but by the smoothing process itself. Another part of the age profile that shows a slight bump is at approximately 64 my in the Indian and Atlantic ocean profiles, these are also recreated in the smoothed half-space cooling model which makes them likely artefacts of the smoothing process as well. The Pacific profiles are somewhat more problematic. There is a small (~0.15 s), but noticeable increase in travel times in the ensemble mean tomographic model with respect to the half-space cooling profile between the ages of 70-110 my. Researchers (Nagihara et al., 1996; Ritzwoller et al., 2004) have identified a possible reheating event in the Pacific Ocean between ~70 Ma and ~100 Ma. We do see an increase in travel times over this time period in the Pacific, consistent with a possible reheating event, but its magnitude is smaller than the estimated 1-σ error bars for this age range in the Pacific. This does not discount the possibility of a reheating event’s occurrence, but due to averaging over the large surface area of the Pacific ocean, as well as the integration of vertical seismic structure between 50 and 200 km depth in our travel time calculations, the finer details of such an event may be averaged out the global tomographic models. 102 Figure 2-17 Mean ensemble tomographic models (blue dots), plotted against mean best fitting ensemble half-space cooling models (red line). Many of the bumps in the tomographic models are recreated in the half-space cooling models post-smoothing, indicating they are the result of the smoothing process rather than actual physical features. In Figure 2-18 we look for the presence of this reheating event using the SAVANI tomographic model sorted by ocean basin between 50 and 300 km depth, sorted into √5 my age bins. We compare these figures to a smoothed half-space cooling model. There does not appear to be a strong reheating event in the Pacific at roughly 70 my, nor are there strong anomalies present in the other ocean basins. This does not preclude the occurrence of such an event, however, it does imply that if such an event has occurred, it is either weak or is not basin-wide in the Pacific Ocean. 103 Figure 2-18 Aspherical V S for the SAVANI tomographic model by ocean basin (column1); half-space cooling oceanic models smoothed using κ and vertical smoothing parameters calculated from SAVANI (column 2); Korenaga, 2015’s study on the thermal budget of the Earth included a map delineating “normal” seafloor from seafloor affected by hotspots or oceanic plateaus. In order to determine if this “abnormal” seafloor has a significant effect on our results we calculated average travel times for all seafloor vs. normal seafloor alone. Our results show very little difference in travel time profiles as a result of removing anomalous seafloor, indicating the effects of hotspots and oceanic plateaus is negligible to our travel time results. It also appears to indicate that the if the 104 reheating event at ~70 my in the Pacific is a robust feature, it is not related to the hotspots and oceanic plateaus, as there is no change in the travel times between the two regionalization. Figure 2-19 All seafloor vs. normal Seafloor (as defined by Korenaga and Korenaga, 2008b). The removal of the anomalous seafloor has a very small effect on the travel time profiles except at the oldest seafloor ages in the Indian Ocean, where the surface areas of the age regions are very small (Figure 2-12) so small fluctuations in travel times may have a disproportionally large effect on the average for these small area age bins. Our results are consistent with the half-space cooling model as the best-fitting thermal model for the oceanic asthenosphere. This is problematic since seafloor topography, greater than 70 million years old, appears shallower (Parsons and Sclater, 1977; Stein and Stein, 1992; Crosby et al., 2006) than that predicted by a simple half-space cooling model. If half-space 105 cooling is taking place to depths of at least 170 km, as indicated in this study, we need an explanation for this topographic flattening while maintaining half-space cooling as the thermal model. This flattening effect has been postulated to be a result of small scale convection. Korenaga, 2015 discusses the “flattening effect” visible on seafloor greater than 70 million years old, and discusses how this deviates from a pure half-space cooling model. He finds that the flattening signal can be explained by half-space cooling modified by small-scale convection fueled by basally heated convection in the presence of radioactive isotopes, which allows the cooling of the lithosphere to deviate from pure half-space cooling. We calculated travel time by age curves for both normal and anomalous seafloor as defined by Korenaga. In our vertical travel time models, we do not see a flattening signal in the integrated travel time models beyond 70 my. The lack of a flattening signal may indicate that this small-scale convection which contributes to the flattening of the seafloor is advected along with the cooling asthenosphere. This coupling of the lithosphere and asthenosphere for times scales on the order of the age of the oceans is also supported by seismic observations of apparent subduction of 100+ 50 km thick asthenospheric mantle at convergent margins (Song and Kawakatsu, 2012; Song and Kawakatsu, 2013). 2.13 Conclusions This study analyzed the correlation between oceanic age-based regionalizations and Earth’s seismic S-velocity and S travel times (from 50-300 km) for a suite of tomographic models. We conclude that the convergence depths of average S-velocity profiles between OAR1 age regions persist to depths >170 km, which is inconsistent with the velocity structure required by the standard plate cooling models and Chablis models, and deeper than existing definitions of the bottom of the LAB (Revenaugh and Jordan, 1991; Beghein et al., 2009; Schmerr, 2012; 106 Kawakatsu et al., 2009; Olugboji et al., 2016; Rychert and Shearer, 2009; Naif et al., 2013; Burgos et al., 2014). We also concluded that our S travel time curves are consistent with a half- space cooling thermal model for Earth’s asthenosphere and the seafloor flattening observed in oceanic bathymetry beyond 70 my is not observed in our oceanic travel time curves . Our results support a coupling of the lithosphere and asthenosphere for times scales on the order of the age of the oceans, and while small-scale convection may still create heterogeneities in the oceanic upper mantle, large-scale flow advects these small-scale heterogeneities along with the plates, allowing the upper part of the asthenosphere to continue cooling linearly on average with the square-root of age >200 my. This indicates that the upper part of the oceanic asthenosphere participates in the cooling that forms the kinematic plates (tectosphere). The dominance of this large-scale horizontal flow may be related to the high stresses associated with its channelization in a thin (~100 km) asthenosphere, as well as the possible focusing of the subtectospheric strain in a low-viscosity channel immediately above the 410-km discontinuity (Toffelmier and Tyburczy, 2007; Youngs and Bercovici, 2009; Karato, 2011). This observed thermal aging of oceanic asthenosphere is inconsistent with the LAB hypothesis, which states that lithospheric plates are decoupled from deeper mantle flow by a shear zone in the upper part of the asthenosphere. The ocean basins show convergence depths much greater than the depth of the oceanic LAB. 107 3 Testing the Isopycnic Hypothesis Using Garnet Lherzolite Xenoliths 3.1 Abstract This study examined a large dataset of garnet lherzolites for consistency with the isopycnic hypothesis. Temperature, pressure, and density values were calculated for this data based on available chemical analyses, and a method was determined to allow density calculation based on the Mg# of olivine present in the rocks. The data was analyzed using a non-parametric Theil method and assuming a linear trend for the data. The result of the Theil test allows us to reject the null hypothesis that the slope of the depth-density line is zero (i.e. density does not change with depth), and an isopycnic line calculated for the 1350° adiabat appears consistent with the slopes estimated by the Theil method within approximately a 95% confidence interval. If the isopycnic hypothesis holds true, it puts constraints on the formation cratonic tectosphere. The tectosphere must have reached approximately isopycnic conditions at the time of large-scale tectonic stabilization, and the tectosphere and crust must have stabilized nearly simultaneously. The isopycnic hypothesis also favors the creation of continental tectosphere by advective thickening rather than plume or multi-stage formation models, consistent with the petrological constraint that the garnet lherzolites were depleted at shallow depths, above the garnet stability field, and advected downward. This implies that the formation of continents consists of the convective flow advectively thickening the tectosphere, interspersed with periods of disruption by small scale double-diffusive instabilities, which in some cases, lead to rifting and drift of the continent. Over time, these mechanisms stabilize the tectosphere in isopycnic balance. 108 3.2 Introduction Continental cratons are characterized by their stability over long periods of time (>2.5 by), as well as their low temperature relative to ocean basins, inferred from heat flux measurements (Nyblade and Pollack, 1993) , seismic data (Dziewonski, 1970; Jordan, 1975) , and xenolith data (Jordan, 1978a). This long period of stability is contrary to the expected behavior of dense, cold materials, which are expected to subside, causing cratons to form basins with thick sediment cover, however, continental cratons tend to be relatively high topographically and have little sediment cover. (Jordan, 1975) noted that body wave travel times averaged about 5s greater for stable continents than for oceans, showing that temperature differences between shields and ocean basins extend deep into the mantle. He then hypothesized that these temperature gradients were stabilized by compositional gradients in the underlying rocks, the isopycnic hypothesis. (Jordan, 1978a) further refined the hypothesis stating that this compositional gradient could be generated by the simple depletion of mantle material in basaltic components (Fe and Al) relative to Mg. The isopycnic hypothesis, which describes the relationship between the density of cratonic mantle rocks and their temperatures as a function of depth is shown in equation (3-1). 𝛥𝜌 𝑖𝑠𝑜 ( 𝑧 )= ρ 0 α 𝛥𝑇 ( 𝑧 ) (3-1) Δρ iso (z) is the density difference as a function of depth, while ρ 0 is the density of undepleted mantle material. α is the coefficient of thermal expansion and ΔT(z) is the temperature difference between the continental geotherm and the oceanic adiabat as a function of depth. The isopycnic hypothesis appears to be supported by some authors, i.e.(Doin and Fleitout, 1996; 109 Shapiro et al., 1999a; Lee, 2003; Carlson et al., 2004; James et al., 2004a; Kaban et al., 2004), while other authors present arguments as to the degree of depletion, i.e.(Gaul et al., 2000; Poudjom et al., 2001b; Kaban et al., 2003; Kelly et al., 2003; Artemieva, 2009; Schutt and Lesher, 2010). The strong version of the isopycnic hypothesis further states that at depths >100 km, isopycnic balance is maintained on length scales less than the maximum tectospheric thickness (Jordan, 1978a; Jordan, 2010). This implies a tectosphere with more depleted (thus less dense) peridotites at shallower depths, which grades to less depleted (therefore denser) peridotites at deeper depths. The validity of the isopycnic hypothesis places important constraints on the structure and evolution of cratons. For instance, the three main theories of craton formation (plume origin, underthrusting and imbrication of oceanic lithosphere, and accretion and orogenic thickening of arcs/continental collisions(Schmitz et al., 2004)) would result in differing density-depth profiles for the resulting cratons (Abbot, 1991; Griffin et al., 2003; Parman et al., 2004; Schoene et al., 2008; Arndt et al., 2009; Lee et al., 2011). Since the tectosphere must have achieved an approximately isopycnic state at the time of large-scale tectonic stabilization, single stage and plume models don’t work (Jordan, 2010). If the isopycnic hypothesis holds, the most likely method of craton formation would be the accretion and orogenic thickening of already depleted material at island arcs and areas of continental collision. It also implies that the tectonic stabilization of the crust and tectosphere occurred at approximately the same time (Pearson, 1999; Artemieva, 2006). Seismic, thermal, gravimetric, and petrologic data have all been utilized to test this hypothesis, both directly and indirectly with varying results. Part 1 of this thesis studied the seismic S-velocity (v S ) profiles for stable continental regions and found that their differential (v S ) profiles with respect to mature oceanic regions converged roughly monotonically with depth 110 below ~100 km and those statistically significant differences in v S between these continental regions and mature oceans were maintained to depths of ~340 km. (Doin et al., 1996) looked at geoid anomalies and observed that the geoid over the cratons was the same as that for young (50 my) oceans. These results supported the isopycnic hypothesis, since, if the cratons were not density compensated, a geoid low would be expected to be observed under the cratons. (Pari and Peltier, 1996) utilized the geoid and free-air gravity constraints to assert that, rather than neutrally buoyant, chemically distinct material, there is a cold, high density down-welling of convecting mantle beneath the Laurentian craton. Alternatively, (Shapiro et al., 1999b) , found that shields and platforms have neither geoid nor free-air gravity signals significantly different from zero. (Kaban et al., 2003; Kaban et al., 2004) used similar data to conclude that the density variations due to temperature were only partially (40%) compensated by compositional differences. Artemieva (2009b) and Deschamps et al., (2002) utilized combined global seismic tomography and thermal models in an attempt to distinguish compositional from thermal differences in the continental lithospheric mantle. Artemieva, (2009b)’s results showed that, beneath cratons, the compositional component of seismic velocity anomalies gradually decreases with depth until the compositional component disappears, consistent with the isopycnic hypothesis’ decreasing depletion with depth. Artemieva (2009b) also argues that, under Archaean cratons, the base of the seismic lithosphere correlates with the base of thermal lithosphere, and the chemical boundary layer beneath cratons. Deschamps et al. (2002) used combined gravity data and tomographic models to conclude that down to 300 km depth, the continent beneath old cratons and stable platforms is cold and depleted in Fe, however, the continental roots are stable due to strongly temperature dependent mantle viscosity. 111 Poudjom et al., (2001b) and Gaul et al., (2000) plotted xenocryst derived density vs. depth and concluded, that in general, Mg# decreases with depth, but the density variation is not smooth, but rather shows a rapid increase in density over a short vertical distance in the range of 180-220 km, deviating from a strict interpretation of the isopycnic hypothesis. Kelly et al., (2003b) utilized a suite of garnet peridotite xenoliths in order to estimate their density and conclude that rather than being neutrally buoyant, as predicted by the isopycnic hypothesis, the cratonic lithosphere is positively buoyant, and this positive buoyancy is offset by dense layers in the cratonic crust or upper mantle which counteract this positive buoyancy leading to overall isopycnic balance. In contrast, Lee, (2003b), in a similar study utilizing bulk-rock Mg#, found that the density of cratonic mantle balances to within error of the negative thermal buoyancy induced by its cooler state relative to asthenospheric mantle. James et al., (2004b) analyzed a large (120) sample of xenoliths from the Kaapvaal craton and concluded that the density of hot primitive mantle does not differ significantly from that of cold cratonic mantle. (Schutt and Lesher, 2006; Schutt and Lesher, 2010) argue that for the Kaapvaal craton, the excess density of the Kaapvaal mantle is not completely compensated by melt depletion and that the cratonic stability is more strongly related to high viscosity, rather than density effects. (Forte and Perry, 2000) inverted geodynamic data and tomography-based flow models to conclude that the approximate equilibrium between thermal and chemical buoyancy contributes to cratonic stability. This lack of a consensus regarding the validity and/or the extent of isopycnic compensation highlights the need for formal testing of this hypothesis. As a first step, this study utilizes a formal statistical analysis in order to test the temperature-density relationship for a large set of garnet lherzolites for several cratons at depths between 100-220 km. This initial test evaluates 112 whether there is a positive trend to the density-depth relationship in sub-cratonic mantle peridotites, vs. the null hypothesis that there is no trend (i.e. no depth-density relationship). Initially, temperatures and depths of equilibration were calculated for the xenoliths, and density values were then obtained for each of the samples. Statistical tests were performed on the depth- density relationships for these samples in order to test their consistency with the isopycnic hypothesis. 3.3 Data Set The xenolith data for this study were compiled from four major sources, (Nimis and Grütter, 2010) ’s Go-Forward data set, (Baptiste et al., 2012) ’s xenolith data, a dataset provided by (McDonough, 2013) , as well as from a list of peridotite bulk compositions from the supplementary material to (Canil and Lee, 2009). This study used 4-phase garnet lherzolites and garnet harzburgites because they represent the dominant peridotites present at the depths of interest for studying the isopycnic hypothesis (Ringwood, 1975; Yoder, H. S., 1976). I avoided spinel-bearing samples, since at the depths of interest (>100 km), garnet is the stable Al phase, and under these conditions, spinel’s presence likely indicates sample disequilibrium. All pyroxenites and eclogites from the data set were also eliminated, as they are likely high density cumulates and not representative of the mantle for the isopycnic hypothesis (Pearson et al., 2003; Carlson et al., 2005). Xenoliths with thermobarometric depths equivalent to < 100 km or 825° C (corresponds to ~ 100 km on the 40 mW geotherm) were eliminated, as this is shallower than the depths of interest. All samples were required to have electron microprobe oxide percentage results for orthopyroxene, clinopyroxene, olivine, and garnet which allowed calculation of pressure and temperature by various methods. Since pure harzburgites lack clinopyroxene and garnet, I was 113 unable to calculate pressure and temperature for them, so my results may be biased towards the higher density lherzolites. However, the data set has a large representation at all depths of interest, and the harzburgites which are present in the dataset fall along roughly the same P-T curve as the lherzolites, so this bias is thought to be small. The xenoliths were all from kimberlite sources, and were also required to have either whole rock oxide percentage data or olivine Mg# data in order to allow density calculations. The kimberlite localities were all required to be on-craton, rather than craton margin in order to attempt to sample true cratonic mantle. In the majority of cases, whole rock oxide information was not available, so this paper details a method utilizing whole rock compositions and olivine Mg#’s to estimate xenolith density. The dataset was filtered using the method of (Nimis and Grütter, 2010). They used a large database of garnet peridotites to examine internal consistency between geothermobarometers. They then used the correlation between various thermobarometers to examine the dataset for “bad” samples. By comparing the TA98, (Taylor, 1998) orthopyroxene (opx)-clinopyroxene (cpx), two mineral thermometer with the BKN90, (Brey et al., 1990) Ca in opx single mineral thermometer, and looking at their correlation, they were able to eliminate samples that may have a departure from opx-cpx equilibrium or have bad chemical analyses. They also used the correlation between the TA98 opx-cpx thermometer with (Nimis and Grütter, 2010)’s cpx-Gt thermometer in order to eliminate samples which may have decoupling of Fe-Mg and Ca-Mg equilibria due to thermal perturbations, redox conditions inconsistent with unperturbed mantle, or poor chemical analysis. Artemieva, (2009b) and Griffin et al., (2009) argue that data from kimberlite suites, often used to test the isopycnic hypothesis, are metasomatized and have artificially reduced degrees of depletion. This method of data selection addresses their concerns 114 that kimberlite xenoliths do not represent undisturbed cratonic mantle by removing those xenoliths which show apparent disequilibrium and metasomatism. The final data set consists of 326 samples from the Kaapvaal, Zimbabwe, Slave, Rae, and Superior cratons. 3.4 Pressure and Temperature Temperature was calculated using several geothermometers, including both single pyroxene (cpx), NT00, (Nimis and Taylor, 2000), and two pyroxene (cpx-opx), BKN90, (Brey et al., 1990), TA98, (Taylor, 1998) methods. Various geobarometers for garnet peridotites, including Al in opx, MC74, (MacGregor, 1974), NG85, (Nickel and Green, 1985), BKN90, (Brey et al., 1990), as well as single cpx based NT00, (Nimis and Taylor, 2000) were used in conjunction with these geothermometers. In order to choose a best temperature/pressure calculation method, I used the xenoliths with existing carbon phases in order to eliminate those methods which put the diamond or graphite phases in the incorrect stability field. I compared the pressure and temperature results for both diamond and graphite bearing xenoliths for each method using the graphite-diamond line of (Kennedy and Kennedy, 1976). Nine thermometer–barometer combinations were examined. Figure 3-1 shows the results of this comparison. The combinations of the TA98- NG85, NT00-MC74 and BKN90-MC74 (circles on Figure 3-1), provide results which fall within the correct graphite-diamond stability fields. The other six methods placed several diamond-bearing xenoliths in the graphite stability field, and were eliminated from consideration. Even though the three best thermobarometers all fit the diamond/graphite stability 115 criteria, there is still a temperature spread of over 100° C possible for a given xenolith, as well as a pressure spread of roughly 8 kbar. Figure 3-1 Comparison of T-P calculation methods for carbon bearing xenoliths. Colored circles denote methods which yielded carbon phases in the correct stability field. Gray symbols denote methods which yielded carbon phases in the incorrect stability field. Solid symbols indicate diamond-bearing xenoliths, and hollow symbols indicate graphite-bearing xenoliths. Figure 3-2 shows the dataset plotted using the 3 consistent thermobarometer combinations. For the remainder of this paper, the NT00-MC74 method will be used; however, the statistical test was performed for all three temperature-pressure method calculations and will be discussed in that section. 116 Figure 3-2 Temperature vs. Pressure for the xenolith dataset for three thermobarometer combinations. Solid circles denote non-sheared xenoliths, while hollow circles denote sheared xenoliths. Even though all three methods yield carbon bearing phases in the correct diamond stability field, there is a possible temperature spread of ~100° C and pressure spread of ~8 kbar. For the chosen thermobarometers, the related depth was calculated from pressure values using the pressure/depth relationship from PREM (Dziewoński and Anderson, 1981). For several samples, more than one set of microprobe data was available and comparison of two analyses on a given sample was possible. For these samples, the differences in pressure and temperature calculated for a given method remained small. For example, three individual analyses of a single xenolith sample yielded a temperature difference of 25°C and pressure differences of less than 2 kbar. 117 3.5 Density Two methods were used to estimate density for the xenoliths. A method using whole-rock oxides to calculate normative density described by (Jordan, 1979) , and a method to using olivine’s magnesium number (Mg#) described later in this paper. The whole-rock oxide method requires 4-phase xenoliths with whole rock oxide percentage information. This method uses oxide percentages, apparent distribution coefficients for garnet lherzolite xenoliths, the Ca-Mg ratio in clinopyroxene, and mineral stoichiometry for phases present in these xenoliths to determine normative xenolith densities. This method yields accurate normative densities (Jordan, 1978a), however, the majority of the samples in the data set did not have whole rock oxide information available. In order to increase the usable density data, a density calculation method using the Mg# of olivine in conjunction with the whole rock oxide method is detailed below. 3.6 Ol Mg# as a Proxy for Density Most of the samples which have temperature and pressure data do not have adequate information available to calculate normative densities using the whole rock method. The olivine Mg# is suggested as a proxy for calculating density. Figure 3-3 shows the olivine Mg#’s for a large set of kimberlite source peridotites plotted against density calculated using the whole-rock oxide method. A 2 nd degree polynomial regression of the data was used to estimate the Mg#/density relationship. This relationship was then used to calculate the density of the xenolith data set. This allowed the calculation of an approximate density from the available olivine Mg#’s for those rocks which had no other source of density information and expanded my usable dataset from 90 to 326 samples. The coefficient of determination (R 2 ) for the linear relationship 118 between whole rock oxide density and olivine Mg# shows that the polynomial equation accounts for 71% of the variance in the olivine magnesium numbers and the root mean square error for the model is 0.015 g/cm 3 . This relationship appears to indicate that Olivine Mg# may be a used as a proxy for density, with the understanding that the error of the density values is increased. Figure 3-3 Ol Mg# vs. density for a suite of kimberlite-derived xenoliths, the data shows a roughly 2 nd degree polynomial trend (with scatter). The equation is the result of a 2 nd degree polynomial regression of the data. The coefficient of determination for the regression (R 2 ) is 0.71 and the root mean square error is 0.0015 g/cm 3 . 3.7 Estimated Geotherms Figure 3-4 below shows the depth-temperature relationships for xenoliths on various cratons using the NT00-MC74 geothermobarometer. The plot shows the “standard” geotherm 119 from (Pollack and Chapman, 1977) at the 40 mW/m 3 level for the upper mantle, as well as the 1350° oceanic adiabat . The Kaapvaal data fit the 40 mW/m 3 geotherm well, while the Slave craton tends to be colder and the Rae craton warmer. For the purposes of this study, I use the 40 mW/m 3 geotherm for all of the isopycnic calculations. Figure 3-4 Temperature-Depth relationships for the xenolith dataset by craton. Solid circles denote unsheared samples, while hollow circles denote sheared samples. Individual cratons show distinct P-T relationships. The solid lines denote the 40 mW geotherm and the 1350°C adiabat respectively. 120 3.8 Oceanic Adiabat In order to calculate the oceanic adiabat, the thermal expansivity of the mantle with respect to depth (z) was calculated using the method outlined by (Afonso et al., 2005), which uses results from mineral physics experiments to calculate thermal expansivity at various pressures. The adiabatic temperature gradient 𝜕 𝑇 𝜕𝑟 was then calculated using the equation from (Turcotte and Schubert, 2002) for 1300°, 1350°, and 1400° adiabats. 𝜕𝑇 𝜕𝑟 = −𝑇𝛼𝑔 𝐶 𝑃 (3-2) where T is the temperature in Kelvin, α is the coefficient of thermal expansion, g is acceleration due to gravity, and C p is specific heat , 3.9 Test of Isopycnicity The isopycnic relationship between density, pressure and temperature is calculated for a given depth using the method of (Jordan, 1988) with ρ 0 = 3.397 g/cm 3 and α calculated from (Afonso et al., 2005). Figure 3-5 shows the data set plotted with isopycnic curves created using the (Pollack and Chapman, 1977)’s 40 mW/m 3 geotherm and adiabats calculated at 1300°, 1350°, and 1400° C. 121 Figure 3-5 Relationship between temperature and Ol Mg# based density for the xenolith suite with isopycnic lines for various adiabats. The simple linear regression fit line is shown in black. Figure 3-6 shows the same data and isopycnic curves plotted as a function of depth. The data visually appears to follow an isopycnic pattern; however, a statistical analysis must be performed to verify this observation. 122 Figure 3-6 Density vs. depth for the xenolith suite, with isopycnic lines for various adiabats. 3.10 Theil Slope Estimation Method Figure 3-7 shows the data plotted with a simple linear regression, the 95% confidence intervals on the values of the slope of the regression, and the 1350° C isopycnic line. Although this estimate appears to fit the isopycnic line fairly well, it assumes a Gaussian distribution of the data. The presence of several high density outliers appears to shallow the slope. In order to minimize the effects of these outliers, the Theil method non-parametric test as outlined by (Hollander and Wolfe, 1973). This method estimates the slope by choosing the median slope among all lines through pairs of two dimensional sample points. For this method, the best fit for the data is assumed to be a straight line (Y i ) , modeled as: 123 𝑌 𝑖 = 𝛼 + 𝛽 𝑥 𝑖 + 𝑒 𝑖 , 𝑖 = 1, … , 𝑛 (3-3) where x’s are the known constants (in this case, the depth of the xenolith), and α and β are unknown parameters, and 𝑒 𝑖 ’s are mutually independent and come from the same continuous population. We perform a distribution-free test for the slope coefficient, β. The null hypothesis (H 0 ) is stated as 𝐻 0 : 𝛽 = 𝛽 0 (specified) (3-4) where we specify β 0 = 0 for a constant density data set. We calculate n number of differences (D i ) between the data points and the line specified by the null hypothesis as: 𝐷 𝑖 = 𝑌 𝑖 − 𝛽 0 𝑥 𝑖 = 1, … , 𝑛 (3-5) and let 𝐶 = ∑ 𝑐 ( 𝐷 𝑗 − 𝐷 𝑖 𝑛 𝑖 <𝑗 ) , (3-6) and 𝑐 ( 𝐷 𝑗 − 𝐷 𝑖 )= { 1 𝑖𝑓 𝐷 𝑗 > 𝐷 𝑖 0 𝑖𝑓 𝐷 𝑗 = 𝐷 𝑖 −1 𝑖𝑓 𝐷 𝑗 < 𝐷 𝑖 . (3-7) In other words, for each pair of subscripts (i,j) with i<j, we score one if D j -D i is positive, and minus one if D j -D i is negative. We add up the ones and minus ones and call the sum C. 124 For a one-sided test of H 0 versus the alternative hypothesis 𝐻 1 : 𝛽 > 𝛽 0 (3-8) at the α level of significance, we use the large sample approximation of 𝐶 ∗ = 𝐶 { 𝑛 ( 𝑛 − 1) ( 2𝑛 + 5) /18} 1/2 (3-9) We then test the null hypothesis using the normal theory approximation for large datasets where we 𝑟𝑒𝑗𝑒𝑐𝑡 𝐻 0 𝑖𝑓 𝐶 ∗ ≥ 𝑧 𝛼 , 𝑎𝑐𝑐𝑒𝑝𝑡 𝐻 0 𝑖𝑓 𝐶 ∗ > 𝑧 𝛼 , (3-10) and z a is the solution of the integral equation of the normal distribution above the desired probability, α. Using a value of α=0.05, β 0 =0, and z a = 1.6449, we calculate C * = 8.9493. According to Equation (3-10), this value of C allows us to reject the null hypothesis that the data set has a slope of 0 (constant density). An estimator for the slope of the data, β, was calculated using the Theil method where we form the N sample-pair slope values 𝑆 𝑖𝑗 = ( 𝑌 𝑖 − 𝑌 𝑗 ) ( 𝑥 𝑗 − 𝑥 𝑖 ) , 𝑓𝑜𝑟 𝑖 < 𝑗 ⁄ (3-11) Then the estimator of β is 𝛽 ̂ = 𝑚𝑒𝑑𝑖𝑎𝑛 { 𝑆 𝑖𝑗 } (3-12) Let S (1) ≤ ... ≤ S (N) denote the ordered values of S ij , For an odd value of N, let N=2k+1, and 𝛽 ̂ = 𝑆 ( 𝑘 +1) (3-13) 125 for an even value of N, let N=2k, and 𝛽 ̂ = 𝑆 ( 𝑘 ) + 𝑆 ( 𝑘 +1) 2 . (3-14) For the data set, we calculate a slope of 𝛽 ̂ = 3.6619 x 10 -4 , which is plotted on Figure 3-8. The Theil method is then also used to calculate a two-sided confidence interval for the slope, β, using the large-sample approximation where 𝐶 𝛼 ≈ 𝑧 ( 𝑎 2 ⁄ ) { 𝑛 ( 𝑛 − 1) ( 2𝑛 + 5) 18 } 1/2 (3-15) and 𝑀 1 = 𝑁 − 𝐶 𝛼 2 (3-16) 𝑀 2 = 𝑁 + 𝐶 𝛼 2 (3-17) The 1-α confidence intervals (β L, β U ) contain the middle 95% of the slopes of the lines determined by the pairs of points. These are given by 𝛽 𝐿 = 𝑆 ( 𝑀 1 ) (3-18) and 𝛽 𝑈 = 𝑆 ( 𝑀 2 +1) (3-19) which results in (β L, β U ) = (2.9208 x 10 -4 , 4.3795 x 10 -4 ), as shown on Figure 3-8. The 1350° C isopycnic line is consistent with the Theil-Sen estimated slope within roughly a 95% confidence interval. By using a non-parametric method to estimate the slope which minimizes the effect of outliers, we get a steeper slope and a better fit to the isopycnic line. 126 The Theil method was also employed using two additional temperature-pressure calculation methods (BKN90-MC74 and TA98-NG85) which resulted in slopes of 𝛽 ̂ = 5.2701 x 10 -4 and 𝛽 ̂ = 7.2595 x 10 -4 , respectively. In both cases the null hypothesis was rejected in favor of the alternative hypothesis that the slope was not equal to zero, however BKN90-MC74 better fits the 1400˚C isopycnic line and the TA98_NG85 Theil estimated slope is too shallow to fit any of the isopycnic lines well. Figure 3-7 Dataset plotted with linear regression line and 95% confidence intervals on the slope of the linear regression. The 1300, 1350, and 1400° C isopycnic lines are also shown. The slope of the regression appears to be artificially shallowed by the presence of several high density outliers. 127 Figure 3-8 Density-Depth relationship for xenolith dataset plotted with the estimated slope and its 95% confidence interval calculated by the Theil method. The green line shows the 1350 °C isopycnic curve, which is consistent with the Theil estimate within the confidence interval. 3.11 Discussion We examined a large dataset of garnet lherzolites for consistency with the isopycnic hypothesis. We carefully filtered our data set using the method of Nimis and Grütter, 2010, to remove samples which may have chemical disequilibrium, bad chemical analysis, or have been metasomatized. We used the diamond stability field to choose consistent three 128 thermabarometers that allowed our diamond and graphite-bearing lherzolites to fall within the correct stability fields. Temperature, pressure, and density values were then calculated for this data based on available chemical analyses, and a method was determined to allow density calculation based on the Mg# of olivine present in the rocks. The data was analyzed using the non-parametric Theil method and assuming a linear trend for the data. The result of the Theil test allows us to reject the null hypothesis that the slope of the depth-density line is zero (i.e. density does not change with depth) for all three thermobarometers, and an isopycnic line calculated for the 1350° C adiabat fits within the 95% confidence intervals of the slopes estimated by the Theil to depths of approximately 180 km for our NT00-MC74 thermobarometer. The isopycnic line is slightly outside of the 95% confidence intervals below this depth, this is likely due to the Theil method’s assumption of a linear data trend, while the isopycnic line is slightly non-linear. The BKN90-MC74 thermobaromater also shows a linear trend similar to that of an isopycnic line, however it better fits an isopycnic line calculated for ~1400˚ C. The TA98-NG85 thermobarometer also has a non-zero slope, but does not fit the isopycnic lines well, it’s slope is much shallower than the isopycnic lines. The depth-density plots for all three thermobarometers are shown in Appendix G. This thermobarometer combination, although consistent by the criteria established earlier in this paper, also has the smallest range of pressures calculated for our xenolith suite. It has no xenoliths present at <40 kbar, while the other two methods have a large number of xenoliths present in the 30-40 kbar pressure range. Given the large xenolith suite, and the unlikeliness that we would not have samples of lherzolites from the shallowest depths, our confidence in this thermobarometer is less than in the other two. 129 The lherzolites in this study appear to follow an linear depth-density relationship close to that predicted by our calculated isopycnic lines for reasonable mantle adiabatic temperature. If this is true, in order to achieve isopycnic balance, the tectosphere must have reached approximately isopycnic conditions at the time of large-scale tectonic stabilization and the crust must have stabilized nearly simultaneously (Griffin et al., 1999; Artemieva, 2006) since the depletion of the lithosphere in basaltic components creates a low density peridotite that resists remixing into the mantle (Jordan, 1978a). This supports the creation of continental tectosphere by advective thickening rather than plume or multi-stage formation models, consistent with the petrological constraint that the garnet lherzolites were depleted at shallow depths, above the garnet stability field, and advected downward (Schutt and Lesher, 2006; Lee et al., 2011) and implies that the formation of continents consists of the convective flow advectively thickening the tectosphere, interspersed with periods of disruption by small scale double-diffusive instabilities, which in some cases, lead to rifting and drift of the continent. Over time, these mechanisms stabilize the tectosphere in isopycnic balance (Jordan, 1978a; Jordan, 2010). 130 4 Appendices 4.1 Appendix A: Calculation of Profile Variance For each of the models in Table 1-1, we discretized the aspherical variation v S (z, ) on a 1º  1º geographic grid, and we calculated the radial functional of interest, F( ), using a depth discretization appropriate to the model resolution. This generated an array {F k } of 180  360 = 64,800 elements, each corresponding to a surface area , where  k is the colatitude at the center of the k th 1º  1º element. The area of the X th GTR1 region is 𝐴 𝑋 = ∑ 𝑎 𝑘 𝑘𝑒𝑋 (A1) where the summation is taken over the K X elements of region X. Under this discretization, we approximated the regional averages in (2) by 𝑓 𝑋 = 1 𝐴 𝑋 ∑ 𝑎 𝑘 𝐹 𝑘 𝑘𝑒𝑋 (A2) For the Fisher-von Mises covariance function, the discretized version of equation (4) is 𝜎 𝑋 2 = 𝘀 𝑋 2 𝐴 𝑋 2 ∑ ∑ 𝑎 𝑘 𝑎 𝑘 ′ 𝑒 𝜅 𝑋 ( 𝑐𝑜𝑠 𝛹 𝑘 𝑘 ′ −1) 𝑘 ′𝑒𝑋 𝑘𝑒𝑋 . (A3) We calculated the local variance using the unbiased estimator 𝘀 ̃ 𝑋 2 = 𝐾 𝑋 ( 𝐾 𝑋 −1) 𝐴 𝑋 ∑ 𝑎 𝑘 ( 𝐹 𝑘 − 𝑓 𝑋 ) 2 𝑘𝑒𝑋 . (A4) Estimates of the dispersion parameters  X were obtained from a relation that we derived between  X 2 and a variance computed from 5º  5º averages of the radial functional. Defining R Xn a k » (p /180) 2 sinq k 131 to be the n th 5º  5º cell of the X th region, we computed the average of F k over the 25 elements of this cell: 𝑓 𝑋 𝑛 = 1 𝑎 𝑋 𝑛 ∑ 𝑎 𝑘 𝐹 𝑘 𝑘𝑒 𝑅 𝑋 𝑛 , where 𝑎 𝑋 𝑛 = ∑ 𝑎 𝑘 𝑘𝑒 𝑅 𝑋 𝑛 . (A5) The variance of this average is the analog of equation (A3): 𝜈 𝑋 𝑛 2 ≡ 〈( 𝐹 𝑘 − 𝑓 𝑋 𝑛 ) 2 〉 = 𝘀 𝑋 2 𝑎 𝑋 𝑛 2 ∑ ∑ 𝑎 𝑘 𝑎 𝑘 ′ 𝑒 𝜅 𝑋 ( 𝑐𝑜𝑠 𝛹 𝑘 𝑘 ′ −1) 𝑘𝑒 𝑅 𝑋 𝑛 𝑘𝑒 𝑅 𝑋 𝑛 . (A6) According to the assumptions associated with equation (A5), the variances are the same for all N X cells in the X th region, so that an unbiased estimator for the variance of the 5º  5º averages is the sample variance, 𝜈 ̃ 𝑋 2 = 𝐾 𝑋 ( 𝐾 𝑋 −1) 𝐴 𝑋 ∑ 𝑎 𝑋 𝑛 ( 𝑓 𝑋 𝑛 − 𝑓 𝑋 ) 2 𝑁 𝑋 𝑛 =1 . (A7) We invoked the Gaussian approximation of equation (12) and represented each 5º  5º cell as a Cartesian square of 1º  1º elements, which reduced (A6) to 𝜈 𝑋 2 = 𝘀 𝑋 2 ( 25) 2 ∑ ∑ e −𝛹 𝑘 𝑘 ′ 2 /2𝜆 𝑋 2 25 𝑘 ′ =1 25 𝑘 =1 (A8) We renumbered the elements of the cells into five rows (i = 1-5) and five columns (j = 1-5) and used the Cartesian approximation , where  kk' is now expressed in degrees. Equation (A8) can thus be written as the product of two identical double sums. Taking the square-root reduced it to 𝜈 𝑋 = 𝘀 𝑋 25 ∑ ∑ 𝑒 −( 𝑖 −𝑖 ′ ) 2 /2𝜆 𝑋 2 5 𝑖 ′ =1 5 𝑖 =1 . (A9) Equation (A9) was rearranged by summing along the diagonals of the 5  5 matrix to obtain n Xn 2 y kk' 2 = (i-i') 2 +(j-j') 2 132 𝜈 𝑋 = 𝘀 𝑋 25 ( 5 + 2 ∙ 4e − 1 2𝜆 𝑋 2 + 2 ∙ 3e − 4 2𝜆 𝑋 2 + 2 ∙ 2e − 9 2𝜆 𝑋 2 + 2 ∙ 1e − 16 2𝜆 𝑋 2 . (A10) The general form for a square domain of M 2 elements is 𝜈 𝑋 = 𝘀 𝑋 𝑀 2 [𝑀 + 2 ∑ ( 𝑀 − 1) 𝑒 −𝑚 2 /2𝜆 𝑋 2 𝑀 −1 𝑚 =1 ]. (A11) When  X << 1, equation (A11) gives  = M , which is the s.e.m. of M 2 completely independent samples, and when  X >> 1, it yields  = , which is the s.e.m. of completely dependent samples. To calculate  X in units of angular degrees, we solved equation (A10) using the estimates from equation (A4) and from equation (A7). Owing to the areal averaging in the estimator (A7), which appropriately down-weights the distorted cells near the poles, the Cartesian relation (A10) is an excellent approximation. Profiles of {  X (z) : X = A, B, C, Q, P, S} were computed for each model, and we used these profiles to estimate the s.e.m. profiles,  X (z), from equation (A3). 4.2 Appendix B: Regional Profiles The figures in this appendix display the regional profiles {f X = /v S (z) : X = A, B, C, Q, P, S; 50 km ≤ z ≤ 2891 km} for each of the global tomographic models in Table 1-1, color-coded by region and organized by model group. The profiles have been normalized to a relative velocity anomaly, expressed in percent, using the spherically symmetric model assumed in the inversion of the aspherical model. All profiles are banded with their 1  uncertainties, computed from equation (A3). e X n X dv S X (z) 133 Figure B-1 GTR1 projections of the relative aspherical variations v S /v S for the Texas models, banded with their 1  uncertainties. The left panel, for the Texas exemplar TX2008, is reproduced in Figure 4. 134 Figure B-2 GTR1 projections of the relative aspherical variations v S /v S for the Harvard models, banded with their 1  uncertainties. The left panel, showing the Harvard exemplar S362ANI, is reproduced in Figure 4. 135 Figure B-3 GTR1 projections of the relative aspherical variations v S /v S for the Scripps models, banded with their 1  uncertainties. The left panel, showing the Scripps exemplar HMSL-S06, is reproduced in Figure 4. 136 Figure B-4 GTR1 projections of the relative aspherical variations v S /v S for the Berkeley models, banded with their 1  uncertainties. The left panel, showing the Berkeley exemplar SAW642ANB, is reproduced in Figure 4. 137 Figure B-5 GTR1 projections of the relative aspherical variations v S /v S for models by Ritsema et al. (2004, 2010) banded with their 1  uncertainties. The left panel, showing the Ritsema exemplar S20RTS, is reproduced in Figure 4. 138 Figure B-6 GTR1 projections of the relative aspherical variations v S /v S for models by Takeuchi (2007, 2012) and Auer (2014) banded with their 1  uncertainties. The bottom panel, showing the Auer exemplar SAVANI, is reproduced in Figure 4. 4.3 Appendix C: Bias and Variance of Regional Profiles If is the true (“real-Earth”) value of the aspherical function at position r = (z, ), then the local model error is F(r) – . For the tomographic models considered here, the forward calculation of the data functionals was based on linearized perturbation theory, and the inversion process can therefore be approximated as a linear filter with noise: ˆ F(r) ˆ F(r) 139 𝐹 ( 𝐫 )= ∫ 𝐺 ( 𝐫 ; 𝐫 ′ ) 𝐹 ̂ ( 𝐫 ′ ) 𝑑 𝐫 ′ + N( 𝐫 ) 𝑽 . (C1) Here G(r; r') is the kernel of the inversion filter, N(r) is additive noise arising from data errors, and the integral is taken over the volume of the mantle. Within a tectonic region X, can be written as the sum of its mean value and a structural variation that averages to zero over X. Applying the projection operator (2) to equation (C1) yields (C2) The first summation involves bi-regional averages of the inversion kernel: . (C3) In a well-conditioned inversion of data sets with good global coverage, G(r; r') will be a localized function with a radial smoothing length comparable to the mean correlation length and an angular smoothing length comparable to the mean correlation angle . For most of the tomographic models in Table 1-1, the latter is substantially smaller than the size of an individual tectonic region (see §6). The angular averaging in (C3) thus ensures that the ratio g XY /g XX will be small for Y  X, allowing us to drop these terms from the first summation in (C2). The second summation involves the structural fluctuations , which average to zero over each region, so the bi-regional averages for Y  X will contribute even less. Ignoring these cross terms ˆ F(r) ˆ f X (z) = A X -1 ˆ F(z,W)dW X ò ˆ E X (r) = ˆ F(r)- ˆ f X (z) f X (z) = g XY (z,z') ˆ f Y (z') dz' ò Y å + 1 A X G(z,W;z',W') ˆ E Y (z',W') dz'dW' + N(z,W) ò Y ò Y å é ë ê ù û ú X ò dW g XY (z,z') º 1 A X A Y G(z,W;z',W') dW Y ò X ò 'dW r(z) l(z) ˆ E Y 140 is the same approximation that justifies the decorrelation condition (equation 11), and it leads to a strictly regional relationship between the model profile and the true Earth structure: 𝑓 𝑋 ( 𝑧 )≈ ∫𝑔 𝑋𝑋 ( 𝑧 , 𝑧 ′ ) 𝑓 ̂ 𝑋 ( 𝑧 ′ ) 𝑑 𝑧 ′ + 𝑛 𝑋 ( 𝑧 ) . (C4) The additive noise in this expression comprises regional averages of the inverse-filtered structural fluctuations and the inversion noise: . (C5) The expected value of the first term is nearly zero by construction, and that of the second term will be zero if the data used in the inversion are unbiased. Therefore, under the reasonable assumption that 〈𝑛 𝑋 ( 𝑧 ) 〉 ≈ 0, the model error can be described by a profile bias, 𝛽 𝑋 ( 𝑧 )≡ 〈𝛷 𝑋 ( 𝑧 )− 𝑓 ̂ 𝑋 ( 𝑧 ) 〉 = 𝜙 𝑋 ( 𝑧 )− 𝑓 ̂ 𝑋 ( 𝑧 )≈ ∫ [𝑔 𝑋𝑋 ( 𝑧 , 𝑧 ′ )− 𝛿 ( 𝑧 − 𝑧 ′ ) ]𝑓 ̂ 𝑋 ( 𝑧 ′ ) 𝑑𝑧 ′ , (C6) and a profile variance, 𝜎 𝑋 2 ( 𝑧 )≡ 〈[𝛷 𝑋 ( 𝑧 )− 𝜙 𝑋 ( 𝑧 ) ] 2 〉 ≈ 〈𝑛 𝑋 2 ( 𝑧 ) 〉. (C7) An estimator of (C7) is the intra-regional sample variance, as computed in Appendix A. Under the decorrelation condition (11), the variance of a differential profile f X-Y  f X – f Y can be approximated as the sum of the profile variances: 𝜎 𝑋 −𝑌 2 ( 𝑧 )≡ 〈[𝛷 𝑋 −𝑌 ( 𝑧 )− 𝜙 𝑋 −𝑌 ( 𝑧 ) ] 2 〉 ≈ 𝜎 𝑋 2 ( 𝑧 )+ 𝜎 𝑌 2 ( 𝑧 ) . (C8) The bias in the differential profile can be written 𝛽 𝑋 −𝑌 ( 𝑧 )≡ 〈[𝛷 𝑋 −𝑌 ( 𝑧 )− 𝑓 ̂ 𝑋 −𝑌 ( 𝑧 ) ] 〉 = 𝜙 𝑋 −𝑌 ( 𝑧 )− 𝑓 ̂ 𝑋 −𝑌 ( 𝑧 )≈ ∫ [𝑔 𝑋 −𝑌 ( 𝑧 , 𝑧 ′ )− 𝛿 ( 𝑧 − 𝑧 ′ ) ][𝑓 ̂ 𝑋 ( 𝑧 ′ )− 𝑓 ̂ 𝑌 ( 𝑧 ′ ) ]𝑑 𝑧 ′ + ∫𝑔 𝑋 +𝑌 ( 𝑧 , 𝑧 ′ ) [𝑓 ̂ 𝑋 ( 𝑧 ′ )+ 𝑓 ̂ 𝑌 ( 𝑧 ′ ) ]𝑑𝑧 (C9) n X (z) = 1 A X G(z,W;z',W')E X (z',W') dz'dW'dW + 1 A X ò X ò X ò N(z,W) dW X ò 141 where g X-Y  (g XX + g YY )/2, (C10a) g X+Y  (g XX – g YY )/2. (C10b) If the smoothing kernels for the two regions are similar, which we expect for tomographic models derived from dense global data sets, then the second integral will be much smaller than the first and can be dropped. We make this approximation in our bias analysis (equation 13). 4.4 Appendix D: A Bias Model for Vertical Smearing We represent the bi-regional average of inversion kernel, g X-Y (z, z'), as a Gaussian function with an upper half-width (z') and define to be the value of this half-width at the true convergence depth . We suppose that between this depth and some minimum depth of consideration, z 0 , the true differential profile can be approximated by a function of the form: , (D1) where  is a positive real number. If , then , and the integration in equation (15) gives 𝛽 ̂ 𝑋 −𝑌 ≡ 𝛽 𝑋 −𝑌 ( 𝑧 ̂ 𝑋 −𝑌 )= 𝜇 𝜈 𝑓 ̂ 𝑋 −𝑌 ( 𝑧 ̂ 𝑋 −𝑌 − 𝛾 ) . (D2) Here is the integral, 𝜇 𝜈 = 1 √ 2𝜋 ∫ 𝑧 𝜈 e −𝑧 2 /2 𝑑𝑧 ≈ 2 𝜈 2 −1 √ 𝜋 Γ ( 𝜈 +1 2 ) 𝑧 ̂ 𝑋 −𝑌 −𝑧 0 0 . (D3) ˆ g ˆ z X-Y ˆ f X-Y (z) = a( ˆ z X-Y -z) n ˆ z X-Y - ˆ g >z 0 a= ˆ g -n ˆ f X-Y (ˆ z X-Y - ˆ g) m n 142 Figure 8 illustrates the geometrical relationships for  = 1. In this linear case, the gamma- function expression in (D3) is an adequate approximation when the integration range is just a few times larger than . The tomographic data sets provide good integral constraints on regional upper-mantle structure, such as the vertical shear-wave travel times. In a well-conditioned inversion, the smoothing kernel should be nearly unimodular in z (Backus and Gilbert, 1968) and therefore preserve the integral 𝛹 ̂ 𝑋 −𝑌 = ∫𝑓 ̂ 𝑋 −𝑌 ( 𝑧 ) 𝑑𝑧 , as expressed in equation (14). Using this constraint to solve for the constant a in equation (D1) gives . (D4) We apply the linear version (  𝜈 = ) to obtain equation (15) and the convergence depths in Figures 13. 4.5 Appendix E: Regional Radial Correlation Functions The intra-regional covariance function can be generalized to unequal depths: . (E1) We define the regional radial correlation function to be the normalized average of V XX over X: (E2) Here  X is the single-point standard deviation, which appears in equation (8), and ˆ g X -Y g X-Y (z,z') ˆ b X -Y = (n+1)m n ˆ g n ˆ Y X -Y ( ˆ z X -Y -z 0 ) n+1 1/ 2π V XX (z,W;z',W') = (F(z,W)-f X (z)) (F(z',W')-f X (z')) , W, W' Î X R XX (z,z') = 1 A X e X (z)e X (z') V XX (z,W;z',W) dW X ò » 1 A X e X (z)e X (z') E X (z,W) E X (z',W) dW X ò . 143 〈𝐸 𝑋 ( 𝐫 ) 𝐸 𝑋 ( 𝐫 ′ ) 〉 = ∫ ∫ 𝐺 ( 𝐫 , 𝐫 ′′ ) 𝑋 𝑋 〈𝐸 ̂ 𝑋 ( 𝐫 ′′ ) 𝐸 ̂ 𝑋 ( 𝐫 ′′′ ) 〉𝐺 ( 𝐫 , 𝐫 ′′′ ) 𝑑 𝐫 ′′ 𝑑 𝐫 ′′′ (E3) is the inversion-filtered, intra-regional covariance function. If the region X is expanded to include the entire sphere S 1 , then R becomes the global radial correlation function introduced by Jordan et al. (1993). The approximation in the second line of (C12) ignores the contribution of the data noise N, which should be small for the regional scales considered here. At these large scales, the angular variations should also average out, allowing us to write the radial correlation function as the double depth integral, 𝑅 𝑋𝑋 ( 𝑧 , 𝑧 ′ ) ~ ∫∫𝑔 𝑋𝑋 ( 𝑧 , 𝑧 ′′ ) 𝑒 𝑋𝑋 ( 𝑧 ′′ , 𝑧 ′′′ ) 𝑔 𝑋𝑋 ( 𝑧 ′ , 𝑧 ′′′ ) 𝑑𝑧 ′′ 𝑑𝑧 ′′′ (E4) This quadratic form involves g XX and the regional average of the unfiltered intra-regional covariance, (E5) The shape of the correlation function near its diagonal peak is approximated by a Gaussian with radial correlation length  X (z): (E6) If similar approximations are made to the other two functions in (E4), , (E7) , (E8) then this integral yields a quadratic relation among the three characteristic scales: ˆ e XX (z,z') º 1 A X 2 ˆ E X (z,W) ˆ E X (z',W') X ò X ò dW'dW R XX (z,z') ~ e -(z-z') 2 /2r X 2 g XX (z,z') ~ e -(z-z') 2 /2g X 2 ˆ e XX (z,z') ~ e -(z-z') 2 /2h 2 144 . (E9) Here  X is the half-width of the inversion smoothing kernel, and is the vertical scale length of the unfiltered intra-regional variation. In well-conditioned tomographic inversions, the average smoothing lengths for any two large regions will be comparable,  X   Y , so that the inter- regional differences in the radial correlation lengths will be primarily due to differences between and . Therefore, (E9) implies . (E10) This inequality depends on a Gaussian approximation to the regionally averaged inversion kernel (E7). If the kernel has negative sidebands of significant amplitude, then the fluctuations in the radial correlation function due to sideband spillage can violate the inequality, as observed in Figure 11. 4.6 Appendix F: Bayesian Interpretation of Convergence Depth Probabilities Our structural model assumes that > 0 for z < and = 0 for z ≥ . For each trial depth z, we consider three sets of hypotheses: H 0 : = 0; H 1 : > 0. (F1) H 0 ': z = ; H 1 ': z < ; H 2 ': z > . (F2) H 0 ": z ≥ ; H 1 ": z < . (F3) According to the structural model, all three sets are complete, r X 2 (z)= 2g X 2 +h X 2 h X h X h Y g(z) £ 1 2 min r X (z),r Y (z) [ ] ˆ f X -Y (z) ˆ z X-Y ˆ f X -Y (z) ˆ z X-Y ˆ f X -Y (z) ˆ f X -Y (z) ˆ z X-Y ˆ z X-Y ˆ z X-Y ˆ z X-Y ˆ z X-Y 145 P(H 0 ) + P(H 1 ) = P(H 0 ') + P(H 1 ') + P(H 2 ') = P(H 0 ") + P(H 1 ") = 1, (F4) and their members are mutually exclusive; e.g., P(H 0 | H 1 ) = P(H 1 | H 0 ) = 0. Each tomographic estimate M n contributes an observed differential profile for which a bias-corrected p-value can be calculated using equation (18): 𝑝 𝑋 −𝑌 ( 𝑧 |𝑀 𝑛 )= 1 √ 2𝜋 𝜎 𝑋 −𝑌 ∫ exp[− ( 𝑓 𝑋 −𝑌 −𝑏 𝑋 −𝑌 −𝜙 𝑋 −𝑌 ) 2 2𝜎 𝑋 −𝑌 ] 0 −∞ 𝑑 𝜙 𝑋 −𝑌 . (F5) Strictly speaking, this p-value is the probability of conditional on H 0 , which we abbreviate as = P(M n | H 0 ). (F6) In the Bayesian analysis of §6.5, we interpret these p-values to be the probability of z ≥ given the n th model or, in the notation used here, P(H 0 " | M n ). The total probability theorem allows us to write P(H 0 " | M n ) = P(H 0 " | H 0 ) P(H 0 | M n ) + P(H 0 " | H 1 ) P(H 1 | M n ) = P(H 0 | M n ), (F7) because P(H 0 " | H 0 ) = 1 and P(H 0 " | H 1 ) = 0. However, the requisite equality of (F6) and (F7) is not obvious; indeed, the mistaken identification of a conditional probability with its inverse leads to the prosecutor’s fallacy and other famous statistical errors (e.g., Nickerson, 2000). We can demonstrate the equality by appealing to Bayes’s theorem, which in this case states P(H 0 | M n ) = . (F8) f X -Y (z |M n ) f X -Y (z |M n ) p X -Y (z |M n ) ˆ z X-Y P(M n |H 0 ) P(H 0 ) P(M n |H 0 ) P(H 0 )+P(M n |H 1 ) P(H 1 ) 146 P(H 0 ) measures our prior belief in H 0 . We might take P(H 0 ) ~ 1 for a large test depth z (say, in the lower mantle), and P(H 0 ) ~ 0 for a small depth (say, above 100 km), but for any interesting depth, i.e., one approaching , a fair observer should weight the complementary hypotheses equally, P(H 0 ) = P(H 1 ) = ½. (F9) The priors thus cancel out of (F8). P(M n | H 1 ) is the probability of observing – if > 0. Assuming all positive values of are equally probable under H 1 and discounting any error in the bias estimation, then this probability is an integral over the likelihood of , which is assumed to be Gaussian: P(M n | H 1 ) = 1 √ 2𝜋 𝜎 𝑋 −𝑌 ∫ exp[− ( 𝑓 𝑋 −𝑌 −𝑏 𝑋 −𝑌 −𝜙 𝑋 −𝑌 ) 2 2𝜎 𝑋 −𝑌 ] ∞ 0+ 𝑑 𝜙 𝑋 −𝑌 . (F10) Therefore, P(M n | H 1 ) = 1 – P(M n | H 0 ), and P(H 0 " | M n ) = Q.E.D. (F11) We have ignored the uncertainties in the estimation of bias and variance in this analysis; e.g., the likelihood functions in (F5) and (F10) are conditional on b X-Y and  X-Y . A formal statistical analysis would treat both as nuisance parameters and integrate their uncertainties to get the appropriate marginal distributions (e.g., Cox and Hinkley, 1974). Instead, we have lumped them into the epistemic uncertainties, which we have characterized through the sensitivity and ensemble analysis of §6.5. As described in the text, the sensitivity of convergence depths to the bias uncertainty is small compared to the aleatory variability as well as the other epistemic uncertainties indicated by the variations in . An analysis based on equation (4) indicates that the uncertainty in the aleatory variance contributes even less. ˆ z X-Y f X -Y (z |M n ) b X -Y (z |M n ) ˆ f X -Y (z) f X -Y = ˆ f X -Y +b X -Y f X -Y p X -Y (z |M n ) f X -Y (z |M n ) 147 There is one final point to be clarified. Our interpretation of the p-value as P(H 0 " | M n ) relies on P(H 0 " | H 0 ) being unity (equation F7), yet the bias estimate has been conditioned on H 0 ', not H 0 " = H 0 '  H 2 '. Clearly, P(H 0 ' | H 0 ) ≠ 1. However, the bias estimate under H 0 ' is always greater than under H 2 '. This property of our bias model, which is graphically apparent in Figure 8, can be expressed as b X-Y (z | H 0 ', M n ) > b X-Y (z | H 2 ', M n ). It implies p X-Y (z | H 0 ', M n ) > p X-Y (z | H 2 ', M n ) (F12) Therefore, it is more accurate to interpret p X-Y (z | M n ) as an upper bound on P(H 0 " | M n ). 4.7 Appendix G: Density vs. Depth for Thermobarometers 148 Figure H-1 Density-Depth relationship for xenolith dataset plotted with the estimated slope and its 95% confidence interval calculated by the Theil method. The green line shows the 1350 °C isopycnic curve, which is consistent with the Theil estimate within the confidence interval. 149 Figure H-2 Density-Depth relationship for xenolith dataset plotted with the estimated slope and its 95% confidence interval calculated by the Theil method. The green line shows the 1350 °C isopycnic curve, which is consistent with the Theil estimate within the confidence interval. 150 Figure H-1 Density-Depth relationship for xenolith dataset plotted with the estimated slope and its 95% confidence interval calculated by the Theil method. The green line shows the 1350 °C isopycnic curve, which is consistent with the Theil estimate within the confidence interval. 151 List of Figures Figure 1-1 Global tomographic models analyzed in this study, listed by research group. Maps are Mollweide equal-area projections showing the relative aspherical variation in isotropic or Voigt- averaged shear-wave speed at 150 km depth; positive variations are in blue and negative variations in red. See Table 1 for model descriptions and references. ......................................... 12 Figure 1-2 The GTR1 global tectonic regionalization on a 5º  5º geographic grid (Jordan, 1981). Area of each region is percentage of the sphere. .......................................................................... 15 Figure 1-3 Minimum value of the horizontal correlation half-width 𝜆 min plotted against horizontal inner scale  0 for the 27 tomographic models in Table 1. Scaling relation 𝜆 min /  0 = ½ (dashed line) is set by the lower limit of model resolution. .......................................................... 20 Figure 1-4 Profiles of relative aspherical variations v S /v S obtained by projecting the group exemplars onto the six tectonic regions of GTR1 (X = A, B, C, Q, P, S). The profiles are plotted as a function of depth from 50 km to the core-mantle boundary. Profiles for each region (dark lines) were calculated from equation (A2), and the 1  uncertainties (shaded bands) were calculated from equation (A3). Color code is the same as in Figure 1-2. Similar plots for the other 16 tomographic models considered here are given in Appendix B. .................................... 22 Figure 1-5 GTR1 projections of the relative aspherical variations v S /v S for the Texas models, banded with their 1  uncertainties. Same layout as Figure B-1, but truncated at a depth of 700 km to highlight upper mantle structure. Dashed line at 300 km approximates the oceanic. ........ 25 152 Figure 1-6 Differential profiles for oceanic region A relative to B (yellow curves) and oceanic region C relative to B (red curves) for the group exemplars. The curves are banded with their two-sided 95% confidence intervals. The model-specific A-B and C-B convergence depths, listed for all models in Table 1-2, are the minimum depths where 95% confidence bands intersect zero (orange line). ......................................................................................................................... 31 Figure 1-7 Differential shear-velocity profiles for the stable-continent region PS relative to the mature-ocean region BC (blue curves) for the group exemplars. The curves are banded with their two-sided 95% confidence intervals. The model-specific PS-BC convergence depths, listed for all models in Table 2, are the minimum depths where the 95% confidence bands intersect zero (red line). ....................................................................................................................................... 32 Figure 1-8 Schema for the bias calculation using the linear model of equation (17). Left panel: The bi-regional inversion kernel 𝑔𝑋 − 𝑌 centered at the actual convergence depth 𝑧𝑋 − 𝑌 is characterized as a Gaussian with an upper half-width 𝛾𝑋 − 𝑌 . Right panel: Vertical smearing of the actual structural difference 𝑓𝑋 − 𝑌 (black line) results in an observed difference 𝑓𝑋 − 𝑌 (red line) with bias 𝛽𝑋 − 𝑌 ≈ 0.4 𝑓𝑋 − 𝑌 ( 𝑧𝑋 − 𝑌 − 𝛾 ) . The bias model in equation (18) assumes that the inversion process conserves the model integral; i.e., that the areas under the red and black curves are the same. The apparent convergence depth at the (1 – ) confidence level is shown as the red dot. ..................................................................................................................... 36 153 Figure 1-9 Regional radial correlation functions for mature ocean basins (region BC), computed from relative aspherical variations v S /v S for the exemplar models. The radial correlation lengths in the uppermost mantle of region PS are significantly larger than those of region BC, indicating that the v S /v S structure beneath the stable continents is more vertically correlated. The strong decorrelation at 660 km for S362ANI is an artifact of the Harvard inversion scheme, which allows a discontinuity in the model perturbation at this depth (Kustowski et al., 2008). ............. 40 Figure 1-10 Regional radial correlation functions for stable continent (region PS), computed from relative aspherical variations v S /v S for the exemplar models. The radial correlation lengths in the uppermost mantle of region PS are significantly larger than those of region BC (see Figure 1-9), indicating that the v S /v S structure beneath the stable continents is more vertically correlated. The relatively high radial correlation coefficient (~0.6) between PS structure in the uppermost mantle and PS structure near the core-mantle boundary is a consistent feature of the model ensemble. ........................................................................................................................... 41 Figure 1-11 Smoothing half-widths (z) of the relative aspherical variations v S /v S for S362ANI. Three estimates are shown:  BC / 2 (solid curve), half-widths of the B-splines used by Kustowski et al. (2008) to represent S362ANI upper-mantle structure (points), and a linear fit to the latter with  0 = 33 km, b = 0.08 (thick dashed line). The thin dashed lines are the ±20% uncertainty bands on the linear model. ............................................................................................................ 42 Figure 1-12 Conditional bias models for the PS-BC differential profiles, computed from equation (1-19) using the parameter values from Table 3. Dashed lines are ±20% uncertainty bands on the S362ANI bias model. .................................................................................................................... 44 154 Figure 1-13 Convergence depths of the differential profiles for A-B (left), C-B (middle), and PS- BC (right) as a function of p-value, computed from equation (17) using the bias parameters of Table 3. Dashed black line is the ensemble mean (Bayesian model average), which combines the aleatory and epistemic uncertainties. We interpret the p-value as the probability that a test depth z is greater than 𝑧𝑋 − 𝑌 . ................................................................................................................ 45 Figure 1-14 Test of the bias-correction procedure using the differential shear velocity profiles from the Faul-Jackson theoretical models (solid lines). These theoretical profiles were smoothed using the (z) values for TX2008 (dashed lines), and the bias-correction procedure was applied to recover the convergence depths (solid points). ............................................................................. 51 Figure 1-15 Comparison of shear-velocity profiles from Pacific waveform inversion studies (lines) with the lower bounds on the A-B and C-B convergence depths determined from the exemplar ensemble of this study (box). East Pacific Rise (EPR) model is by Gu et al. (2005), PA5 model is by Gaherty et al. (1996), and Old Pacific and N. Fiji Basin models are by Xu and Wiens (1997). The Pacific convergence depths from these regional models (stippled) are consistent with the tomographic lower bounds, which are labeled with their confidence levels. 58 Figure 1-16 Comparison of the convergence-depth probabilities from the tomographic ensemble (left side) with the Faul and Jackson (2005) models of upper mantle temperature (middle) and shear velocity (right). The PS average of the original Faul-Jackson model (dashed blue lines) does not converge with the oceanic averages in the upper mantle, so a linear correction was added to force convergence at 400 km (solid blue lines). Colored dots on the shear-velocity profiles show the convergence depths recovered from the Faul-Jackson models in Figure 1-14. The probable depth range for the G discontinuity, 55-70 km, is from Schmerr, (2012)’s SS precursor data. ............................................................................................................................... 59 155 Figure 2-1 Selected tomographic models at 150 km depth........................................................... 76 Figure 2-2- Regionalization in 1˚x1˚ and 2.5My age bins based on oceanic ages from (Müller et al., 2008) with marginal basins excluded...................................................................................... 78 Figure 2-3- Shear-velocity profiles and 1-sigma error bars obtained by projecting the models onto OAR1. The shear velocities increase as the square root of crustal age with apparent convergence depths greater than 200 km. ..................................................................................... 79 Figure 2-4- Differential shear-velocity profiles and 1-sigma error bars obtained by projecting the models onto OAR1’s combined A, B, and C regions. The shear velocities increase as the square root of crustal age with apparent convergence depths greater than 200 km. ................................ 81 Figure 2-5 Regional radial correlation functions for mature ocean basins (region BC), computed from relative aspherical variations vS/vS for the model ensemble. 0.78 contours are outlined in black. ............................................................................................................................................. 83 Figure 2-6- Convergence depths of the differential profiles for A-B (left), C-B (right) as a function of p-value, computed from using the method described in Part 1 and bias parameters of Table 2-2. The black line is the ensemble mean (Bayesian model average), which combines the aleatory and epistemic uncertainties. We interpret the p-value as the probability that a test depth z is greater then 𝑧𝑋 − 𝑌 . .............................................................................................................. 85 Figure 2-7 - OAR1 region ensemble and individual aspherical travel times with 1σ error bars .. 88 Figure 2-8 - Regionalization of individual ocean basins .............................................................. 89 Figure 2-9- Ensemble and individual travel times by ocean basin with 1 σ error bar .................. 91 Figure 2-10- Travel time curves by OAR age bins and ocean basin for ensemble and individual tomographic models. ..................................................................................................................... 92 156 Figure 2-11- Travel time curves by 5 million year age bins and ocean basin for ensemble and individual tomographic models..................................................................................................... 94 Figure 2-12- Relative surface area contributions for OAR and √𝟓 my age bins by ocean basin. Bins are centered at mean area-weighted age location within each bin. ...................................... 95 Figure 2-13 Method used to create the “ideal” travel time model. A strict square root of age ocean velocity model is created, then the continental velocities are derived from the tomographic model, and the resulting raw model is smoothed to create a smoothed ideal travel time model. . 96 Figure 2-14- Best fitting ensemble smoothed tomographic models (red solid line) compared to initial tomographic model (blue solid line) by ocean basin. The HSC input model pre-smoothing is shown as a red dashed line. The √𝟓 my tomographic age bins are shown as blue dots. .......... 97 Figure 2-15 95% error ellipses on the slope and intercept of the best-fitting half-space cooling model for the individual models in the ensemble. ........................................................................ 98 Figure 2-16 Initial fit of HSC travel time curve, assuming a purely linear starting HSC model, compared to the travel time curve for the SAVANI model’s combined PacIndAtl region, assuming a half-space cooling starting model which varies for each ocean basin. This results in a slightly non-linear starting HSC model, and better fits of our smoothed model to the tomographic model. The results for the remainder of the models are plotted in Appendix G. ...................... 100 Figure 2-17 Mean ensemble tomographic models (blue dots), plotted against mean best fitting ensemble half-space cooling models (red line). Many of the bumps in the tomographic models are recreated in the half-space cooling models post-smoothing, indicating they are the result of the smoothing process rather than actual physical features. ....................................................... 102 157 Figure 2-18 Aspherical V S for the SAVANI tomographic model by ocean basin (column1); half- space cooling oceanic models smoothed using κ and vertical smoothing parameters calculated from SAVANI (column 2); ......................................................................................................... 103 Figure 2-19 All seafloor vs. normal Seafloor (as defined by Korenaga and Korenaga, 2008b). The removal of the anomalous seafloor has a very small effect on the travel time profiles except at the oldest seafloor ages in the Indian Ocean, where the surface areas of the age regions are very small (Figure 2-12) so small fluctuations in travel times may have a disproportionally large effect on the average for these small area age bins. .................................................................... 104 Figure 3-1 Comparison of T-P calculation methods for carbon bearing xenoliths. Colored circles denote methods which yielded carbon phases in the correct stability field. Gray symbols denote methods which yielded carbon phases in the incorrect stability field. Solid symbols indicate diamond-bearing xenoliths, and hollow symbols indicate graphite-bearing xenoliths............... 115 Figure 3-2 Temperature vs. Pressure for the xenolith dataset for three thermobarometer combinations. Solid circles denote non-sheared xenoliths, while hollow circles denote sheared xenoliths. Even though all three methods yield carbon bearing phases in the correct diamond stability field, there is a possible temperature spread of ~100° C and pressure spread of ~8 kbar. ..................................................................................................................................................... 116 Figure 3-3 Ol Mg# vs. density for a suite of kimberlite-derived xenoliths, the data shows a roughly 2 nd degree polynomial trend (with scatter). The equation is the result of a 2 nd degree polynomial regression of the data. The coefficient of determination for the regression (R 2 ) is 0.71 and the root mean square error is 0.0015 g/cm 3 . ................................................................. 118 158 Figure 3-4 Temperature-Depth relationships for the xenolith dataset by craton. Solid circles denote unsheared samples, while hollow circles denote sheared samples. Individual cratons show distinct P-T relationships. The solid lines denote the 40 mW geotherm and the 1350°C adiabat respectively. .................................................................................................................... 119 Figure 3-5 Relationship between temperature and Ol Mg# based density for the xenolith suite with isopycnic lines for various adiabats. The simple linear regression fit line is shown in black. ..................................................................................................................................................... 121 Figure 3-6 Density vs. depth for the xenolith suite, with isopycnic lines for various adiabats. . 122 Figure 3-7 Dataset plotted with linear regression line and 95% confidence intervals on the slope of the linear regression. The 1300, 1350, and 1400° C isopycnic lines are also shown. The slope of the regression appears to be artificially shallowed by the presence of several high density outliers......................................................................................................................................... 126 Figure 3-8 Density-Depth relationship for xenolith dataset plotted with the estimated slope and its 95% confidence interval calculated by the Theil method. The green line shows the 1350 °C isopycnic curve, which is consistent with the Theil estimate within the confidence interval. ... 127 159 References Abbot, D., 1991, The case for accretion of the tectosphere by buoyant subduction: Geophysical Research Letters, v. 18, no. 4, p. 585–588. Afonso, J.C., Ranalli, G., and Fernàndez, M., 2005, Thermal expansivity and elastic properties of the lithospheric mantle: results from mineral physics of composites: Physics of the Earth and Planetary Interiors, v. 149, no. 3-4, p. 279–306, doi: 10.1016/j.pepi.2004.10.003. Arndt, N.T., Coltice, N., Helmstaedt, H., and Gregoire, M., 2009, Origin of Archean subcontinental lithospheric mantle: Some petrological constraints: Lithos, v. 109, no. 1-2, p. 61–71, doi: 10.1016/j.lithos.2008.10.019. Artemieva, I.M., 2006, Global 1 ° × 1 ° thermal model TC1 for the continental lithosphere : Implications for lithosphere secular evolution: Tectonophysics, v. 416, p. 245–277, doi: 10.1016/j.tecto.2005.11.022. Artemieva, I.M., 2009, The continental lithosphere: Reconciling thermal, seismic, and petrologic data: Lithos, v. 109, no. 1-2, p. 23–46, doi: 10.1016/j.lithos.2008.09.015. Artemieva, I.M., Thybo, H., Kaban, M.K., Survey, U.S.G., and Park, M., 2006, Deep Europe today : geophysical synthesis of the upper mantle structure and lithospheric processes over 3 . 5 Ga: European Lithosphere Dynamics, v. 32, p. 11–41. Audet, P., and Bürgmann, R., 2011, Dominant role of tectonic inheritance in supercontinent cycles: Nature Geoscience, v. 4, no. 3, p. 184–187, doi: 10.1038/ngeo1080. Auer, L., Boschi, L., Becker, T.W., and Giardini, D., 2014, Savani : A variable resolution whole- mantle model of anisotropic shear velocity variations based: Journal of Geophysical Research: Solid Earth, v. 119, p. 1–29, doi: 10.1002/2013JB010773.Received. Backus, G., and Gilbert, F., 1968, The Resolving Power of Gross Earth Data: Geophysical Journal International, v. 16, no. 2, p. 169–205, doi: 10.1111/j.1365-246X.1968.tb00216.x. Bagley, B., and Revenaugh, J., 2008, Upper mantle seismic shear discontinuities of the Pacific: v. 113, no. December, p. 1–9, doi: 10.1029/2008JB005692. Baptiste, V., Tommasi, A., and Demouchy, S., 2012, Deformation and hydration of the lithospheric mantle beneath the Kaapvaal craton, South Africa: Lithos, v. 149, p. 31–50, doi: 10.1016/j.lithos.2012.05.001. Barrell, J., 1914, The strength of the Earth’s crust: The Journal of Geology, v. 22, no. 4, p. 289– 314. Becker, T.W., and Boschi, L., 2002, A comparison of tomographic and geodynamic mantle models: 2Geochemistry, Geophysics, Geosystems, v. 3. Becker, T.W., Kustowski, B., and Ekström, G., 2008, Radial seismic anisotropy as a constraint for upper mantle rheology: v. 267, p. 213–227, doi: 10.1016/j.epsl.2007.11.038. Beghein, C., Yuan, K., Schmerr, N., and Xing, Z., 2009, Changes in seismic anisotropy shed light on the nature of the Gutenberg discontinuity: Science, v. 343, no. 2014, p. 1237–1249. Bercovici, D., and Karato, S., 2003, Whole-mantle convection and the transition-zone water filter: Nature, v. 425, p. 39–44. Bommer, J.J., and Scherbaum, F., 2008, The Use and Misuse of Logic Trees in Probabalistic Seismic Hazard Analysis: Earthquake Spectra, v. 24, no. 4, p. 997–1009, doi: 10.1193/1.2977755. 160 Bostock, M.G., Eaton, D.W., and Snyder, D.B., 2010, Teleseismic studies of the Canadian landmass: Lithoprobe and its legacy: Canadian Journal of Earth Sciences, v. 47.4, p. 445. Brey, G.P., Kohler, T., and Nickel, K.G., 1990, Geothermobarometry in four-phase lherzolites I . Experimental results from 10 to 60 kb: Journal of Petrology, v. 31, p. 1313–1352. Budnitz, R.J., Apostolakis, G., Boore, D.M., Cluff, L.S., Coppersmith, K.J., Cornell, C.A., and Morris, P.A., 1997, Recommendations for Probabilistic Seismic Hazard Analysis : Guidance on Uncertainty and Use of Experts: Power, v. 1, no. 6, p. 998–1006. Burgos, G., Montagner, J., Beucler, E., Capdeville, Y., Mocquet, A., and Drilleau, M., 2014, Journal of Geophysical Research Oceanic lithosphere-asthenosphere boundary from surface wave dispersion data: Journal of Geophysical Research: Solid Earth, v. 119, p. 1079–1093, doi: 10.1002/2013JB010528.Abstract. Burov, E.B., and Diament, M., 1995, The effective elastic thickness (Te) of continental lithosphere: What does it really mean? Journal of Geophysical Research, v. 100, no. B3, p. 3905–3927. Cammarano, F., and Romanowicz, B.A., 2007, Insights into the nature of the transition zone from physically constrained inversion of long-period seismic data.: Proceedings of the National Academy of Sciences of the United States of America, v. 104, no. 22, p. 9139– 9144. Canil, D., and Lee, C.A., 2009, Were deep cratonic mantle roots hydrated in Archean oceans ? Geology, v. 37, no. 7, p. 667–670, doi: 10.1130/G25610A.1. Cao, Q., van der Hilst, R.D., de Hoop, M. V, and Shim, S.-H., 2011, Seismic imaging of transition zone discontinuities suggests hot mantle west of Hawaii.: Science (New York, N.Y.), v. 332, no. 6033, p. 1068–1071, doi: 10.1126/science.1202731. Carlson, R.W., James, D.E., Boyd, F.R., Schutt, D., and Bell, D.R., 2004, Xenolith constraints on seismic velocities in the upper mantle beneath southern Africa: Geochemistry, Geophysics, Geosystems, v. 5, no. 1, p. n/a–n/a, doi: 10.1029/2003GC000551. Carlson, R.W., Pearson, D.G., and James, D.E., 2005, Physical, chemical, and chronological characteristics of continental mantle: Reviews of Geophysics, v. 43, no. RG1001, p. 1–24, doi: 10.1029/2004RG000156.1.INTRODUCTION. Conrad, C.P., and Lithgow-Bertelloni, C., 2006, Influence of continental roots and asthenosphere on plate-mantle coupling: Geophysical Research Letters, v. 33, no. 5, p. L05312, doi: 10.1029/2005GL025621. Courtier, A.M., and Revenaugh, J., 2007, Deep upper-mantle melting beneath the Tasman and Coral Seas detected with multiple ScS reverberations: Earth and Planetary Science Letters, v. 259, no. 1-2, p. 66–76, doi: 10.1016/j.epsl.2007.04.027. Cox, D.R., and Hinkley, D. V., 1974, Theoretical Statistics: Chapman and Hall, London. Crosby, A.G., McKenzie, D., and Sclater, J.G., 2006, The relationship between depth, age and gravity in the oceans: Geophysical Journal International, v. 166, no. 2, p. 553–573, doi: 10.1111/j.1365-246X.2006.03015.x. Dalton, C. a., Ekström, G., and Dziewonski, A.M., 2009, Global seismological shear velocity and attenuation: A comparison with experimental observations: Earth and Planetary Science Letters, v. 284, no. 1-2, p. 65–75, doi: 10.1016/j.epsl.2009.04.009. Dalton, C.A., and Faul, U.H., 2010, The oceanic and cratonic upper mantle: Clues from joint interpretation of global velocity and attenuation models: Lithos, v. 120, no. 1-2, p. 160–172, doi: 10.1016/j.lithos.2010.08.020. 161 Dalton, C.A., Langmuir, C.H., and Gale, A., 2014, Geophysical and geochemical evidence for deep temperature variations beneath mid-ocean ridges: Science, v. 344, no. April, p. 80–83. Debayle, E., and Ricard, Y., 2012, A global shear velocity model of the upper mantle from fundamental and higher Rayleigh mode measurements: Journal of Geophysical Research, v. 117, no. B10, p. B10308, doi: 10.1029/2012JB009288. Deschamps, F., Trampert, J., and Snieder, R., 2002, Anomalies of temperature and iron in the uppermost mantle inferred from gravity data and tomographic models: Physics of the Earth and Planetary Interiors, v. 129, p. 245–264. Doin, M.P., and Fleitout, L., 1996, Thermal evolution of the oceanic lithosphere: an alternative view: Earth and Planetary Science Letters, v. 142, no. 1-2, p. 121–136, doi: 10.1016/0012- 821X(96)00082-9. Doin, M., Fleitout, L., and McKenzie, D., 1996, crust: Journal of Geophysical Research, v. 101, no. B7, p. 16119–16135. Dziewonski, A.M., 1970, On Regional Differences in Dispersion of Mantle Rayleigh Waves: Geophys. J R. astr. Soc., v. 22, p. 289–325. Dziewoński, A.M., and Anderson, D.L., 1981, Preliminary reference Earth model: Physics of the Earth and Planetary Interiors, v. 25, p. 297–356. Eaton, D.W., Darbyshire, F.A., Evans, R.L., Grütter, H.S., Jones, A.G., and Yuan, X., 2009, The elusive lithosphere-asthenosphere boundary (LAB) beneath cratons: Lithos, v. 109, no. 1-2, p. 1–22, doi: 10.1016/j.lithos.2008.05.009. Ekström, G., and Dziewoński, A.M., 1998, The unique anisotropy of the Pacific upper mantle: Nature, v. 394, no. July, p. 168–173. Faul, U.H., and Jackson, I., 2005, The seismological signature of temperature and grain size variations in the upper mantle: Earth and Planetary Science Letters, v. 234, p. 119–134. Ferreira, a. M.G., Woodhouse, J.H., Visser, K., and Trampert, J., 2010, On the robustness of global radially anisotropic surface wave tomography: Journal of Geophysical Research: Solid Earth, v. 115, no. 4, p. 1–17, doi: 10.1029/2009JB006716. Field, E.H., Dawson, T.E., Felzer, K.R., Frankel, A.D., Gupta, V., Jordan, T.H., Parsons, T., Petersen, M.D., Stein, R.S., Weldon II, R.J., and Wills, C.J., 2009, Uniform California Earthquake Rupture Forecast: Bulletin of the Seismological Society of America, v. 99, p. 2053–2107, doi: 10.1785/0120080049. Fischer, K.M., Ford, H.A., Abt, D.L., and Rychert, C.A., 2010, The Lithosphere-Asthenosphere Boundary: Annual Review of Earth and Planetary Sciences, v. 38, no. 1, p. 551–575, doi: 10.1146/annurev-earth-040809-152438. Fishwick, S., Kennett, B.L.N., and Reading, A.M., 2005, Contrasts in lithospheric structure within the Australian craton—insights from surface wave tomography: Earth and Planetary Science Letters, v. 231, no. 3-4, p. 163–176, doi: 10.1016/j.epsl.2005.01.009. Forte, A.M., and Perry, H.K.C., 2000, Geodynamic evidence for a chemically depleted continental tectosphere.: Science (New York, N.Y.), v. 290, no. 5498, p. 1940–4. Fouch, M.J., James, D.E., Vandecar, J.C., and van der Lee, S., 2004, Mantle seismic structure beneath the Kaapvaal and Zimbabwe Cratons: South African Journal of Geology, v. 107, p. 33–44. French, S., Lekic, V., and Romanowicz, B., 2013, Waveform tomography reveals channeled flow at the base of the oceanic asthenosphere.: Science (New York, N.Y.), v. 342, no. 6155, p. 227–30, doi: 10.1126/science.1241514. 162 Gaherty, J.B., and Jordan, T.H., 1995, Lehmann discontinuity as the base of an anisotropic layer beneath continents: Science, v. 268, no. 5216, p. 1468–1471. Gaherty, J.B., Jordan, T.H., and Gee, L.S., 1996, Seismic structure of the upper mantle in a central Pacific corridor: Journal of Geophysical Research, v. 101, no. B10, p. 22291–22309. Gaherty, J.B., Kato, M., and Jordan, T.H., 1999, Seismological structure of the upper mantle: a regional comparison of seismic layering: Physics of the Earth and Planetary Interiors, v. 110, p. 21–41. Gaul, O.F., Griffin, W.L., O’Reilly, S.Y., and Pearson, N.J., 2000, Mapping olivine composition in the lithospheric mantle: Earth and Planetary Science Letters, v. 182, p. 223–235. Goodman, S.N., 1999, Toward Evidence-Based Medical Statistics. 1: The: Ann Intern Med., v. 130, p. 995–1004. Grand, S.P., 1994, Mantle shear structure beneath the Americas and surrounding oceans: Journal of Geophysical Research, v. 99, no. B6, p. 11591–11621. Grand, S.P., 2002, Mantle shear-wave tomography and the fate of subducted slabs: Philosophical Transactions: Mathematical, Physical and Engineering Sciences, v. 360, no. 1800, p. 2475– 2491. Grand, S.P., and Helmberger, D. V, 1984, Upper mantle shear structure of North America: Geophys. J R. astr. Soc., v. 76, p. 399–438. Grand, S.P., van der Hilst, R.D., and Widyantoro, S., 1997, Global seismic tomography: A snapshop of convection in the Earth: GSA Today, v. 7, no. 4, p. 1–7. Griffin, W.L., Doyle, B.J., Ryan, C.G., Pearson, N.J., O’Reilly, S.Y., Davies, R.M., Kivi, K., Achterbergh, E.V.A.N., and Natapov, L.M., 1999, Layered Mantle Lithosphere in the Lac de Gras Area , Slave Craton : Composition , Structure and Origin: v. 40, no. 5, p. 705–727. Griffin, W.., O’Reilly, S.., Abe, N., Aulbach, S., Davies, R.., Pearson, N.., Doyle, B.., and Kivi, K., 2003, The origin and evolution of Archean lithospheric mantle: Precambrian Research, v. 127, no. 1-3, p. 19–41, doi: 10.1016/S0301-9268(03)00180-3. Griffin, W.L., O’Reilly, S.Y., Afonso, J.C., and Begg, G.C., 2009, The composition and evolution of lithospheric mantle: A re-evaluation and its tectonic implications: Journal of Petrology, v. 50, no. 7, p. 1185–1204, doi: 10.1093/petrology/egn033. Gu, Y.J., 2005, Upper mantle structure beneath the eastern Pacific Ocean ridges: Journal of Geophysical Research, v. 110, no. B6, p. 1–18, doi: 10.1029/2004JB003381. Gu, Y.J., Dziewoński, A.M., and Ekström, G., 2003, Simultaneous inversion for mantle shear velocity and topography of transition zone discontinuities: Geophysical Journal International, v. 154, no. 2, p. 559–583. Gu, Y.J., Dziewoński, A.M., Su, W.J., and Ekström, G., 2001, Models of the mantle shear velocity and discontinuities in the pattern of lateral heterogeneities: Journal of Geophysical Research, v. 106, no. B6, p. 11169–11199. Gung, Y., Panning, M.P., and Romanowicz, B.A., 2003, Global anisotropy and the thickness of continents: Nature, v. 422, no. April, p. 707–711, doi: 10.1038/nature01557.1. Hirth, G., and Kohlstedt, D.L., 1996, melt extraction and the evolution of the lithosphere: v. 144, p. 93–108. Hoeting, J.A., Madigan, D., Raftery, A.E., and Volinsky, C.T., 1999, Bayesian Model Averaging : A Tutorial: Statistical Science, v. 14, no. 4, p. 382–401. Höink, T., Lenardic, A., and Richards, M., 2012, Depth-dependent viscosity and mantle stress amplification: Implications for the role of the asthenosphere in maintaining plate tectonics: 163 Geophysical Journal International, v. 191, no. 1, p. 30–41, doi: 10.1111/j.1365- 246X.2012.05621.x. Hollander, M., and Wolfe, D., 1973, Nonparametric Statistical Methods: John Wiley & Sons, Inc., New York. Houseman, G.A., Mckenzie, D.P., and Molnar, P., 1981, Convective instability of a thickened boundary layer and its relevance for the thermal evolution of continental convergent belts: Journal of Geophysical Research, v. 86, p. 6115–6132. Houser, C., Masters, G., Shearer, P.M., and Laske, G., 2008, Shear and compressional velocity models of the mantle from cluster analysis of long-period waveforms: Geophysical Journal International, v. 174, no. 1, p. 195–212, doi: 10.1111/j.1365-246X.2008.03763.x. Huang, J., 2003, Controls on sublithospheric small-scale convection: Journal of Geophysical Research, v. 108, no. B8, p. 1–13, doi: 10.1029/2003JB002456. van Hunen, J., 2003, The effect of shearing on the onset and vigor of small-scale convection in a Newtonian rheology: Geophysical Research Letters, v. 30, no. 19, p. 1–4, doi: 10.1029/2003GL018101. Ishii, M., and 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, v. 145, no. 1, p. 77–96. James, D.E., Boyd, F.R., Schutt, D., Bell, D.R., and Carlson, R.W., 2004a, Xenolith constraints on seismic velocities in the upper mantle beneath southern Africa: Geochemistry, Geophysics, Geosystems, v. 5, no. 1, p. n/a–n/a, doi: 10.1029/2003GC000551. James, D.E., Boyd, F.R., Schutt, D.L., Bell, D.R., and Carlson, R.W., 2004b, Xenolith constraints on seismic velocities in the upper mantle beneath southern Africa: Geochemistry, Geophysics, Geosystems, v. 5, no. 1, p. n/a–n/a, doi: 10.1029/2003GC000551. James, D.E., Fouch, M.J., Vandecar, J.C., and van der Lee, S., 2001, Tectospheric structure beneath southern Africa: Geophysical Research Letters, v. 28, no. 13, p. 2485–2488. Jasbinsek, J.J., Dueker, K.G., and Hansen, S.M., 2010, Characterizing the 410 km discontinuity low-velocity layer beneath the LA RISTRA array in the North American Southwest: Geochemistry, Geophysics, Geosystems, v. 11, no. 3, p. n/a–n/a, doi: 10.1029/2009GC002836. Jordan, T.H., 1978a, Composition and development of the continental tectosphere: Nature, v. 274, p. 544–548. Jordan, T.H., 1981a, Continents as a chemical boundary layer: Philosophical transactions. Series A, Mathematical, physical, and engineering sciences, v. 301, no. 1461, p. 359–373. Jordan, T.H., 1981b, Global tectonic regionalization for seismological data analysis: Bulletin of the Seismological Society of America, v. 71, no. 4, p. 1131–1141. Jordan, T.H., 2010, Implications of the Strong Isopycnic Hypothesis for Cratonic Evolution, in Out of Africa,. Jordan, T.H., 1978b, Implications of the Strong Isopycnic Hypothesis for the Evolution of Continental Cratons: , no. 1, p. 90089. Jordan, T.H., 1979, Mineralogies, densities and seismic velocities of garnet lherzolites and their geophysical implications (F. R. Boyd & H. O. A. Meyer, Eds.): The Mantle Sample: Inclusion in Kimberlites and Other Volcanics, v. 16, p. 1–14, doi: doi:10.1029/SP016p0001. Jordan, T.H., 1988, Structure and formation of the continental tectosphere: Journal of Petrology, , no. Special Lithosphere Issue, p. 11–37. 164 Jordan, T.H., 1975, The continental tectosphere: Reviews of geophysics and space physics, v. 13, no. 3, p. 1–12. Jordan, T.H., and Gaherty, J.B., 1995, Stochastic modeling of small-scale, anisotropic structures in the continental upper mantle:. Jordan, T.H., and Paulson, E.M., 2013, Convergence depths of tectonic regions from an ensemble of global tomographic models: Journal of Geophysical Research: Solid Earth, v. 118, no. August, p. 4196–4225, doi: 10.1002/jgrb.50263. Jordan, T.H., Puster, P., Glatzmaier, G.A., and Tackley, P.J., 1993, Comparisons between seismic Earth structures and mantle flow models based on radial correlation functions: Science, v. 261, no. 5127, p. 1427–1431. Kaban, M.K., Schwintzer, P., Artemieva, I.M., and Mooney, W.D., 2003, Density of the continental roots : compositional and thermal contributions: Earth and Planetary Science Letters, v. 209, no. 1-2, p. 53–69, doi: 10.1016/S0012-821X(03)00072-4. Kaban, M.K., Schwintzer, P., and Reigber, C., 2004, A new isostatic model of the lithosphere and gravity field: Journal of Geodesy, v. 78, no. 6, p. 368–385, doi: 10.1007/s00190-004- 0401-6. Karato, S., 2012, On the origin of the asthenosphere: Earth and Planetary Science Letters, v. 321- 322, p. 95–103, doi: 10.1016/j.epsl.2012.01.001. Karato, S., 2011, Water distribution across the mantle transition zone and its implications for global material circulation: Earth and Planetary Science Letters, v. 301, no. 3-4, p. 413–423, doi: 10.1016/j.epsl.2010.11.038. Katzman, R., Zhao, L., and Jordan, T.H., 1998, High-resolution, two-dimensional vertical tomography of the central Pacific mantle using ScS reverberations and frequency-dependent travel times: Journal of Geophysical Research, v. 103, no. B8, p. 17933–17971. Kawakatsu, H., Kumar, P., Takei, Y., Shinohara, M., Kanazawa, T., Araki, E., and Suyehiro, K., 2009, Seismic Evidence for Sharp Lithosphere-Asthenosphere Boundaries of Oceanic Plates: Science, v. 324, no. 5926, p. 499–502. Kearey, P., Klepeis, K.A., and Vine, F.J., 2009, Global Tectonics: Wiley-Blackwell. Kelly, R.K., Kelemen, P.B., and Jull, M., 2003, Buoyancy of the continental upper mantle: Geochemistry, Geophysics, Geosystems, v. 4, no. 2, p. n/a–n/a, doi: 10.1029/2002GC000399. Kennedy, C.S., and Kennedy, G.C., 1976, The equilibrium boundary between graphite and diamond: Journal of Geophysical Research, v. 81, no. 14, p. 2467–2470, doi: doi:10.1029/JB081i014p02467. Khan, A., Boschi, L., and Connolly, J.A.D., 2011, Mapping the Earth’s thermochemical and anisotropic structure using global surface wave data: Journal of Geophysical Research, v. 116, no. B1, p. B01301, doi: 10.1029/2010JB007828. Korenaga, J., 2015, Seafloor topography and the thermal budget of Earth: Geological Society of America Special Papers, v. 514, p. 167–185. Korenaga, J., and Jordan, T.H., 2004, Physics of multiscale convection in Earth’s mantle: Evolution of sublithospheric convection: Journal of Geophysical Research, v. 109, no. B1, p. 1–11, doi: 10.1029/2003JB002464. Korenaga, J., and Jordan, T.H., 2003, Physics of multiscale convection in Earth’s mantle: Onset of sublithospheric convection: Journal of Geophysical Research, v. 108, no. B7, p. 2333, doi: 10.1029/2003JB002464. Korenaga, T., and Korenaga, J., 2008a, Subsidence of normal oceanic lithosphere, apparent 165 thermal expansivity, and seafloor flattening: Earth and Planetary Science Letters, v. 268, no. 1-2, p. 41–51, doi: 10.1016/j.epsl.2007.12.022. Korenaga, T., and Korenaga, J., 2008b, Subsidence of normal oceanic lithosphere, apparent thermal expansivity, and seafloor flattening.pdf:. Kustowski, B., Ekström, G., and Dziewoński, A.M., 2008, Anisotropic shear-wave velocity structure of the Earth’s mantle: A global model: Journal of Geophysical Research, v. 113, no. B6, p. B06306, doi: 10.1029/2007JB005169. Landuyt, W., and Ierley, G., 2012, Linear stability analysis of the onset of sublithospheric convection: Geophysical Journal International, v. 189, no. 1, p. 19–28, doi: 10.1111/j.1365- 246X.2011.05341.x. Leahy, G.M., and Bercovici, D., 2007, On the dynamics of a hydrous melt layer above the transition zone: Journal of Geophysical Research, v. 112, no. B7, p. B07401, doi: 10.1029/2006JB004631. Lee, C.-T.A., 2003, Compositional variation of density and seismic velocities in natural peridotites at STP conditions: Implications for seismic imaging of compositional heterogeneities in the upper mantle: Journal of Geophysical Research, v. 108, no. B9, p. 2441, doi: 10.1029/2003JB002413. Lee, C.-T.A., Luffi, P., and Chin, E.J., 2011, Building and Destroying Continental Mantle: Annual Review of Earth and Planetary Sciences, v. 39, no. 1, p. 59–90, doi: 10.1146/annurev-earth-040610-133505. Van der Lee, S., and Frederiksen, A., 2005, Surface wave tomography applied to the North American upper mantle (A. Levander & G. Nolet, Eds.):. Lekić, V., Panning, M.P., and Romanowicz, B.A., 2010, A simple method for improving crustal corrections in waveform tomography: Geophysical Journal International, v. 182, p. 265– 278, doi: 10.1111/j.1365-246X.2010.04602.x. Lekić, V., and Romanowicz, B.A., 2011, Inferring upper-mantle structure by full waveform tomography with the spectral element method: Geophysical Journal International, v. 185, p. 799–831, doi: 10.1111/j.1365-246X.2011.04969.x. Lerner-Lam, A.L., and Jordan, T.H., 1987, How thick are the continents? Journal of Geophysical Research, v. 92, no. B13, p. 7–14. MacGregor, I.D., 1974, The system MgO-Al2O3-SiO2: Solubility of Al2O3 in enstatite for spinel and garnet peridotite compositions: American Mineralogist, v. 59, p. 110–119. Masters, G., Johnson, S., Laske, G., and Bolton, H., 1996, A shear-velocity model of the mantle: v. 354, no. 1711, p. 1385–1411. Masters, G., Laske, G., Bolton, H., and Dziewonski, A.M., 2000, The relative behavior of shear velocity, bulk sound speed, and compressional velocity in the mantle: Implications for chemical and thermal structure, in Karato, S.-I., Forte, A.M., Liebermann, R.C., Masters, G., and Stixrude, L. eds., Mineral Physics and Tomography from the Atomic to the Global Scale, AGU, Washington, D. C., p. 63–87. McDonough, W.F., 2013, Garnet lherzolite data compilation provided at 2013 CIDER conference:. Megnin, C., and Romanowicz, B.A., 2000, The three-dimensional shear velocity structure of the mantle from the inversion of body , surface and higher-mode waveforms: Geophysical Journal International, v. 143, p. 709–728. Montagner, J.P., 2002, Upper mantle low anisotropy channels below the Pacific Plate: Earth and Planetary Science Letters, v. 202, no. 2, p. 263–274, doi: 10.1016/S0012-821X(02)00791-4. 166 Montagner, J.-P., and Tanimoto, T., 1991, Global Upper Mantle Tomography of Seismic Velocities and Anisotropies: Journal of Geophysical Research, v. 96, no. B12, p. 20337– 20351. Montelli, R., Nolet, G., Dahlen, F. a., and Masters, G., 2006, A catalogue of deep mantle plumes: New results from finite-frequency tomography: Geochemistry, Geophysics, Geosystems, v. 7, no. 11, p. n/a–n/a, doi: 10.1029/2006GC001248. Müller, R.D., Sdrolias, M., Gaina, C., and Roest, W.R., 2008, Age, spreading rates, and spreading asymmetry of the world’s ocean crust: Geochemistry, Geophysics, Geosystems, v. 9, no. 4, doi: 10.1029/2007GC001743. Nagihara, S., Lister, C.R.B., and Sclater, J.G., 1996, Reheating of old oceanic lithosphere : Deductions from observations: Earth and Planetary Science Letters, v. 139, p. 91–104. Naif, S., Key, K., Constable, S., and Evans, R.L., 2013, Melt-rich channel observed at the lithosphere-asthenosphere boundary.: Nature, v. 495, no. 7441, p. 356–9, doi: 10.1038/nature11939. Nettles, M., and Dziewoński, A.M., 2008, Radially anisotropic shear velocity structure of the upper mantle globally and beneath North America: Journal of Geophysical Research, v. 113, no. B2, p. B02303, doi: 10.1029/2006JB004819. Nickel, K.G., and Green, D.H., 1985, Empirical geothermobarometry for garnet peridotites and implications for the nature of the lithosphere, kimberlites and diamonds: Earth and Planetary Science Letters, v. 73, no. 1, p. 158–170, doi: 10.1016/0012-821X(85)90043-3. Nimis, P., and Grütter, H.S., 2010, Internally consistent geothermometers for garnet peridotites and pyroxenites: Contributions to Mineralogy and Petrology, v. 159, no. 3, p. 411–427, doi: 10.1007/s00410-009-0464-8. Nimis, P., and Taylor, W.R., 2000, Single clinopyroxene thermobarometry for garnet peridotites. Part I. Calibration and testing of a Cr-in-Cpx barometer and an enstatite-in-Cpx thermometer: Contributions to Mineralogy and Petrology, v. 139, no. 5, p. 541–554, doi: 10.1007/s004100000156. Nyblade, A.A., and Pollack, H.N., 1993, A Global Analysis of Heat Flow From Precambrian Terrains ’ Implications for the Thermal Structure of Archcan and Proterozoic Limosphere: Journal of Geophysical Research: Solid Earth, v. 98, no. B7, p. 12207–12218. Okal, E.A., 1977, The effect of intrinsic oceanic upper-mantle heterogeneity on regionalization of long-period Rayleigh-wave phase velocities: , p. 357–370. Olugboji, T.M., Park, J., Karato, S., and Shinohara, M., 2016, Geochemistry, Geophysics, Geosystems: Geochem Geophys Geosyst, v. 17, p. 1265–1282, doi: 10.1002/2015GC006214.Received. Panning, M.P., Lekić, V., and Romanowicz, B.A., 2010, Importance of crustal corrections in the development of a new global model of radial anisotropy: Journal of Geophysical Research, v. 115, no. B12, p. B12325, doi: 10.1029/2010JB007520. Panning, M.P., and Romanowicz, B.A., 2006, A three-dimensional radially anisotropic model of shear velocity in the whole mantle: Geophysical Journal International, v. 167, no. 1, p. 361– 379, doi: 10.1111/j.1365-246X.2006.03100.x. Pari, G., and Peltier, W.R., 1996, The free-air gravity constraint on subcontinental mantle dynamics Abstract . An outstanding geophysical issue concerns that have been tomographically imaged beneath upper mantle and transition these to geophysically more realistic chemical models of sub .: Journal of Geophysical Research, v. 101, no. B12, p. 28105–28132. 167 Parker, R.L., and Oldenburg, D.W., 1973, Thermal model of ocean ridges: Nature Physical Science, v. 242, p. 137–139. Parman, S.W., Grove, T.L., Dann, J.C., de Wit, M.J., Wit, M.J.D.E., and Wit, M.J.D.E., 2004, A subduction origin for komatiites and cratonic lithospheric mantle: South African Journal of Geology, v. 107, no. 1-2, p. 107–118. Parsons, B., and McKenzie, D., 1978, Mantle convection and the thermal structure for plates: Journal of Geophysical Research, v. 83, no. B9, p. 4485–4496. Parsons, B., and Sclater, J.G., 1977, An analysis of the variation of ocean floor bathymetry and heat flow with age: Journal of Geophysical Research, v. 82, no. 6, p. 803–827. Pearson, D.., 1999, The age of continental roots: Lithos, v. 48, no. 1-4, p. 171–194, doi: 10.1016/S0024-4937(99)00026-2. Pearson, D.G., Canil, D., and Shirey, S.B., 2003, Mantle Samples Included in Volcanic Rocks : Xenoliths and Diamonds: Treatise on Geochemistry, v. 2. Pedersen, H. a., Bruneton, M., and Maupin, V., 2006, Lithospheric and sublithospheric anisotropy beneath the Baltic shield from surface-wave array analysis: Earth and Planetary Science Letters, v. 244, no. 3-4, p. 590–605, doi: 10.1016/j.epsl.2006.02.009. Pollack, H.N., and Chapman, D.S., 1977, Mantle heat flow: Earth and Planetary Science Letters, v. 34, p. 174–184. Poudjom, Y.H., Reilly, S.Y.O., Gri, W.L., and Morgan, P., 2001a, The density structure of subcontinental lithosphere through time: Earth and Planetary Science Letters, v. 184, no. 3- 4, p. 605–621. Poudjom, Y.H., Reilly, S.Y.O., Gri, W.L., Morgan, P., O’Reilly, S.Y., Griffin, W.L., and Morgan, P., 2001b, The density structure of subcontinental lithosphere through time: Earth and Planetary Science Letters, v. 184, no. 3-4, p. 605–621. Priestley, K., McKenzie, D., and Debayle, E., 2006, The state of the upper mantle beneath southern Africa: Tectonophysics, v. 416, no. 1-4, p. 101–112, doi: 10.1016/j.tecto.2005.11.024. Priestly et al 2006 The state of the upper mantle beneath southern Africa.pdf, 2006,. Puster, P., and Jordan, T.H., 1997, How stratified is mantle convection? Journal of Geophysical Research, v. 102, no. B4, p. 7625–7646. Puster, P., and Jordan, T.H., 1994, Stochastic analysis of mantle convection experiments using two-point correlation functions: Geophysical Research Letters, v. 21, no. 4, p. 305–308. Revenaugh, J., and Jordan, T.H., 1991, Mantle Layering From ScS Reverberations 3 . The Upper Mantle concerns ML2 uses the refe { ence frame: Journal of Geophysical Research, v. 96, no. B12, p. 19781–19810. Revenaugh, J., and Sipkin, S.A., 1994, Seismic evidence for sillicate melt atop the 410-km mantle discontinuity: Nature, v. 369, p. 474–476. Ringwood, A.E., 1975, Composition and Petrology of the Earth’s Mantle: McGraw-Hill, New York. Ritsema, J., Deuss, A., van Heijst, H.-J., and Woodhouse, J.H., 2010, 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, v. 184, no. 3, p. 1223–1236, doi: 10.1111/j.1365-246X.2010.04884.x. Ritsema, J., and Van Heijst, H.J., 2000, Seismic imaging of structural heterogeneity in Earth’s mantle: evidence for large-scale mantle flow.: Science progress, v. 83 ( Pt 3), p. 243–259. Ritsema, J., van Heijst, H.-J., and Woodhouse, J.H., 2004, Global transition zone tomography: 168 Journal of Geophysical Research, v. 109, no. B2, p. B02302, doi: 10.1029/2003JB002610. Ritzwoller, M.H., Shapiro, N.M., and Zhong, S.-J., 2004, Cooling history of the Pacific lithosphere: Earth and Planetary Science Letters, v. 226, no. 1-2, p. 69–84, doi: 10.1016/j.epsl.2004.07.032. Rychert, C. a, and Shearer, P.M., 2009, A global view of the lithosphere-asthenosphere boundary.: Science (New York, N.Y.), v. 324, no. 5926, p. 495–498, doi: 10.1126/science.1169754. Rychert, C. a., and Shearer, P.M., 2011, Imaging the lithosphere-asthenosphere boundary beneath the Pacific using SS waveform modeling: Journal of Geophysical Research, v. 116, no. B7, p. B07307, doi: 10.1029/2010JB008070. Saltzer, R.L., Gaherty, J.B., and Jordan, T.H., 2000, How are vertical shear wave splitting measurements affected by variations in the orientation of azimuthal anisotropy with depth ? Geophysical Journal International, v. 141, p. 374–390. Schaeffer, A.J., and Lebedev, S., 2013, Global shear speed structure of the upper mantle and transition zone: Geophysical Journal International, v. 194, no. 1, p. 417–449, doi: 10.1093/gji/ggt095. Schmerr, N.C., 2012, The Gutenberg discontinuity: melt at the lithosphere-asthenosphere boundary.: Science (New York, N.Y.), v. 335, no. 6075, p. 1480–3, doi: 10.1126/science.1215433. Schmitz, M.D., Bowring, S.A., de Wit, M.J., and Gartz, V., 2004, Subduction and terrane collision stabilize the western Kaapvaal craton tectosphere 2.9 billion years ago: Earth and Planetary Science Letters, v. 222, no. 2, p. 363–376, doi: 10.1016/j.epsl.2004.03.036. Schoene, B., de Wit, M.J., and Bowring, S. a., 2008, Mesoarchean assembly and stabilization of the eastern Kaapvaal craton: A structural-thermochronological perspective: Tectonics, v. 27, no. 5, p. 1–27, doi: 10.1029/2008TC002267. Schutt, D.L., and Lesher, C.E., 2010, Compositional trends among Kaapvaal Craton garnet peridotite xenoliths and their effects on seismic velocity and density: Earth and Planetary Science Letters, v. 300, no. 3-4, p. 367–373, doi: 10.1016/j.epsl.2010.10.018. Schutt, D.L., and Lesher, C.E., 2006, Effects of melt depletion on the density and seismic velocity of garnet and spinel lherzolite: Journal of Geophysical Research, v. 111, no. B5, p. B05401, doi: 10.1029/2003JB002950. Sclater, J.G., and Francheteau, J., 1970, The implications of terrestrial heat flow observations on current tectonic and geochemical models of the crust and upper mantle of the Earth: Geophysical Journal International, v. 20, no. 5, p. 509–542. Sclater, J.G., Parsons, B., and Jaupart, C., 1981, Oceans and continents: Similarities and differences in the mechanisms of heat loss: Journal of Geophysical Research, v. 86, no. B12, p. 11535, doi: 10.1029/JB086iB12p11535. Shapiro, S.S., Hager, B.H., and Jordan, T.H., 1999a, Stability and dynamics of the continental tectosphere: , no. February, p. 115–133. Shapiro, S.S., Hager, B.H., and Jordan, T.H., 1999b, The continental tectosphere and Earth’s long-wavelength gravity field: Lithos, v. 48, no. 1-4, p. 135–152, doi: 10.1016/S0024- 4937(99)00027-4. Silver, P.G., Gao, S.S., and Liu, K.H., 2001, Mantle deformation beneath southern Africa: v. 28, no. 13, p. 2493–2496. Simmons, N. a., Forte, A.M., Boschi, L., and Grand, S.P., 2010, GyPSuM: A joint tomographic model of mantle density and seismic wave speeds: Journal of Geophysical Research, v. 115, 169 no. B12, p. B12310, doi: 10.1029/2010JB007631. Simmons, N.A., Forte, A.M., and Grand, S.P., 2009, Joint seismic, geodynamic and mineral physical constraints on three-dimensional mantle heterogeneity: implications for the relative importance of thermal versus compositional heterogeneity: Geophysical Journal International, v. 177, no. 3, p. 1284–1304, doi: 10.1111/j.1365-246X.2009.04133.x. Simmons, N. a., Forte, A.M., and Grand, S.P., 2007, Thermochemical structure and dynamics of the African superplume: Geophysical Research Letters, v. 34, no. 2, p. L02301, doi: 10.1029/2006GL028009. Song, T.A., Helmberger, D. V, and Grand, S.P., 2004, Low-velocity zone atop the 410-km seismic discontinuity in the northwestern U ...:. Song, T.-R.A., and Kawakatsu, H., 2013, Subduction of oceanic asthenosphere: A critical appraisal in central Alaska: Earth and Planetary Science Letters, v. 367, p. 82–94, doi: 10.1016/j.epsl.2013.02.010. Song, T.-R.A., and Kawakatsu, H., 2012, Subduction of oceanic asthenosphere: Evidence from sub-slab seismic anisotropy: Geophysical Research Letters, v. 39, no. 17, doi: 10.1029/2012GL052639. Stein, C.A., and Stein, S., 1992, A model for the global variation in oceanic depth and heat flow with lithospheric age: Nature, v. 359, p. 123–129. Su, W., and Dziewoński, A.M., 1997, Simultaneous inversion for 3-D variations in shear and bulk velocity in the mantle: Physics of the Earth and Planetary Interiors, v. 100, no. 1-4, p. 135–156. Takeo, A., Kawakatsu, H., Isse, T., Nishida, K., Sugioka, H., Ito, A., Shiobara, H., and Suetsugu, D., 2016, Journal of Geophysical Research : Solid Earth Seismic azimuthal anisotropy in the oceanic lithosphere and asthenosphere from broadband surface wave analysis of OBS array records at 60 Ma seafloor: , p. 1927–1947, doi: 10.1002/2015JB012429.Received. Takeuchi, N., 2012, Detection of ridge-like structures in the Pacific Large Low-Shear-Velocity Province: Earth and Planetary Science Letters, v. 319-320, p. 55–64, doi: 10.1016/j.epsl.2011.12.024. Takeuchi, N., 2007, Whole mantle SH velocity model constrained by waveform inversion based on three-dimensional Born kernels: Geophysical Journal International, v. 169, no. 3, p. 1153–1163, doi: 10.1111/j.1365-246X.2007.03405.x. Tan, Y., and Helmberger, D. V., 2007, Trans-Pacific upper mantle shear velocity structure: Journal of Geophysical Research: Solid Earth, v. 112, no. 8, p. 1–20, doi: 10.1029/2006JB004853. Tarantola, A., 2005, Inverse Problem Theory and Methods for Model Parameter Estimation: Society of Industrial and Applied Mathematics. Tauzin, B., Debayle, E., and Wittlinger, G., 2010, Seismic evidence for a global low-velocity layer within the Earth’s upper mantle: Nature Geoscience, v. 3, no. 10, p. 718–721, doi: 10.1038/ngeo969. Taylor, W.R., 1998, An experimental test of some geothermometer and geobarometer formulations for upper mantle peridotites with application to the thermobarometry of fertile lherzolite and garnet websterite: v. 172, no. 2/3, p. 381–408, doi: 10.1127/njma/172/1998/381. Toffelmier, D.A., and Tyburczy, J.A., 2007, Electromagnetic detection of a 410-km-deep melt layer in the southwestern United States.: Nature, v. 447, no. 7147, p. 991–4, doi: 10.1038/nature05922. 170 Turcotte, D.L., and Oxburgh, E.R., 1967, Finite amplitude convective cells and continental drift: Journal of Fluid Mechanics, v. 28, p. 29–42. Turcotte, D.L., and Schubert, G., 2002, Geodynamics, in Geodynamics, The Press Syndicate of the University of Cambridge, Cambridge, UK, p. 187. Vinnik, L., 2003, Super-deep low-velocity layer beneath the Arabian plate: Geophysical Research Letters, v. 30, no. 7, p. 1415, doi: 10.1029/2002GL016590. Vinnik, L., and Farra, V., 2007, Low S velocity atop the 410-km discontinuity and mantle plumes: Earth and Planetary Science Letters, v. 262, no. 3-4, p. 398–412, doi: 10.1016/j.epsl.2007.07.051. Vinnik, L., and Farra, V., 2002, Subcratonic low-velocity layer and flood basalts: v. 29, no. 4, p. 2–5. Watson, G.S., 1982, Distributions on the circle and sphere: Journal of applied probability, v. 19, p. 265–280. Watts, A.B., 2001, Isostasy and flexure of the lithosphere:. Wittlinger, G., and Farra, V., 2007, Converted waves reveal a thick and layered tectosphere beneath the Kalahari super-craton: Earth and Planetary Science Letters, v. 254, no. 3-4, p. 404–415, doi: 10.1016/j.epsl.2006.11.048. Woodhouse, J.H., and Dziewoński, A.M., 1984, Mapping the upper mantle: three-dimensional modeling of Earth structure by inversion of seismic waveforms: Journal of Geophysical Research, v. 89, no. B7, p. 5953–5986. Xu, Y., and Wiens, D. a., 1997, Upper mantle structure of the southwest Pacific from regional waveform inversion: Journal of Geophysical Research, v. 102, no. B12, p. 27439, doi: 10.1029/97JB02564. Yoder, H. S., J., 1976, Generation of Basaltic Magmas: National Academy of Sciences, Washington, D. C. Youngs, B., and Bercovici, D., 2009, Stability of a compressible hydrous melt layer above the transition zone: Earth and Planetary Science Letters, v. 278, no. 1-2, p. 78–86, doi: 10.1016/j.epsl.2008.11.024. Yuan, H., and Romanowicz, B.A., 2010, Lithospheric layering in the North American craton.: Nature, v. 466, no. 7310, p. 1063–8, doi: 10.1038/nature09332.

Abstract (if available)

Abstract The objective of this thesis is to study the structure of the lithosphere and underlying asthenosphere as well as the correlation of this structure with its overlying tectonic regions. We use several approaches to examining this structure. The first part involves examining the average aspherical velocity structure beneath the Earth’s tectonic regions and determining convergence depths below which the velocity profiles of these regions become indistinguishable from each other. The second part involves an in-depth examination of the Earth’s oceanic regions and their underlying thermal structure. We examine the convergence depths of the oceanic regions for a finer tectonic ocean age regionalization (OAR1) in order to determine related convergence depths, as well as average travel time models for each region in order to examine the thermal structure of that region while minimizing vertical smearing bias which is introduced when examining a region based simply on the underlying velocity structure. The third part uses garnet lherzolite xenoliths to statistically test whether we can reject the null hypothesis that the xenolith’s density changes with depth in accordance with the isopycnic hypothesis, which has implications regarding the structure of cratonic continental lithosphere and asthenosphere. All of these approaches allow us to examine the LAB hypothesis as it applies to both oceanic and continental regions, as well as comparing models of the thermal structure of the oceanic asthenosphere.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Structural clustering analysis of CVMS-4.26: a 3D seismic velocity model for southern California

PDF

Tectonic and kinematic evolution of the Alborán domain during late orogenic extension

PDF

Stress-strain characterization of complex seismic sources

PDF

Constraints on tectonic collisions from seismic time series data

PDF

Structural, emplacement, and tectonic history of spatially and temporally linked volcanic, hypabyssal, and plutonic units in the Beartrap Lake area: eastern Sierra Nevada, CA

PDF

Quantifying ground deformation of large magnitude earthquakes using optical imgaging systems

PDF

Deformation patterns and subduction behavior of continental lithosphere: effects of rheology, surface processes, and material parameters

PDF

Multi-scale imaging of the fault zone velocity structure: double-difference tomography, inversion of fault zone headwaves, and fault zone sensitivity kernels

PDF

Applying automated techniques to large seismic datasets for systematic analyses of phases, source, and structure

PDF

Structural and thermobarometric constraints on the exhumation of the northern Snake Range metamorphic core complex, Nevada

PDF

The role of slab geometry and rheology in inducing net rotation and dynamic topography

PDF

On the hunt for cryptic thrusts: hot rocks, low geothermal gradients, and elusive structures in the eastern Great Basin, Nevada

PDF

The impact of faulting and folding on earthquake cycles in collisional orogens

PDF

Concentration and size partitioning of trace metals in surface waters of the global ocean and storm runoff

PDF

Mapping the Moho in southern California using P receiver functions

PDF

Evolution of lithospheric architecture in arc orogens

PDF

Detection and modeling of slow slip events as creep instabilities beneath major fault zones

PDF

Spatial and temporal evolution of magmatic systems in continental arcs: a case study of dynamic arc behaviors in the Mesozoic Sierra Nevada, California

PDF

Orogenic burial and slow unroofing of the Humber Arm Allochthon constrained by zircon and apatite (U-Th)/He thermochronology, Newfoundland Appalachians

PDF

Integration and validation of deterministic earthquake simulations in probabilistic seismic hazard analysis

Asset Metadata

Creator Paulson, Elizabeth M. (author)

Core Title Deep structure of Earth's tectonic regions

Contributor Electronically uploaded by the author (provenance)

School College of Letters, Arts and Sciences

Degree Doctor of Philosophy

Degree Program Geological Sciences

Publication Date 07/27/2016

Defense Date 06/10/2016

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag asthenosphere,continental,convergence depths,isopycnic,LAB,lithosphere,OAI-PMH Harvest,oceanic,tectonic correlation,thermal

Format application/pdf (imt)

Language English

Advisor Jordan, Thomas H. (committee chair), Lynett, Patrick (committee member), Sammis, Charles (committee member)

Creator Email empaulso@usc.edu,empaulson7@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c40-282116

Unique identifier UC11279611

Identifier etd-PaulsonEli-4640.pdf (filename),usctheses-c40-282116 (legacy record id)

Legacy Identifier etd-PaulsonEli-4640.pdf

Dmrecord 282116

Document Type Dissertation

Format application/pdf (imt)

Rights Paulson, Elizabeth M.

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA