Volume 126, Issue 1 e2020JB020545
Research Article
Free Access

Evidence for the Innermost Inner Core: Robust Parameter Search for Radially Varying Anisotropy Using the Neighborhood Algorithm

J. Stephenson

Corresponding Author

J. Stephenson

Research School of Earth Sciences, The Australian National University, Canberra, ACT, Australia

Correspondence to:

J. Stephenson,

[email protected]

Search for more papers by this author
H Tkalčić

H Tkalčić

Research School of Earth Sciences, The Australian National University, Canberra, ACT, Australia

Search for more papers by this author
M. Sambridge

M. Sambridge

Research School of Earth Sciences, The Australian National University, Canberra, ACT, Australia

Search for more papers by this author
First published: 07 December 2020
Citations: 18

Abstract

The model of cylindrical anisotropy in the inner core (IC) states that seismic rays traveling parallel to the Earth's rotational axis travel faster than those parallel to the equator. There have been continuing discrepancies in estimates of the strength and orientation of anisotropy, with some evidence suggesting that such a model may not be supported by available data. Here, we scrutinize the radial dependence of anisotropy within the IC, where the nature of anisotropy has been shown to change anywhere between a 300 and 800 km radius. We use recent travel time data from the International Seismological Centre in conjunction with the neighborhood algorithm to provide a robust means of testing this idea, through examination of an ensemble of models that satisfactorily fit the data. This can be done with no explicit regularization and without the need for subjective choices associated with binning of phase data. In addition, uncertainty bounds are calculated for anisotropic parameters using a likelihood ratio approach. We find evidence to suggest that commonly employed spatial averaging (binning) methods may be detrimental to obtaining reliable results. We conclude that there is no significant change in the strength of anisotropy with depth in the IC. Instead, we find a change in the slow direction of anisotropy to 54° within the innermost IC at an ∼650 km radius with fast direction parallel to the Earth's rotational axis.

Key Points

  • The “neighborhood algorithm” provides a robust methodology for testing layered anisotropy in the Earth's inner core

  • There is no significant change observed in the strength of anisotropy with depth in the inner core

  • The innermost inner core is defined as a gradual change in the slow propagation direction of anisotropy

1 Introduction and Motivation

The Earth's solid inner core (IC) remains one of the most enigmatic parts of our planet, despite numerous studies in the fields of seismology, geodynamics, mineral physics, and materials science. Making up only 1% of the Earth's total volume with a temperature of over 5,000°C, understanding the structure, current dynamics, and evolution of the IC is vital to understanding the Earth's thermal history (see Tkalčić [2017] for an in-depth discussion). In particular, this provides insights into the geodynamo, without which life would not exist as we know it today. Since its discovery, new models for the structure of the IC are being proposed with increasing complexity in terms of both seismic velocity and attenuation.

Discussions on IC anisotropy arose following a study by Poupinet et al. (1983), who observed that a laterally homogeneous velocity model does not fit the global PKP differential travel time data. Later, Morelli et al. (1986) noticed that absolute PKIKP phases travel up to 3 s faster parallel to the Earth's rotational axis. Anisotropy, rather than isotropic heterogeneity was confirmed in the same year by Woodhouse et al. (1986), who used the hypothesis of IC anisotropy to explain the previously enigmatic observation of the IC sensitive normal modes (e.g., 10S2). The seismic body waves sensitive to the IC structure and the corresponding travel time data used predominantly in IC research are illustrated in Figure 1.

Details are in the caption following the image

(a) Example seismogram showing the arrivals of PKPdf (PKIKP), PKPbc, and PKPab for an epicentral distance of 152° for an event from the South Sandwich Islands recorded in Alaska. (b) Theoretical travel time curves for the three IC phases in terms of epicentral distance. (c) Ray paths of PKPdf, PKPbc, and PKPab (labeled) through the Earth for an epicentral distance of 152°. ξ denotes the angle between the tangent to the PKPdf ray and the Earth's rotational axis. The inner core boundary is denoted by ICB. IC, inner core; ICB, IC boundary.

Research observing anisotropy variations in the upper part of the IC has been abundant (Cao & Romanowicz, 2004; Creager, 1992; Garcia, 2002; Mattesini et al., 2010; Romanowicz et al., 1996; Shearer, 1994; Tanaka & Hamaguchi, 1997; Waszek & Deuss, 2011; Song & Helmberger, 1998), but since 2002, research has also focused on the central part of the IC—the innermost IC (IMIC). This followed observations suggesting a change in the anisotropic elastic constants at a radius of 300 km within the IC (Ishii & Dziewoński, 2002). A shallower transition of elastic properties may be attributed to changes in the solidification fabric of the IC, which is symbiotically coupled to convection systems in the outer core (OC) (Song & Helmberger, 1998; Stroujkova & Cormier, 2004). See supplementary information (Figure S1) for an illustration on the initial observation of an IMIC.

Calvet et al. (2006) were the first to scrutinize the IMIC hypothesis. They showed that their data can be explained equally well by three families of models exhibiting varying IMIC structures, which they partially attribute to poor IC sampling of PKIKP. Other definitions of the IMIC place the discontinuity at 400 (Beghein & Trampert, 2003), 500 (Cao & Romanowicz, 2007), and 430 km (Calvet et al., 2006). Furthermore, Frost and Romanowicz, (2019) infer an IMIC radius at 750 km, but state that a sharp change in anisotropy is not required by the data. In addition to the strength of anisotropy, Ishii and Dziewoński (2002) also inferred a change in slow direction in their IMIC at 45° to the polar axis. A similar result was obtained by Cao and Romanowicz (2007), with an IMIC slow direction between 50° and 55°. Note that along with anisotropy, a variation in the IC attenuation with depth has also been observed to decrease after 300 km from the ICB (Li & Cormier, 2002). Contradicting results have been obtained by some studies which interpreted prominent features that resemble PKIIKP and PKIKPPKIKP waves in global correlograms as the corresponding “reconstructed” body waves (Wang et al., 2015; Wang & Song, 2018). They found the fast direction of anisotropy to be quasi-parallel to the equator; however, S. Wang and Tkalčić (2020) have recently demonstrated that those features in the correlogram are not reconstructed PKIIKP and PKIKPPKIKP body waves.

For IMIC studies, authors utilize the data set acquired from the International Seismological Centre (ISC)—a database of seismic arrivals gathered from seismological institutions globally (Bondár & Storchak, 2011). The catalog's merit is the large available number of PKIKP travel time measurements spanning many decades. New location and quality control measures are constantly being implemented to improve the integrity of the data. The main criticism of the data lies in the large spread of travel time residuals—a result of the data being the product of different phase picking skills, variability in instrument response, or clock errors. Often, the higher amplitude phase PKPbc is picked in error, a feature abundant throughout the travel time catalog. Although 20% of the data is reviewed, this problem still prevails, and it is at the discretion of the researcher to discard data they feel is unreliable. Nonetheless, the numerous data in the ISC catalog have been shown to provide a solid foundation for IC research, especially at the very center (e.g., Ishii & Dziewonski, 2003; Shearer, 1994).

Spatial binning has become a standard practice in IC research, used to facilitate visualization or to suppress noise and thereby improve inversion recovery of the IC structure (Calvet et al.,2006; Creager, 1999; Huang et al., 2015; Ishii & Dziewoński, 2002; Ishii et al., 2002; Mattesini et al. 2010; Shearer 1994; Song, 1996; Stroujkova & Cormier, 2004). The binning process varies dramatically between studies, yielding multiple different findings. Calvet et al. (2006) have shown that there is a stark contrast in resolved IC anisotropy parameters between geometric and geographical averaging methods. However, it is unclear whether these differences are within the uncertainty limits on parameters, especially as uncertainty may be underestimated due to regularization in the inversion algorithm. Overly strict data selection criteria can also have an impact on the results, leading to different models created from the same seismic data (Su & Dziewonski, 1995).

The prevailing issue of data sparsity in both global and azimuthal coverage is a limitation of all studies and creates difficulty in moving forward with IC research. As a consequence, additional uncertainty is introduced by the need to determine appropriate geometric bin sizes in order to address data sparsity issues. Ishii and Dziewoński (2002) observed an IMIC at 300 km using 60 geometrically averaged bins for the whole IC. They noted an anomalous parabolic deviation from the line of constant anisotropy at epicentral distances between 170 and 180°. Our study shows that by employing a direct parameter search approach, spatial binning may be avoided altogether. Furthermore, the IMIC can be defined by a change in the nature and slow direction of anisotropy, and not a change in strength.

2 Theory

Following the convention of Morelli et al. (1986) and Ishii and Dziewoński (2002), travel time deviations from propagation through a cylindrically anisotropic media can be represented by the following:
urn:x-wiley:21699313:media:jgrb54590:jgrb54590-math-0001(1)
where ξ is the angle between the Earth's rotational axis and the vector tangent of the bottoming point of the PKIKP ray, and σ and ε are two anisotropic parameters which act as controls for p wave velocity and can be inverted when modeling for cylindrical anisotropy. They relate to the five independent elastic coefficients as defined by Love (1927). A baseline shift term, urn:x-wiley:21699313:media:jgrb54590:jgrb54590-math-0002, is also introduced to account for uncertainties in the 1D Earth reference model. This term is often small to negligible. In this study, residuals are plotted as a function of urn:x-wiley:21699313:media:jgrb54590:jgrb54590-math-0003 as it enables the anisotropy curve to be viewed conveniently as a parabola.
Relative perturbations in wave speed are related to the PKIKP absolute, or differential travel time residual, urn:x-wiley:21699313:media:jgrb54590:jgrb54590-math-0004, by integrating over raypath, s as
urn:x-wiley:21699313:media:jgrb54590:jgrb54590-math-0005(2)
Combining Equations (1) and (2), we can pose the following equation representing one layer of IC anisotropy. Parameters σ and ε can be solved for as a function of urn:x-wiley:21699313:media:jgrb54590:jgrb54590-math-0006 as
urn:x-wiley:21699313:media:jgrb54590:jgrb54590-math-0007(3)
This formulation can be extended to two layers, where the integral limit, H, is also unknown as follows:
urn:x-wiley:21699313:media:jgrb54590:jgrb54590-math-0008(4)
where σ1 and ε1 are the anisotropy parameters for the top layer (outer IC [OIC]) and σ2 and ε2 the bottom layer (IMIC). This model format gives the free parameter, H, to be the limit of the integrals, representing the radius to change in anisotropy. Because of this unknown parameter, the system of equations is now nonlinear. ICB is the IC boundary, set at 1,225.7 km, and B the bottoming point of each ray in the data set. For details on the formulation of the anisotropy equations and calculation of the angle ξ, see Tkalčić (2017).

3 Data and Processing

3.1 Absolute PKIKP Travel Time Data Set

To investigate the presence of an IMIC and address the effect of geometric binning, this study utilizes absolute PKIKP travel times obtained from the Reviewed ISC catalog (Bondár & Storchak, 2011). Here, absolute PKIKP travel times are downloaded between the years 2001 and 2013 to give a total of 183,257 rays. Note that Ishii and Dziewonski (2003) use ISC data from 1964 to 1994 and, although a large time period, there has since been considerable improvements made to the location algorithm, particularly for depth resolution. From 2001 onwards, earthquake locations utilized additional arrival times from S, Sg, Sn, and Sb arrivals.

To avoid complex phases associated with the PKP triplication and contamination from precursors, all data with an epicentral distance <150° are discarded. To further remove outliers and reduce noise, data are processed in 10 epicentral distance bins and the distribution of each bin is analyzed. In instances where the data are contaminated with PKPbc mis-picks, a bimodal distribution is prevalent between 150° and 159°, so clusters representing mis-picks are discarded. Additional scatter is reduced by cutting the data at 3 standard deviations on either side of the mean in each bin. This preserves 99% of the data, while reducing the overall spread. The data were analyzed to make sure that a good number of data exist in all epicentral distance ranges and propagation angles. Most discarded data were from epicentral distances less than 165° at intermediate values of ξ, as these geometries typically have a larger number of data. We noted no dependence on the bottoming location in the IC for the discarded data. Despite careful processing, the data are still limited in spatial distribution, especially at polar antipodes. After full processing, the complete data set contains 33,859 PKIKP ray paths through the IC, giving over 1,000 paths in the range of 172°–180° (bottoming depths in the data range from 20 to 1,000 km radius); however, there are no true antipodal paths sampling the center of the IC. Critically, of these 1,000 antipodal paths, less than 30 are considered polar (ξ < 35°), which could have repercussions for resolving the strength of the anisotropy in the IC below 300 km radius. The data are shown spatially in Figure S2.

All data are corrected for the Earth's ellipticity (Engdahl et al., 1998, and references therein). Figure 2 shows the final data set compared to the original data. A mantle correction is made post data-reduction for all data using the model of Simmons et al. (2012). A more detailed discussion is provided in the Section 1, Figures S3 and S4.

Details are in the caption following the image

Absolute travel time residuals relative to ak135 plotted against the angle ξ for (a) unprocessed data and (b) processed data according to the method described in the main text, from the ISC catalog between 2001 and 2012. ISC, International Seismological Centre.

3.2 Estimating Data Noise

The complete distribution and magnitude of uncertainties in the ISC data set are unknown. However, using a degree of informed judgment, uncertainties can be approximately quantified to aid interpretation and propagate errors into models of IC anisotropy. Uncertainties exist in terms of theoretical error—our reliance on accurate source locations and a velocity model of the Earth or the assumption of ray theory, and experimental error—analyst mis-picks, or instrument timing. Table 1 shows a consolidated summary of known uncertainties in the ISC data. For further information on the bias and uncertainty in the ISC catalog, see Rohm et al. (1999).

Table 1. Uncertainties in the ISC Absolute PKIKP Travel Time Data
Uncertainty Value Details
Phase misidentification Up to 5 s <5% of data, more likely at distances <155°
Inconsistent picks/timing errors ±0.5–1 s All data, >1 s for 8% of ISC P arrivals Rohm et al. (1999)
Repeat measurements ±0.2 s Extrapolated to all data
Crust/mantle structure ±0.4 s Higher for shallow earthquakes
CMB topography ±0.5–1 s All data
IC heterogeneity Up to 3 s Western hemisphere polar paths
Derived from synthetic testing ±1.5 s All data—random
  • Abbreviations: IC, inner core; ISC, International Seismological Centre.

To visualize the effect of noise on the data, synthetic data is created using the true depth, epicentral distance, and ξ geometry of the ISC data for varying models of anisotropy. Random Gaussian noise with standard deviations up to 2.0 s is added to the synthetic data (see Figure 3), where anisotropy follows Ishii and Dziewoński (2002)'s 1.8% in the OIC and 3.7% in the IMIC). Within the IMIC data (blue crosses), there is a clear parabolic shape with a peak at around urn:x-wiley:21699313:media:jgrb54590:jgrb54590-math-0009, which is an effect of a large negative σ value of −19.6. This effect is also visible, to some extent, in the real data. As noise levels are increased, these details are obscured in the data; however, the synthetic data start to appear similar to real data. In particular, there is a distinctive “fish tail” shape (e.g., panels f, h, and j) where travel time residuals are smaller between urn:x-wiley:21699313:media:jgrb54590:jgrb54590-math-0010  = 0.1 and urn:x-wiley:21699313:media:jgrb54590:jgrb54590-math-0011  = 0.2. This feature can be seen in the real data (panel b) and is likely an effect of the spatial distribution as well as the fact that there are fewer observations in this range. The synthetic data most resembles the shape and distribution of the real data at a noise level of 1.5 s standard deviation. Therefore, a reasonable model for observational uncertainty is a Gaussian with a standard deviation of 1.5 s.

Details are in the caption following the image

Comparison of real and synthetic data for the PKIKP ISC absolute travel time data set. The OIC has 1.8% and the IMIC has 3.7% cylindrical anisotropy with IMIC slow direction aligned at ∼45°to the rotational axis. The IMIC radius is set to 300 km. For all panels, gray crosses show absolute travel time residual data for the bulk of the IC, blue crosses show data for rays that penetrate the IMIC, and yellow crosses are the subset of SSI to Alaska ray paths. The left set of panels show data plotted as a function of epicentral distance and the right set of panels show the data plotted as a function of cos2ξ. Panels (a and b) show real data, (c and d) show synthetic data with minimal Gaussian noise (0.01 s standard deviation), (e and f) 0.5 s, (g and h) 1.0 s, (i and j) 1.5 s, (k and l) 2.0 s standard deviation noise, respectively. IC, inner core; IMIC, innermost IC; ISC, International Seismological Centre; OIC outer IC; SSI, South Sandwich Islands.

4 Methods: The Neighborhood Algorithm

To invert for anisotropic parameters, and nonlinear parameter H, the derivative-free direct-search algorithm, the neighborhood algorithm (NA) is implemented (Rickwood & Sambridge, 2006; Sambridge, 1999). The NA has been used once before in IC anisotropy studies for normal mode data (Beghein & Trampert, 2003), and has been widely applied to other nonlinear inverse problems in geophysics (e.g., Pachhai et al., 2016). Parameter values are obtained by constructing Voronoi cells in the parameter space (Okabe et al., 2000) around lower misfit models, with each iteration guiding the concentration of sampling in the next. The algorithm is relatively straightforward to control, with three user-determined tuning parameters, which are the starting number of models, the number of models created per iteration, and the number of lower misfit models chosen to guide subsequent iterations. Uniform random walks are performed within the chosen Voronoi cells through multidimensional parameter space. By selecting a large number of initial models and then choosing a high number of lowest misfit models the algorithm becomes highly exploratory of model space, which is useful in cases where several local minima may exist. The result of the algorithm is an ensemble of models which fit the data to varying degrees, requiring no explicit regularization. Parameter bounds are chosen based on the maximum plausible values of anisotropy from mineral physics studies. ε values are set at ±0.05 (5% anisotropy), σ values ±20, and H is set between 200 and 1,100 km. A more detailed discussion on the choice of these prior bounds is provided in the supplement (Section S2. During the algorithm, the L2 norm misfit function is calculated, and minimized with each iteration (see supplement Section 3).

5 Experiments With Synthetic Data

Although numerous, IC data from the ISC catalog are noisy and sparse in terms of overall IC sampling. Synthetic tests are carried out for several known IC parameterizations to investigate resolvability at varying noise levels, anisotropy contrasts, and bin sizes. Synthetic testing also provides an opportunity to tune the NA search parameters for optimal convergence, allowing the algorithm to thoroughly sample the parameter space. In general, the parameters were set to be highly exploratory, with a large number of cells being resampled at each iteration. This reduces the tendency of the algorithm to become “stuck,” due to the “choking” phenomenon identified in Sambridge (2001), which is a possibility with such a low-dimensional parameter space. After a large number of tests on differing synthetic data, the optimal NA search parameters were chosen as 1,000 initial samples, 1,000 samples, and 1,000 cells resampled. The algorithm ran for 50 iterations and convergence was met when the misfit value no longer decreased with subsequent iterations. In addition, multiple initial random seeds were given to the algorithm to vary the search starting position. This made little difference to the results obtained as long at the algorithm had enough iterations. All synthetic tests were carried out on data with added random Gaussian noise with a standard deviation of 1.5 s, consistent with the noise level in the true data. Figure 4 shows the results of three contrasting synthetic tests for the five parameters, using the NA. Data have not been binned in these experiments.

Details are in the caption following the image

Results from running the NA inversion on three contrasting, synthetic data sets. The upper row shows the tradeoffs in 3D between OIC anisotropic parameters ε1 and σ1 with H (IMIC radius - vertical axis). The lower row shows the tradeoffs in 3D between IMIC anisotropic parameters ε2 and σ2 with H. (a and d) are for a homogeneously anisotropic IC at 3%, (b and e) for the IMIC model of Ishii and Dziewoński (2002), and (c and f) for an IC decreasing in anisotropy from 3.2% to 1.8% at a radius of 500 km. In each case the vertical axis is H. The colored points each represent a single model output from the NA at three χ2 confidence levels: green—68%, yellow—95%, and pink—99%. All points have also been projected onto the corresponding 2D axis for clarity (white points). The black triangle marks the values of synthetic input parameters and the black circle, the best fit model from the NA. A triangle is not present for (a and d) as the synthetic model is a single layer. IC, inner core; IMIC, innermost IC; OIC outer IC; NA, neighborhood algorithm.

Let us first examine the case of an anisotropic IC with only one layer of anisotropy, which remains constant at 3% with increasing depth within the IC. The results are presented in Figure 4a (OIC parameters) and Figure 4d (for IMIC parameters). The clouds of colored points within the 3D plots represent the best solutions from the NA at the 99%, 95%, and 68% confidence intervals (see supplement Section 3 for calculation of confidence intervals). For clarity, 2D versions of this plot are also provided in the supplementary material (Figures S6–S8). Comparing the size of model ensembles, it is clear that the OIC is better resolved than the IMIC. Here, the NA has placed H at the very edge of the parameter space at a 200 km radius, converging to an IC with a single layer of anisotropy. Here, ε and σ values in the OIC are also close to the true synthetic parameters in the isotropic input, meaning that they have been correctly recovered. Although a two-layer IC is parameterized, if the IC is best fit with one layer of anisotropy, the algorithm is able to converge to a radius at the very edge of the parameter space.

The second synthetic test shown uses the IMIC parameterization of Ishii and Dziewoński (2002) (Figures 4b and 4e). Here, the ensemble of models at each confidence level is smaller for the OIC parameters, meaning that the solution for the OIC is more constrained. Anisotropic parameters from the lowest misfit model closely approximate the synthetic parameters. Additionally, the location of H is close to the true value of 300 km, despite only 3.4% of rays in the data set penetrating to this depth. This result is encouraging for using the real data to identify an IMIC.

Our final test is shown in Figures 4c and 4f. The synthetic input here was a 3.2% anisotropic OIC and a weaker anisotropic IMIC at 1.8%. This particular test was chosen to assess the ability of the IMIC structure to be resolved by the NA when over-printed by a stronger OIC. In comparison to the previous example, the location of H is more poorly defined. The large ensembles of models at each confidence level suggest that the boundary has equal likelihood of being placed anywhere below 600 km, or above 900 km, the latter value suggesting that a homogenously anisotropic IC may be preferred by the data. The ensembles of models for the IMIC are equally large as the previous test, even though the radius is higher at 500 km with around 11% of rays penetrating to this depth. However, the lowest misfit model is still close to the true model, with the exception of ε2, which deviates more. For all three synthetic cases, the total noise level found using the procedure outlined in the supplement Section 3 was 1.5 s for all three cases. In this synthetic test, involving no additional sources of theoretical error, the correct data noise level was recovered.

Our synthetic test lends encouragement to the suitability of the NA for this particular problem, especially with regard to its performance in resolving known anisotropic parameters. Despite the uncertainties being large, the NA converges on the lowest misfit values close to known inputs in all cases. However, the algorithm and data together do not have the ability to distinguish a sharp discontinuity in anisotropy, due to the noise level or address any smaller velocity heterogeneities within the IC.

5.1 Resolvability of H

Model recoverability may be limited for a deep IMIC radius (small H) due to poor spatial sampling and limited azimuthal coverage at these depths. Therefore, several experiments with a varying H are carried out for assessing the recoverability of each depth. For the first test, synthetic data with weak contrast between OIC and IMIC was chosen, that being the least well resolved model in previous testing. The results from this test are shown in Figures 5a and 5b, along with uncertainty estimates based on the output ensemble of models from the NA. The results show that when the radius is at 300 km, the uncertainty on the IMIC parameters is the highest. This is not surprising due to the lack of antipodal, polar data sampling the IMIC at a 300 km radius. Likewise, Beghein and Trampert (2003) note that their anisotropic parameters from normal mode studies have large errors deeper within the IC, which they also attribute to sampling issues. The results further show that OIC parameters have more uncertainty when H is close to the ICB (at 1,100 km), which is also foreseeable due to rays spending only about 100 km in the OIC. In spite of large uncertainties, the best fit model parameters are all close to the true input parameters of the synthetic data. Furthermore, anisotropic parameters are most constrained at larger IMIC radii (>650 km). Therefore, results using real data resolving the IMIC at less than 650 km radius must be interpreted with caution, especially for ε2, where values within the 68% confidence level ensemble vary dramatically, each providing a different picture of the IMIC. As may be hypothesized, synthetic tests for a larger contrast between two layers of anisotropy in the IC show similar results, but tighter regions of uncertainty.

Details are in the caption following the image

Top panels: Results of an NA inversion of synthetic data with H (IMIC radius) fixed at five radii within the IC for (a) ε and (b) σ parameters in each layer. Each dashed line shows the best-fit parameter from the NA for each radius inverted for. Solid gray lines show the synthetic input parameter values. The bounds around the best fit parameter values (indicated by red- and blue-shaded area) show the ensemble of models at the 68% confidence interval for the location of the best fit IMIC radius. Input anisotropic parameters are ε1 = 3.48, σ1 = −4.95, ε2 = 1.82, and σ2 = −6.75. Bottom panels: Relative accuracy between nonbinned to binned data for a varying number of bins. Values are relative to the accuracy of resolved true parameters of nonbinned data. The gray zero line represents the accuracy of the nonbinned data, thus points above represent less accurate resolved values and below, more accurate (closer to true input parameters). Results are shown for two synthetic models: (c) ε1 = 3.48, σ1 = −4.95, ε2 = 1.82, σ2 = −6.75, H = 500 km and (d) ε1 = 1.80, σ1 = −0.67, ε2 = 3.70, σ2 = −19.7, H = 300 km. IC, inner core; IMIC, innermost IC; NA, neighborhood algorithm.

5.2 Binning of Synthetic Data: To Bin or Not to Bin?

To test the level of binning necessary for robust results, binning is carried out in a simplistic manner and performed on synthetic data with a noise level of 1.5 s. Each ray in the data is divided into bins based on the epicentral distance followed by being divided into subsets based on values of cos2ξ. Within each bin, an average earthquake depth was also calculated from the events of rays in that bin. This process resulted in each bin representing a summary ray for a particular epicentral distance, ξ, bottoming depth, and event depth. The mean and standard deviation travel time residuals in each bin are calculated and travel times from ak135 can then be obtained using this new summary ray geometry. In terms of testing the recovery of known parameters for varying bin size, initial numbers of six epicentral distance bins, and 10 cos2ξ bins are used to replicate the study of Ishii and Dziewoński (2002). This equates to the epicentral distance of bins spanning 5° and cos2ξ bins of 0.1 in size. Following this, the number of either epicentral distance or cos2ξ bins is doubled systematically, which doubles the total number of summary rays.

To compare the binned results with results from data that have not been binned, accuracy tests are carried out providing a quantitative assessment on whether binning enhances or degrades recovery of the IMIC structure. The following equation is used to assess the accuracy of the binned solutions in comparison to nonbinned data (Beiranvand et al., 2017):
urn:x-wiley:21699313:media:jgrb54590:jgrb54590-math-0012(5)
where the relative accuracy of the solution, xacc, is expressed as the difference between the known solution, x, and the binned solution urn:x-wiley:21699313:media:jgrb54590:jgrb54590-math-0018, normalized for the nonbinned solution x0. The logarithm of the ratios is taken to improve the readability of the results. The relative accuracy, xacc, is shown for all five parameters for both a strong (Figure 5c) and a weak contrast in anisotropy (Figure 5d). The gray line represents the nonbinned data—values below this line indicate that the solution for binned data is more accurate, and above, less accurate. It also follows that points plotting on the gray line are of equal accuracy to the nonbinned results, therefore in these instances, binning the data has made little difference to the solution. When the results are presented in this framework, the difference between the two test cases is stark. For a weak decrease in anisotropy (Figure 5c) the solution most accurate to the known parameters is the data set with the most bins, a feature that holds true for all five model parameters. The opposite is the case for the strong contrast in anisotropy (Figure 5d), where IMIC parameters are more accurate with a larger number of bins, but OIC parameters are less accurate. The reason for this is not clear, but it can be assumed that the IMIC parameters are better resolved with a larger number of bins—as the number of bins increases so does the number of summary rays penetrating the IMIC. It is clear from Figure 5 that there is no binning case that improves the solution for all IC parameters—there are no cases where all values plot below the gray line. In addition, tests show that the least accurate model is using only 60 bins, suggesting that here data are being over-averaged and further uncertainty is being introduced into the solution. Therefore, we conclude that binning is not necessary and is potentially detrimental to observing the IMIC structure.

5.3 Uncertainty Analysis

Due to the large ensemble of models at each confidence interval from NA inversion, there is difficulty in fully quantifying uncertainty, especially on H. To further address the model robustness, we use a complementary linearized least squares (LLS) method to provide additional uncertainty information.

Consider the radius between two anisotropic layers of the IC, H, to be fixed at a certain known depth. In this case, the travel time integrals urn:x-wiley:21699313:media:jgrb54590:jgrb54590-math-0013 and urn:x-wiley:21699313:media:jgrb54590:jgrb54590-math-0014 in Equation 4, can be computed for the IMIC and OIC, respectively, as the integral limits are known and fixed. Here, the system becomes linear and can be solved using the standard methods of discrete parameter estimation for H between the bounds 200 ≤ H ≤ 1,100. In this inversion, data uncertainties are assumed normally distributed, statistically independent with a uniform variance, which corresponds to a diagonal data covariance equal to the identity matrix multiplied by a factor of 1.52 s—the determined noise level of the ISC data (Figure 3). The results of the inversion give the best fit for four anisotropy parameters at each position of H and their associated least squares misfit.

In the linearized inversion problem, the poor geometric sampling of the IMIC at certain H values may mean that the problem is ill posed. To quantify this, the condition number κ is calculated for the model covariance (Cm) corresponding to each linear parameter estimation problem at fixed H. The condition number is defined as the ratio of the maximum, λmax, to minimum, λmin, eigenvalue of Cm
urn:x-wiley:21699313:media:jgrb54590:jgrb54590-math-0015(6)
Smaller condition numbers are generally preferred, signifying a more well-posed inversion (Aster et al., 2019). While the model covariance matrix provides a linear estimate of uncertainty of the anisotropic parameters for a given value of IMIC radius, H, we also require some estimate of error in H itself. This can be facilitated with the likelihood ratio test (Neyman & Pearson, 1928). If we define the log-likelihood of the least squares solution of the linear system for a given H as (m1H), then the log-likelihood ratio, LR, is written as
urn:x-wiley:21699313:media:jgrb54590:jgrb54590-math-0016(7)
where (m0, HO) is the global minimum, that is, the best data fit of anisotropic parameters and IMIC radius, H, combined. The convenience of this expression is that if the sample size is large, the statistic −2logLR follows an χ2 distribution with one degree of freedom (Wilks, 1938). This permits confidence intervals for H to be calculated and displayed around the global minimum, HO.

Multiple synthetic tests were carried out using this approach to test the feasibility and to further investigate the effects of the H position on the resulting model and uncertainty. Synthetic results show that the higher the condition number of the corresponding model covariance matrix of the anisotropic parameters, the less constrained the position of H. However, the location of the global minimum is always located within 30 km of the true synthetic input (Figure 6). The only exception is for synthetic data with H = 1,100 km, which is unsurprizing as it is at the edge of the parameter space. Furthermore, the inversion is most constrained (low condition number), with a boundary between 800 and 1,000 km, likely due to the abundance of data and good azimuthal coverage in this range. These observations are in good agreement with observations of the NA synthetic tests and lend confidence to both recovery of the IMIC radius and its uncertainty from the available data.

Details are in the caption following the image

Synthetic tests showing constraint on the IMIC radius, H, for five synthetic models with a weak contrast in anisotropy. (a) The condition number (κ) of the model covariance matrix at the least squares solution for each (h) (b–f) The likelihood ratio (10) against H for five synthetic tests with the true value of H indicated by a black arrow. (g) Results from the real data. Colors represent the confidence intervals at 99% (purple), 95% (beige), and 68% (green) determined using the likelihood ratio test as described in the main text. IMIC, innermost IC.

6 Results and Discussion

Up to this point, we have used synthetic data to show that the NA is a robust methodology for identifying the presence of an IMIC using the available geometry of the ISC data set. We have also shown that binning is not necessary and, in some cases, detrimental to observing the structure in the deepest parts of the IC. As we have assumed a two-layer IC in our model, the results discussed in this section do not preclude other structures in the IC, such as hemispherical variations in anisotropy (e.g., Lythgoe et al., 2014). Thus, the NA inversion was carried out on subsets of quasi-western and quasi-eastern hemisphere data independently, but the global coverage was not deemed optimal for reliable interpretation of the results based on additional synthetic testing of these subsets. The lack of data, especially at the polar antipodes may, in part, account for the inconsistencies in the hemispherical structure in the IC when the depth dependence of anisotropy is also considered (see Irving & Deuss, 2011; Lythgoe et al., 2014; Tanaka, 2012). We also applied an anisotropy correction to the upper 400 km of the IC based on the independent higher quality PKPbc-PKIKP differential travel time data of Tkalčić et al. (2002), Leykam et al. (2010), Young et al. (2013) and additional data collected for this study. Values for correction were obtained from a simple least squares fit to the PKPbc-PKIKP data (Figures S10 and S11). The correction parameters were calculated both with and without paths from the South Sandwich Islands (SSI). This tested the influence of anomalous paths from the SSI to Alaska on the IMIC, which may be due to a subducted slab or other heterogeneity in the mantle (Breger et al., 2000; Frost et al., 2020). The correction had a small effect on the strength of anisotropy in the OIC (as it removed most of the anisotropy) but had little effect on values for the IMIC, or the level of uncertainty around parameters. However, the position of H with the correction was higher at ∼800 km. This higher position of H may be due to an artificial boundary being created at 400 km from the ICB, due to the correction applied. Therefore, given our primary aim is investigation of the deep IC, the preferred result does not include this correction. In contrast, applying a mantle correction had a significant effect on reducing the spread of models within each ensemble, but negligible effect on resolved parameters, including H.

Full results for a two layered IC from the NA inversion can be found in Figures 7 and 8 and Table 2. Noticeably, at the 68% uncertainty level, the confidence bounds around ε parameters in both layers significantly overlap, suggesting either a weak or gradational change in anisotropy or a high uncertainty in parameter values. Taking into consideration this range of solutions represents an average of the two hypothetical anisotropic hemispheres, our strength of anisotropy has good agreement with previous work using IC body waves. This includes all three models of Calvet et al. (2006) (1.00%–1.56%), the 320–620 km IC radius range of Su and Dziewonski (1995) (1.6%), and IMIC models of Ishii and Dziewoński (2002) (2%) and Ishii and Dziewonski (2003) (1.8%). Our result has a much lower magnitude of anisotropy than the IC model of Lythgoe et al. (2014), where they observed 1.79% for the 750–1,220 km range, increasing to 3.31% for 320–620 km. However, better agreement is seen with Lythgoe et al.'s (2014) OIC anisotropy value for the eastern hemisphere alone (1.16%).

Details are in the caption following the image

ISC PKIKP travel time data with the final model and three published models shown for comparison, as a function of cos2ξ. Real data are shown as gray circles and models as colored crosses with corresponding models indicated on the legend. Published models that include three or more layers have not been included, nor the models that do not state the exact values of σ for each layer. ISC, International Seismological Centre.

Details are in the caption following the image

 Ensemble of models produced by the NA for all ISC data. In all cases, the circle indicates the location of the best-fit model and ensembles are grouped into three confidence levels as indicated in the legend using a noise level of 1.5 s. Left panels: (a) information on algorithm convergence as described in the legend, (b–e) tradeoffs between anisotropic parameters plotted against the IMIC boundary radius, H, and (f–k) tradeoffs between different sets of anisotropic parameters. ISC, International Seismological Centre; NA, neighborhood algorithm.

Table 2. Parameters for the Preferred Model of Two-Layered Cylindrical Anisotropy in the IC
Parameter Lowest misfit value Solution LSRQ Uncertainty from NA (68%) Uncertainty from LLS (95%)
ε1 1.30% 1.45% 0.3%, 2.2% 1.19%, 1.50%
σ1 −1.11% −1.07% −2.4% −3.32, −0.83%
ε2 2.19% 1.99% −0.9%, 4 % 1.45%, 2.51%
σ2 −7.72% −7.54% −12, −3% −8.41, −6.68%
H 657 km 655 km 560, 880 km 650, 675 km
Slow axis (IMIC) 54° 53° 50°, 68° 52°, 54°
  • Note. Uncertainty is at the 68% (NA) and 95% (LSRQ) confidence intervals with a Gaussian noise level of 1.5 s. Uncertainty on the angle of slow propagation is based on the uncertainty bounds of the urn:x-wiley:21699313:media:jgrb54590:jgrb54590-math-0017 parameter.
  • Abbreviations: IMIC, innermost IC; LLS: linearized least squares; NA, neighborhood algorithm.

Values for the IMIC have a significantly larger region of uncertainty compared with the OIC, as is also observed in synthetic testing. In this layer, ensembles at the 68% confidence level show strengths of anisotropy to be between −0.9% and 4%. Interestingly, the negative value at the lower end of the range indicates that the fast direction of propagation could be aligned with the Earth's equator, which may mean that the IC structure at these depths is more complex than simple cylindrical anisotropy or too obscured by noise. However, over 75% of the models at this uncertainty level have a positive sign. It therefore follows that from these results we cannot draw conclusions on the nature of anisotropy in the IMIC, or whether there is a definite change in the absolute strength at 657 km. This uncertainty is also reflected significantly within previous IMIC studies, where values vary from −1.15% (Calvet et al., 2006) to 5.53% (Lythgoe et al., 2014).

Despite the uncertainty described above, further insights can be gained by closer scrutiny of the σ parameters in each layer, which are far more constrained. In particular, from the NA, the resulting IMIC parameter σ2 is relatively large and negative (−7.72%) with uncertainty between −12 and −3. Calculating the zero root of the derivative of the anisotropy function enables us to obtain the slow direction of propagation (54°) in the IC along with its uncertainty (50°–68°) (see Table 2). In the OIC, the slow direction of propagation is most likely to be parallel with the equator as we observe no maxima in the fuction of the best fit model. There is, however, a small maxima in the anisotropy function at around 50° at the upper uncertainty bounds of the anisotropy function, but this evidence for a slow direction deviating from equatorial in the OIC is weak. This is compelling evidence to suggest that the slow direction undergoes a change in direction to ∼54° relative to the Earth's rotational axis in the IMIC. This change occurs at ∼650 km. Our result is further supported by the abundant data at intermediate angles through the IMIC.

Despite a high uncertainty, our results for the IMIC agree well with prior work, with the majority of prior studies showing a large, negative σ value for the IMIC below a radius of 600 km (Calvet et al., 2006; Frost & Romanowicz, 2019; Garcia & Souriau, 2000; Ishii & Dziewoński, 2002; 2003; Lythgoe et al., 2014). However, Beghein and Trampert (2003) argue for a negative ε in the IMIC, which we only observe at the lower end of the ε2 uncertainty bound. In addition, many experimental and theoretical mineral physics studies support our results for the IMIC, where a slow direction is predicted through both the bcc and hcp phase of iron between 40 and 50° from the fast direction (e.g., Belonoshko et al., 2008; Steinle-Neumann et al., 2001). An increase in the angle with fast and slow directions is also observed with increasing temperature in some studies (e.g., Steinle-Neumann et al. 2001), meaning that an IMIC with a different form of anisotropy is physically plausible.

In terms of the depth to a change in the nature of anisotropy, H, the NA uncertainty range between 560 and 880 km (see Figures  8b–8e) encompasses the values of previous studies (Beghein & Trampert, 2003; Calvet et al., 2006; Cao & Romanowicz, 2007; Frost & Romanowicz, 2019; Garcia & Souriau, 2000; Ishii et al., 2002; Lythgoe et al., 2014) with the exception of Su and Dziewonski (1995) and Ishii and Dziewonski (2003), who both placed H deeper at around 300 km. In both the cases, these studies chose to bin data, which likely had an effect on the result. Cormier and Stroujkova (2005) found little evidence for a sharp IMIC transition and suggested that a structural change within the IC is more likely to be gradual. A gradual change between the OIC and IMIC could account for our large uncertainties in depth to the IMIC boundary.

The results for the ISC data using the complementary linearized inversion scheme show good agreement with the NA parameter search results (Table 1), with the exception of the ε2 IMIC parameter. This parameter was the least well resolved in all synthetic tests for both the NA and LLS, so a disagreement here between the methods is not unexpected. However, the agreement between the σ2 parameters for both the inversion methods and its relatively small uncertainty provides additional support for a change in the fast direction of propagation. Furthermore, the uncertainty on H using this method is constrained tightly around the inversion global minima along with being located at a position of reasonable inversion accuracy (Figure 6g). Uncertainties for this method are likely smaller due to an underestimation of the noise levels or our assumption of a Gaussian error distribution. However, the result in Figure 6g shows a local minima at just over 800 km. If a larger noise value is assumed in the data covariance matrix, this minima is included in the 68% and 95% uncertainty bounds. This minima also explains the spread of H values at each confidence level observed in the NA result. These additional results using complementary methods provide further compelling evidence for a change in the nature of anisotropy at around 600 km. In addition, the direction of slow propagation is also in agreement with the NA result, also with tighter uncertainty bounds.

6.1 Conclusions and Future Directions

In summary, this work has shown that an IC with two layers of differing cylindrical anisotropy is a plausible inference for the IC if a two-layer parameterization with cylindrical anisotropy in each layer is explicitly imposed. This is a significant outcome given the noise level in the data and the sparsity of the ISC data set. Through a rigorous treatment of uncertainty, we concluded that binning the ISC data set is not a favored practice to observing the IMIC structure. While there is no strong evidence for a change in the overall anisotropy with depth, the large, negative IMIC σ2 parameter suggests a change in the slow direction of anisotropy at ∼54° to the Earth's rotational axis. This change is likely gradual somewhere between 600 and 640 km, but we cannot completely rule out the change occurring higher up within the IC. This model has good agreement with recent studies on IMIC anisotropy. In particular, Frost and Romanowicz (2019) observed a slow direction at 60° in the IC anisotropy at radii below 900 km. Their observed effect becomes higher in magnitude with depth, which is likely why we also observe a clear difference in the nature of the IC below 600 km, and note a change below 800 km when the top 400 km of anisotropy is removed. Our results do not rule out further complexity in the structure of the IC or the plausibility of hemispherical variations in anisotropy. Note also that in the inversion, cylindrical anisotropy has also been defined with the symmetry axis parallel to the ERA. If the axis was allowed to vary spatially in the inversion, it is likely that we will observe a tilt in the fast direction in some areas of the IC. However, this does still not overcome issues in data sparsity at polar antipodes.

Based on our results, we can also confidently rule out the existence of anisotropy near the Earth's center with a fast axis parallel to the equator, as proposed in recent studies interpreting the features in global correlograms as reconstructed body waves (Wang et al., 2015; Wang and Song 2018). Although global correlation studies are increasing and have the potential to resolve the long-standing puzzle of IC anisotropy, care should be taken when interpreting correlation features that resemble body waves sensitive to the Earth's IC. Using different approaches, it has been shown recently that the formation of correlation features does not equal the “reconstruction” of body waves between pairs of receivers (e.g., Phạm et al,. 2018; Poli et al. 2017). Moreover, Wang and Tkalčić (2020) present observational proof that the correlation features in global correlograms sensitive to the IC are constructed by cross-terms of a multitude of body waves, whose relationship with the Earth's structure is complex. Therefore, the previous inference that anisotropy near the Earth's center is characterized by a fast axis parallel to the equator was likely a result of oversimplification of the ray paths associated with the correlation features in question.

Unfortunately, given the difficulties in calculating the behavior of hcp and bcc iron at temperatures and pressures of the IC, it is difficult to attribute our results solely to a phase change of iron, although some studies suggest that the bcc phase is responsible for the low attenuation observed in the deeper parts of the IC (Belonoshko et al., 2019). However, disagreement of IC models based on seismic data with hcp iron studies have been attributed to further phase changes at depth or impurities deep within the IC, including the existence of multiple phases of iron (e.g., Beghein & Trampert, 2003; Mattesini et al., 2013; Tkalčić, 2010). There is also evidence from numerical modeling to suggest that a structurally different IMIC may be due to delayed nucleation during IC growth (Lasbleis et al., 2020).

At present, including in this work, we are limited by the distribution of global earthquakes and receivers, especially at the polar antipodes, even in a data set as extensive as the ISC. Polar data for the OIC is dominated by anomalously fast paths from the SSI to Alaska, and deep polar antipodes rely upon noisy data from Antarctic stations. This is the most prevalent issue to date for research on IC anisotropy. However, there is hope that this situation will change with the advent of studies exploiting the Earth's correlation wavefield (e.g., Phạm et al., 2018; Tkalčić & Phạm, 2018). This approach promises the ability to derive deep Earth seismic structure constraints uninhibited by the limited distribution of sources and receivers, potentially in sparse locations globally, significantly increasing the resolution and understanding of the IC and hence Earth's evolution.

Acknowledgments

J. Stephenson is grateful for the support and guidance of the co-authors leading up to the submission of this manuscript—work which formed a part of her PhD research. The authors especially thank Vernon Cormier and an anonymous reviewer for their constructive reviews, which resulted in significant improvements to the manuscript. J. Stephenson was supported by an Australian National University research scholarship.

    Data Availability Statement

    Data were downloaded directly from the server of the International Seismological Centre (ISC) http://www.isc.ac.uk/iscbulletin/search/. The Neighborhood Algorithm code was provided by M. Sambridge (http://www.iearth.org.au/codes/NA/).