Constraining lowermost mantle anisotropy with body waves: a synthetic modelling study

SUMMARY

Different mechanisms have been proposed as explanations for seismic anisotropy at the base of the mantle, including crystallographic preferred orientation of various minerals (bridgmanite, post-perovskite and ferropericlase) and shape preferred orientation of elastically distinct materials such as partial melt. Investigations of the mechanism for D" anisotropy usually yield ambiguous results, as seismic observations rarely (if ever) uniquely constrain a mechanism or orientation and usually rely on significant assumptions to infer flow patterns in the deep mantle. Observations of shear wave splitting and polarities of SdS and PdP reflections off the D" discontinuity are among our best tools for probing D" anisotropy; however, currently available data sets cannot constrain one unique scenario among those suggested by the mineral physics literature. In this work, we determine via a forward modelling approach what combinations of body wave phases (e.g. SKS, SKKS and ScS) are required to uniquely constrain a mechanism for D" anisotropy. We test nine models based on single-crystal and polycrystalline elastic tensors provided by mineral physics studies. Our modelling predicts fast shear wave splitting directions for SKS, SKKS and ScS phases, as well as polarities of P- and S-wave reflections off the D" interface, for a range of propagation directions, via solution of the Christoffel equation. We run tests using randomly selected synthetic data sets based on a given starting model, controlling the total number of measurements, the azimuthal distribution, and the type of seismic phases. For each synthetic data set, we search over all possible elastic tensors and orientations to determine which are consistent with the synthetic data. Overall, we find it difficult to uniquely constrain the mechanism for anisotropy with a typical number of seismic anisotropy measurements (based on currently available studies) with only one measurement technique (SKS, SKKS, ScS or reflection polarities). However, data sets that include SKS, SKKS and ScS measurements or a combination of shear wave splitting and reflection polarity measurements increase the probability of uniquely constraining the starting model and its orientation. Based on these findings, we identify specific regions (i.e. North America, northwestern Pacific and Australia) of the lowermost mantle with sufficient ray path coverage for a combination of measurement techniques.

1 INTRODUCTION

Mantle convection finds its surface expression in plate tectonics and represents a crucial dynamic process in the deep Earth. Despite its importance, the pattern of mantle convection and the forces that drive mantle flow remain imperfectly understood. This is particularly true for the deepest mantle: flow at the base of the mantle likely influences (and/or is influenced by) structures such as large low shear velocity provinces (LLSVPs). Subducting slabs likely penetrate into the lower mantle and hot mantle plumes generate from or near the LLSVPs, indicating a strong connection between the surface and deep mantle processes (e.g. Garnero et al.2016).

Observations of seismic anisotropy have the potential to illuminate mantle flow, due to the relationship between strain due to mantle convection and seismic anisotropy via lattice preferred orientation (LPO) or shape preferred orientation (SPO) mechanisms. The presence of anisotropy in the D" layer at the base of the mantle has been known for several decades (e.g. Lay & Helmberger 1983) from the analysis of body wave phases (as summarized in Nowacki et al.2011). At this point a relatively small fraction (Fig. 1) of the core mantle boundary region has been explored for D" anisotropy using body waves. Fig. 1 shows a map, updated from Nowacki et al. (2011) illustrating the geographical coverage of previous studies (including recent work by Cottaar & Romanowicz 2013; Lynner & Long 2014a; Ford et al.2015; Long & Lynner 2015; Simmons et al.2015; Creasy et al.2017; Deng et al.2017; Thomas et al.20112011). Despite these observations, however, we still do not fully understand the anisotropy in these regions. Several different models for D’ anisotropy have been proposed, including those that invoke LPO of bridgmanite (Br), post-perovskite (Ppv), or ferropericlase (Fp) and those that invoke SPO of partial melt (see Nowacki et al.2011 for a review). The mechanisms responsible for D" anisotropy, the dominant slip systems involved, the orientation of the anisotropic fabric, and the implications for mantle flow geometries thus remain poorly understood.

A variety of body waves has been used to study anisotropy in the deepest mantle. Specifically, direct S, ScS and Sdiff have been used to observe lowermost mantle anisotropy by measuring shear wave splitting (e.g. Wookey et al.2005a; Cottaar & Romanowicz 2013; Ford et al.2006; Thomas et al.2007). Combinations of phases, such as SKS-SKKS (e.g. Wang & Wen 2007; Long 2009) or S-ScS (e.g. Wookey et al.2005a; Nowacki et al.2010), are often useful to isolate the lowermost mantle contribution to splitting. Thomas et al. (2011) used an array analysis technique to observe reflected P and S waves off the D" discontinuity; the azimuthal dependence of the polarity of D" reflections SdS and PdP contains information about lowermost mantle anisotropy. While body wave observations have been extensively used to study anisotropy at the base of the mantle, such studies suffer from the fundamental limitation of small azimuthal coverage; most studies are essentially restricted to a single ray path, which means that the geometry of anisotropy cannot be tightly constrained.

Several recent studies of deep mantle anisotropy have ameliorated this limitation by targeting regions of D" that are sampled by body waves over multiple azimuths (pink regions in Fig. 1). These include studies of the lowermost mantle beneath Siberia (Wookey & Kendall 2008; Thomas et al.2011), North America (Nowacki et al.2010), the Afar region of Africa (Ford et al.2015) and Australia and New Zealand (Creasy et al.2017). In some cases, one can test whether the observations could clearly distinguish among different mechanisms for anisotropy. For example, Ford et al. (2015) and Creasy et al. (2017) carried out forward modelling of ScS, SKS and SKKS splitting data sets over multiple azimuths to test whether a unique mechanism for anisotropy and/or a unique orientation of an assumed mechanism could be identified. In each of these studies it was found that LPO of Ppv matches the observations, but other mechanisms (such as LPO of Br or Fp) were also consistent with the data. None of the studies summarized in Fig. 1 has successfully identified a uniquely constrained mechanism or orientation for anisotropy. Motivated by this, we attempt here to understand what observations are needed to distinguish the various possible models for D’ anisotropy.

The goal of this study is to understand what combination of body wave data sets (SKS, SKKS, ScS and reflection polarities) are necessary to uniquely constrain the mechanism and geometry of anisotropy in the lowermost mantle using observations of shear wave splitting and D" reflection polarities. Such an understanding will aid in the design of future observational studies to maximize the chances of uniquely constraining a mechanism. We are interested in understanding the characteristics of data sets that are best suited to constrain the details of D" anisotropy, including the number of measurements needed, the optimal azimuthal coverage, and the optimal combinations of body wave phases. We address two specific questions: (1) What types of data sets (potentially including SKS, SKKS and/or ScS splitting, and/or reflection polarities) are needed to uniquely identify the causative mechanism for anisotropy (e.g. LPO of Ppv, Br, Fp or SPO of partial melt) and (2) if we assume that the mechanism for anisotropy is known to be LPO of Ppv, what type of data sets are needed to uniquely constrain the orientation of the anisotropy?

We carry out forward modelling tests for a suite of synthetic body wave data. Our approach to forward modelling of synthetic data sets follows our previous work on observations of shear wave splitting in D’ (Ford et al.2015; Creasy et al.2017) and also incorporates measurements of D" polarities of P- and S-wave reflections (Thomas et al.2011). Our approach is to test a variety of candidate elastic tensors that describe various mechanisms for lowermost mantle anisotropy. For each model, we randomly generate more than 5000 unique synthetic data sets (for SKS, SKKS and ScS shear wave splitting, plus PdP and SdS polarities) with a certain set of characteristics (e.g. number and type of measurements, as described below) and a random azimuthal distribution. For each set of random ray paths, we compute a set of predicted ‘observations’ of shear wave splitting and/or reflection polarities using a ray theoretical approach. We then attempt to determine what characteristics of body wave data sets are optimal for uniquely constraining anisotropy in the lowermost mantle.

2 METHODS

2.1 Candidate models for D’ anisotropy

We first consider which plausible models for D’ anisotropy should be tested. The lower mantle is likely composed of pyrolite (e.g. Lee et al.2004), a model composition that consists of ∼76 mol per cent of bridgmanite (Br: MgSiO₃), ∼17 mol per cent of periclase (Fp: (Mg,Fe)O), and ∼7 mol per cent of calcium perovskite (Capv: CaSiO₃). In the D" layer at the base of the mantle, we expect a phase change of Br to post-perovskite (Ppv: MgSiO₃) (e.g. Murakami et al.2004). Based on ab initio calculations and laboratory experiments, Br, Fp and Ppv all have strong single-crystal anisotropy, with Fp being the weakest mineral and the most anisotropic (as summarized in Nowacki et al.2011), although it is less abundant than Br/Ppv. This suggests that LPO development in any of the dominant lowermost mantle minerals may contribute to the observed anisotropy, as long as deformation is taking place in the dislocation creep regime (e.g. McNamara et al.2001). Another possible mechanism is aligned pockets of an elastically distinct material such as partial melt in configurations such as disks, tubes or sheets, creating shape preferred orientation (SPO, e.g. Kendall & Silver 1998, Tables 1 and 2).

We test a suite of models that describe single-crystal elasticity of lowermost mantle materials derived from ab initio calculations, following our previous modelling work (Ford et al.2015; Creasy et al.2017). This approach assumes that an aggregate will have the same anisotropic geometry (although not strength) as a single crystal. In addition to the single-crystal models, we test one model (for Fp LPO) based on deformation experiments (Long et al.2006) and models that invoke the SPO (shape-preferred orientation) of partial melt (Table 1), with elastic constants calculated using an implementation of effective medium theory within the MSAT toolbox (Walker & Wookey 2012).

Finally, our last candidate model approximates a textured Ppv aggregate and is derived from a 3-D, global mantle flow field calculation in combination with a visco-plastic self-consistent model LPO development in Ppv (Walker et al.2011). We determined a representative elastic tensor for Ppv texture development in high-strain simple shear by querying the TX2008.V1.P010 model of Walker et al. (2011), which combined a lower mantle viscosity model from Mitrovica & Forte (2004) with a mantle density model from Simmons et al. (2009). We only considered the case in which slip on the (010) plane dominates; this is the most likely slip plane for Ppv based on experiments (Yamazaki et al.2006; Walte et al.2009), modelling (Goryaeva et al.2017), and observations of D" anisotropy (Ford et al.2015; Creasy et al.2017; Thomas et al. 2011). To obtain a representative average tensor for simple shear, we identified a 15° by 15° geographical region of the global flow (beneath the northern Atlantic Ocean) that is dominated by strong horizontal shear. We then extracted and averaged the 16 elastic tensors (the model calculated tensors every 5°) from the resulting TX2008.V1.P010 elasticity predictions in this region.

2.2 Computation of reflection polarities and fast splitting directions

Given the full suite of candidate models for elasticity in D’ to be used in our study (Table 1), we implement methods for predicting various types of body wave observations for these scenarios. We calculated predicted shear wave splitting fast directions for SKS, SKKS and ScS phases (Fig. 2) over a range of azimuths (every 5°) and inclinations for each of these models (Table 1) by solving the Christoffel equation using the MSAT toolkit of Walker & Wookey (2012). The three different phases propagate at different inclination angles: ∼55°, 35°, 0° from the horizontal, respectively.

We then calculated the reflection polarities of SdS and PdP and the corresponding predicted shear wave splitting fast directions (Fig. 3) over a range of azimuths (every 5°) and inclinations for each of these models (Tables 1 and 2). Table 2 summarizes the models used to generate predictions of D" reflection polarities (SdS and PdP), including the assumed slip system, based on the methodology of Thomas et al. (2011). These models were constructed by assuming horizontal simple shear at the base of the mantle, where the dominant slip direction aligns parallel to the CMB, the slip plane is assumed to be horizontal, and 12 per cent of the aligned single crystals are mixed linearly with its isotropic equivalent. This choice of 12 per cent alignment was based on the previous work of Thomas et al. (2011), and yields reasonable anisotropic strengths; since we focus on reflection polarities and not amplitudes, however, this choice of value is not critical. We assume that the aligned grains are sub-parallel with the slip direction and the slip plane is subparallel to the CMB and the remaining grains are randomly oriented for Models A, B and C (Fig. 3). We tested three models (Models A [Ppv], C [Br] and D [Fp] in Table 2) in which the D" discontinuity represents a change in alignment of the mineral grains from an isotropic (above the discontinuity) to an anisotropic (below the discontinuity) regime. In Model B, the D’ discontinuity is an isotropic phase transformation from Br to anisotropic Ppv. The predicted values for reflection polarities for each model are shown in Fig. 3 and were calculated using Guest & Kendall (1993) from the velocity perturbation and reflection coefficients at the interface between an isotropic and anisotropic layer with respect to azimuth from the dominant slip direction and epicentral distance (Thomas et al.2011).

Our approach to calculating predicted shear wave splitting parameters and reflection polarities for our synthetic models makes several simplifying assumptions. First, we only directly model shear wave splitting due to lowermost mantle anisotropy, and ignore any potential contributions from the upper mantle. Our approach therefore assumes that any upper mantle contribution (in real data) has been correctly accounted for; we further assume that the bulk of the lower mantle is isotropic (Meade et al.1995). We do not explicitly consider how incorrect upper mantle corrections could bias the resulting D’ observations, which is beyond the scope of our study. Secondly, we rely on ray theory and do not consider finite frequency wave effects in our modelling. Ray theoretical predictions are generally adequate for homogenous regions of D’ (e.g. Nowacki & Wookey 2016), although they may break down for laterally heterogeneous anisotropic models (heterogeneities varying over hundreds of km). Third, in our modelling we approximate the propagation directions for SKS and SKKS with average inclination angles for these phases, and for ScS we assume that propagation is horizontal through the D’ layer. This assumption follows previous work (Nowacki et al.2010; Ford et al.2015; Creasy et al.2017). In the Earth and at the relevant epicentral distances, ScS can be inclined from the horizontal up to ∼15°, but this assumption has only a modest effect on the predicted splitting parameters. We assume the three different phases (SKS, SKKS and ScS) propagate at different inclination angles: ∼55°, 35°, 0° from the horizontal, respectively. Inclination angles are based on a straight-line approximation, calculated using TauP (Crotwell et al.1999) based on the PREM earth model (Dziewonski & Anderson 1981) for distances of 90°–120° for SKS/SKKS and 60°–80° for ScS with an event at a depth 10 km. We use these average propagation angles for SKS and SKKS in our modelling for simplicity, although for real data they can vary by 10° to 20° from these average values.

2.3 Modelling approach and strategy

Our goal is to conduct a series of stochastic forward modelling simulations to test whether we can uniquely constrain a given starting model (an elastic tensor) and its orientation using a data set with a given set of characteristics (e.g. number and type of measurements, azimuthal distribution). Our forward modelling framework follows, who modelled a shear wave splitting data set that samples the lowermost mantle beneath the Afar peninsula along the edge of the African LLSVP. We did not consider delay times in our modelling. Individual delay time measurements contain larger error bars, which limit the utility of using the relative traveltimes in a data set as a discriminant. The complete trade-off between fabric strength and layer thickness also limits the utility of using absolute traveltimes as a constraint.

For each of our modelling experiments, we first choose a starting model and orientation from the possibilities listed in Table 1. As an example, we first consider a horizontally aligned elastic tensor of Ppv with [100] and [010] axes parallel to the CMB, which we will use to illustrate our approach in several of the following figures. Secondly, we randomly identify a set of ray paths of SKS, SKKS and/or ScS sometimes in combination with SdS and PdP reflection polarities. Thirdly, we calculate the predicted fast-axis directions and/or reflection polarities of SdS and PdP for each ray path, as described in Section 2.2.

In the fourth step, we model the synthetic data set by applying the same forward modelling technique that we typically use for real data (Ford et al.2015). Specifically, we treat the synthetic observations as though the actual model used to generate them was not known, and test all possible models listed in Tables 1 and 2 in all possible orientations (every 5°) to identify models/orientations that are consistent with the synthetic data set. A candidate model/orientation is discarded if the predicted and ‘observed’ fast splitting directions differ by more than 20° or if the predicted reflection polarities are opposite to those of the ‘observations’. We apply this 20° cut-off for the splitting observations, based on methods and reasonable estimates of errors in previous shear wave splitting studies (see Ford et al.2015). While this misfit criterion is appropriate for measurement errors, it does not take into account effects such as inaccurate upper mantle corrections for actual D’ anisotropy observations or the possible finite-frequency effects of complex structure. Explicit consideration of these effects in D’ anisotropy studies is a subject of ongoing research. For each candidate model/orientation that was considered an acceptable fit to the synthetic data, we calculated a total misfit value for the fast polarization directions only based on a residual sum of squares approach, following Ford et al. (2015). Each fast direction misfit is normalized by the maximum residual of 90° and summed by using the residual sum of squares, in which we calculate the square of the difference between the observation and data prediction.

The fifth and final step in our modelling strategy is to repeat the entire process a large number (M) of times for random ray path configurations. All of these steps are illustrated in Fig. 4. In each iteration, we randomly choose a new azimuthal distribution of ray paths for a new synthetic data set with varying characteristics (such as the number and type of observations, described in more detail below). We report our results by considering what percentage of the M iterations could uniquely identify the starting model. Each individual iteration was designated as ‘uniquely constrained’ if it successfully identified the correct starting model and could completely rule out any other candidate model. However, if there was at least one other anisotropy configuration (any candidate elastic tensor model, in any orientation) which was found to be consistent with the synthetic observations, then that iteration was designated ‘not uniquely constrained.’ Therefore, all our model results are characterized through a per cent-uniquely constrained value, which identifies what percentage of the M simulations could uniquely constrain the starting model. The actual values of these ‘ per cent-uniquely constrained’ estimates are strongly dependent on our modelling choices, and the estimates could change with different modelling assumptions. However, these percentages can be compared across our suite of numerical simulations, since our assumptions are consistent across the various tests.

We tested different combinations of N and SKS number to gain insight into how many measurements, and in what combination, are typically needed to uniquely constrain the anisotropy. For angular dispersion, we calculated the value of R for each of the M iterations carried out in each test; then, we queried the large number of simulations to understand how the azimuthal distribution of the synthetic data affected its ability to constrain anisotropy.

2.4 Distinguishing the mechanism and orientation of anisotropy

For the first round of tests, we sought to understand how many shear wave splitting measurements, and in what combination (as described by the SKS number), are generally needed to uniquely constrain the mechanism for anisotropy. That is, we tested whether synthetic data sets could be shown to be consistent only with the correct starting model (e.g. Ppv, as opposed to other models listed in Table 1), and with no other candidate mechanism. For this round of tests, we used the single crystal models in Table 1 as starting models, each in several different orientations. The LPO model of Ppv was only used for the second round of tests. We defined the starting model orientation via the rotation angle about the [100] axis from the horizontal (note that Fig. 2 only shows an example with a horizontal [100] and [010] direction). We arbitrarily tested each single crystal model at three different orientations based on the rotation angle about the [100] axis from the horizontal: 0°, 45°, and 90°. For the LPO model of Ppv, we only test the original orientation for the starting model and do not test a rotated version of the elastic tensor since this model is based on a region in the Walker et al. (2011) with horizontal shear. In our initial round of tests, we focused only on cases in which shear wave splitting observations of SKS, SKKS, and ScS (for varying N, SKS number, and R values) were used to constrain the models. In later tests, we explored scenarios in which reflection measurements were combined with shear wave splitting data in order to estimate the improvement obtained by combining different data types.

We also carried out a series of tests whose goal was to constrain the orientation of the elastic tensor for the case in which the mechanism for anisotropy is known (or assumed). For this line of inquiry, we focused on Ppv as a test case; we did not test other mechanisms in this part of the study. The choice to focus on Ppv was made for simplicity and because Ppv is often invoked as the preferred mechanisms for anisotropy in D" (Wookey et al.2005b; Nowacki et al.2010; Thomas et al.2011a; Ford & Long 2015; Ford et al.2015; Creasy et al.2017). We consider both single-crystal Ppv tensors and elastic tensors derived from texture modelling, as discussed above. As in our first series of tests, we initially focus on synthetic data sets that only contain shear wave splitting observations, and then examine cases that also include reflection measurements.

Lastly, in addition to the two major lines of inquiry we address in our modelling (what kind of data sets are needed to constrain the mechanism and orientation) of lowermost mantle anisotropy, we performed two practical tests using horizontal Ppv as a starting model. First, we carried out a test of how many iterations (that is, values of M) are needed for our forward algorithm to converge on an estimate of the probability of identifying unique models. Second, we tested the addition of Gaussian noise to the shear wave splitting predictions, in order to understand how well real, noisy data sets might perform. The results of these practical tests are described below. While seismic data can deviate from a Gaussian distribution (Groos & Ritter 2009), we only consider Gaussian distributed noise here since in an ideal case, seismic noise is Gaussian distributed (Bendat & Piersol 2011).

3 RESULTS

3.1 Illustrative examples: Model runs for a ppv starting model

To illustrate the process and results of our modelling, we discuss here the results from a test that attempts to constrain the starting model, as well as one iteration of a test that attempts to constrain the orientation. In both cases, we use synthetic shear wave splitting data only. For these examples, as in all of our tests, we follow the five steps of our method outlined above (Fig. 4): (1) choose a starting model and orientation, (2) choose the number of observations and the SKS number to randomly generate a distribution of ray paths, (3) calculate the predicted fast polarization directions (of SKS, SKKS, and ScS) and reflection polarities (for SdS and PdP) for the synthetic data set for the chosen starting model, (4) conduct a forward modelling search over all possible orientations for all possible candidate models to eliminate all models/orientations that do not fit the ‘observations’ using a misfit cut-off. Then, if all other models and orientations can be eliminated by applying the misfit cut-off, this set of synthetic ray paths are able to uniquely constrain the starting model and designated as ‘uniquely constrained.’ The fifth step would be repeating this process M number of times but for this illustrative example, M = 1.

Our illustrative example is shown in Fig. 5. For this example, we chose a starting model of non-rotated Ppv (in this case, the [100] and [010] crystallographic axes are parallel to the CMB) (Fig. 5a). In all of our single-crystal elasticity tests, we do not assume a dominant slip system; rather, we invoke a starting orientation in the geographic reference frame identified by the angle of the mineralogical axes. This particular example involves 9 splitting observations, 6 of which are SK(K)S (that is, an SKS number of 2/3). The randomly generated azimuthal distribution of these chosen phases is shown in Fig. 5(b). The predicted fast polarization directions for our chosen model and ray configuration, plotted in a ray-centered reference frame, are shown in Fig. 5(c). A search over all possible candidate models and orientations (rotating every 5°) shows that there is no other model, other than the correct starting model (Ppv), that can match each of the synthetic fast splitting directions to within 20° (our pre-defined misfit cut-off). Put another way, for every possible combination of starting model and orientation (other than the correct, known starting model), at least one predicted fast splitting orientation differed from that in the data set by more than the 20° misfit cut-off. Since this particular configuration of observations could uniquely identify the starting model and no other models, it is designated ‘uniquely constrained.’

This particular example illustrates a single iteration (M = 1) of our testing, but the power of our approach lies in repeating this a large number of times to understand what percentage of randomly generated synthetic data sets have the ability to uniquely constrain the starting model. In order to understand how many iterations are needed to converge on an estimate of this probability, we conducted an ‘iteration test’ for our horizontal Ppv starting model, as shown in Fig. 6. For this test, we used 9 shear wave splitting measurements (N = 9) and an SKS number of 2/3, as in the example shown in Fig. 5, and ran a large number of iterations (M = 50,000), each involving a new, random distribution of propagation azimuths. After each successive iteration, we calculated the percentage (of M iterations) for which the synthetic data set was able to uniquely constrain the starting model, as shown in Fig. 6. For this starting model, after a large number of iterations, we found that 41 per cent of all iterations could uniquely constrain the starting model, while for the other 59 per cent of the ray path configurations, there was another model/orientation that could simulate the synthetic data. Our running estimate of how likely a data set with nine splitting observations (6 SK(K)S, 3 ScS) converges on an average value of 41 per cent after approximately 1000 iterations (Fig. 5). Based on this iteration test, we have chosen to run each of our numerical experiments for M = 5000 iterations, balancing computational cost and the need for our estimates to converge.

Next, to illustrate our process for testing whether synthetic data can identify a unique starting orientation, we show in Fig. 6 two examples of searching for the correct starting orientation for the same horizontal Ppv starting model as in Fig. 4. For this example, we chose two different ray path configurations, one with N = 8 observations (5 SKS + SKKS and 3 ScS; Fig. 7a) and one with N = 4 (3 SKS + SKKS and 1 ScS; Fig. 7b). We assume that the mechanism for anisotropy is known to be Ppv and that the elastic constants are known, and search over all possible orientations to test whether there are additional configurations (other than the known starting orientation) that can reproduce the synthetic observations.

Figs 7(a) and (b) show all possible orientations that satisfy this suitability criterion for each of our two examples (N = 8 and 4, respectively), with each orientation colour-coded by its calculated misfit value (eq. 1). Following Ford et al. (2015), we search for local minima of misfit within the 3-D rotation space. For our N = 8 case (Fig. 7a), the set of 8 measurements could uniquely identify the starting orientation, and would be designated as ‘uniquely constrained.’ However, for our N = 4 case (Fig. 7b), we identified two other possible orientations (that is, the known correct starting orientation, plus two others). Therefore, for this particular ray path configuration, the solution is designated ‘not uniquely constrained.’ We note, however, that the orientation with the lowest misfit value (magenta dot in Fig. 7b) is, in fact, the correct starting orientation.

Finally, we illustrate an example calculation that includes Gaussian noise in the synthetic observations (Fig. 7c). This test relies on the same horizontal Ppv starting model, and uses the same ray path configuration (N = 4) as the test shown in Fig. 7(b). The only difference is that when the predicted shear wave splitting fast directions are calculated based on the starting model and ray path distribution, we add Gaussian noise to the fast splitting direction ‘observations,’ with a maximum error excursion of 20° and a standard deviation of 9°. Fig. 7(c) reveals that the case with Gaussian noise produced the same two possible sets orientations as fitting the data, but now the solution with the minimum misfit is not associated with the correct solution.

3.2 Results: constraining the anisotropy mechanism

Building on the illustrative examples discussed in Section 3.1, we now explore the results of a large number of simulations with different starting models and ray path configurations. We first address the question of what kind of data sets are needed to distinguish among the various models listed in Figs 2 and 3. For this suite of numerical experiments, we examined a variety of starting models and orientations, as well as a variety of ray path configurations (as defined by the number of splitting measurements, the SKS number, and the angular dispersion of the ray path azimuths). The results of these experiments are shown in Fig. 8. We first examine those model runs that only included shear wave splitting data, shown in the nine panels of Fig. 8(a).

We initially focus on the mechanism and orientation of the starting model (Fig. 8a, left-hand panels), and explore how the probability of uniquely constraining the mechanism varies as a function of the number of measurements. For each of the models considered, the probability of identifying the unique starting model increases with the number of measurements as expected. Typically, there is a sharp increase in the probability for N values between six and nine measurements. In all cases, approximately nine measurements are needed in order to have a ∼50 per cent chance of constraining the starting model, while a high number of splitting measurements (N ≈ 15) is needed for the probability to reach ∼90 per cent. For comparison, the data sets of Ford et al. (2015) and Creasy et al. (2017) contained between four and eight splitting measurements. The starting model with the highest success rate at constraining the mechanism is Br, as opposed to Fp and Ppv.

The probability of constraining the starting mechanism depends on the orientation of the starting model; as shown in Fig. 8(a), we tested orientations with a horizontal [100] crystallographic axis, 45° rotated about the [100], and 90° rotated about the [100] axis. Interestingly, for Ppv it is easier to uniquely constrain the starting model in the 90° case; in contrast, for Br the chances are highest for the horizontal case, and for Fp the chances are substantially higher for the tilted case. The reason for this result for Ppv can be discerned by examining the predicted splitting patterns in Fig. 2. For the horizontal case, predicted fast splitting directions for ScS do not vary with azimuth; however, if the Ppv tensor is rotated by 90° about the [100] axis, there is significant variation in fast directions with azimuth. With greater variability in the predicted fast polarization directions (lower angular dispersion), there is a higher probability of constraining that model for a given number of ScS observations. A similar principle is at work for Fp: ScS fast directions do not vary with azimuth for either horizontally or vertically aligned Fp, but in the tilted case, variability is present. Generally, the anisotropy scenarios that yield higher chances of uniquely constraining the starting model have lower mean angular dispersion values of the predicted fast-axis directions (Fig. S2). Models that have little variation in fast-axis directions with azimuth, such as non-rotated Fp, are more difficult to uniquely constrain (Figs S2c and 8a).

We also examined how the balance between SKS + SKKS versus ScS phases in the synthetic data set affected the ability of the synthetic ‘observations’ to uniquely constrain the starting model (Fig. 8a, middle panels). For these experiments, we varied the SKS number from 0 (all ScS measurements) to 1 (all SKS + SKKS measurements) for a fixed value of N = 9. For very high or low values of SKS number we find a low probability of uniquely constraining the starting model with substantially higher probabilities for intermediate SKS numbers. The optimal ratio of SK(K)S phases to total measurements differs slightly for different starting models, but in general an SKS number between 0.5 and 0.8 maximizes the chances of constraining the anisotropic mechanism. In all cases, a combination of ScS and SK(K)S shear wave splitting observations, instead of splitting measurements for just one phase type, will drastically improve the probability of constraining the starting model.

Additionally, we explored the importance of how the angular distribution of the synthetic ray paths affected the ability to constrain the starting model, finding only a weak effect (Fig. 8a, right-hand panels). As expected, data sets with a wide angular distribution (R < 0.2) have the largest probability of uniquely constraining the starting model in all cases. At very large values of angular dispersion (R > 0.8), for which the ray paths are clustered over a narrow range of azimuths, the splitting ‘observations’ are sampling similar parts of the elastic tensor. Because of this, data sets that are tightly clustered in azimuth cannot capture the symmetry of the tensor and cannot distinguish among different candidate mechanisms for anisotropy. For intermediate values of R, the dependence on R is not strong.

Finally, we explored the value of combining shear wave splitting and reflection polarity measurements when trying to uniquely constrain an anisotropic model. Fig. 8(b) shows the results of adding a single reflection polarity measurement (that is, a measurement of PdP and SdS polarities for a single ray path) to a data set of shear wave splitting measurements. For this test, we considered a smaller number (four) of potential candidate models (as shown in Fig. 3), so the probabilities of uniquely constraining the anisotropy mechanism are generally higher than in our other tests. For this test, we chose a ray path configuration involving an SKS number of 0.67 and varied the number of shear wave splitting measurements from 0 to 15. We used a starting model A in Fig. 3 (anisotropy due to Ppv), and tested configurations that involved both shear wave splitting measurements and one additional set of reflection polarity measurements (both PdP and SdS) at a single azimuth. This test (Fig. 8b) demonstrates that despite the fact that reflectivity measurements provide only binary information (positive or negative polarities), the incorporation of a different data type into the test increases the probability of uniquely constraining the starting model. In some cases, this increase is substantial; specifically, for data sets containing between four and eight shear wave splitting measurements. The addition of reflection polarity data can increase the probability of constraining the starting model by ∼10–18 per cent (right panel of Fig. 8b).

3.3 Results: constraining the anisotropy orientation

The tests shown in Fig. 8 illustrate the ability of shear wave splitting and reflection polarity data to constrain the anisotropic mechanism if the algorithm is allowed to consider a range of possible models. We now turn our attention to tests in which we assume that the mechanism that creates the anisotropy, as well as the elastic constants associated with that mechanism, are known, but the orientation of the elastic tensor is not known. In general, this is an easier problem than uniquely constraining the starting model, as the observations need not distinguish among different candidate elastic tensors, only among different possible orientations. In practical terms, this type of modelling exercise would be suitable for data sets that sample a region of the lowermost mantle whose mineralogy and temperature conditions can be constrained using independent observations or models (for example, seismic velocities in combination with mineral physics constraints).

For this set of tests, we first consider single-crystal Ppv in three different configurations: (1) [100] and [010] axes oriented in the horizontal plane, (2) a 90° rotation about the [100] axis and (3) randomly chosen orientations. For the third configuration, we randomly identified nine different, unique starting orientations. These randomly generated orientations were used for each of the ∼5000 iterations in this scenario. As with the tests discussed in Section 3.2, we tested a variety of ray path configurations with a range of N (number of measurements), SKS number and examined how our results varied with the angular dispersion characteristics of the synthetic ray paths. The results of our single-crystal Ppv tests are shown in the top row of Fig. 9. The results for our collection of nine random starting orientations are shown in detail in Fig. S3.

As expected, our tests demonstrate that uniquely constraining the orientation of the starting model is much easier and requires fewer measurements than uniquely constraining the starting model/mechanism (Fig. 7). In general, a ∼50 per cent probability of correctly retrieving the anisotropy is achieved with as few as six to nine splitting measurements (top left panel of Fig. 9a). The orientation of the starting model does affect the likelihood of uniquely identifying the anisotropy orientation. With our randomly generated starting orientations, the probability of constraining the starting orientation varies (Fig. S3), but on average randomly oriented starting models do slightly worse compared to the results shown in Fig. 9. For nine measurements, the randomly orientated models on average find the correct orientation in 65 per cent of all simulations, compared to Fig. 9, where a non-rotated and Ppv rotated by 90° can constrain on average 75 per cent of the simulations. As with our previous tests, it is clear that a mixture of SK(K)S and ScS shear wave splitting measurements provide the highest likelihood of constraining the starting orientation, although the optimal mix of ScS and SK(K)S depends on the starting model orientation. Our tests confirm that data sets that contain only ScS measurements (that is, SKS number of zero) cannot constrain the azimuth of the Ppv elastic tensor if its [100] axis is horizontal, due to the lack of variability in predicted fast polarization direction (Fig. 2). The dependence of our results on angular dispersion of the propagation azimuths (right-hand panels of Fig. 9a) are similar to those for the case in which we attempted to retrieve the starting model; in general, a wide distribution of azimuths will increase the probability of uniquely constraining the orientation of Ppv, while data sets whose propagation azimuths are tightly clustered are less ideal. The same is generally true for the random starting models, despite some small excursions from the overall trend (Fig. S3). These small excursions or ‘bumps’ in the curves are artifacts, and are related to stochastic variations in the distribution of the predicted fast splitting directions for different models.

Next, we considered elasticity models that explicitly take into account texture development in a polycrystalline aggregate, in addition to the single-crystal elastic tensors that are the main focus of our study. While there are many uncertainties in texture models for Ppv at lowermost mantle conditions, these models may be more representative of a realistic texture of aligned Ppv mineral grains. We only considered one case, invoking dominant slip on the (010) plane. Somewhat surprisingly, we found that for modelled Ppv LPO, there is a much lower probability of constraining the orientation of the elastic tensor than for test cases that used a single crystal elastic tensor (Fig. 8a). We investigated possible reasons for this, and found that in contrast to the single-crystal models, for the textured Ppv model it is fairly common for the algorithm to identify what we term as ‘unstable’ solutions, which are illustrated in Fig. S4. In this situation, a certain orientation might fit the observations, but adjacent orientations (in which the elastic tensor is rotated by 5°) do not. This is in contrast to the behaviour of single-crystal elastic models (Fig. 2), in which the best-fitting orientations are adjacent to other solutions that also fit the data (in other words, the misfit values vary smoothly as a function of rotation angles of the candidate tensors). In addition, the presence of unstable solutions is highly dependent on our use of the misfit criterion of 20°. Fig. S4 shows results for a range of misfit cut-off values, and demonstrates that these unstable solutions disappear with the application of more conservative misfit criteria.

We define a ‘stable’ solution as one in which, if the elastic tensor is rotated slightly (∼5° in any direction), the rotated elastic tensor would still yield an acceptable fit to the synthetic data. In contrast, an ‘unstable’ solution is one that has no adjacent orientations that yield an acceptable fit to the data. For the case of the textured Ppv model, the algorithm generally identifies many ‘unstable’ orientations (Fig. 9); again, this is in contrast to the generally ‘stable’ orientations identified for single-crystal Ppv (Fig. 7). In order to illustrate the effects of these unstable solutions, we applied a sensitivity cut-off to our textured Ppv simulations (Fig. 8a, second row) to illustrate the effects of removing all unstable solutions. If we consider only stable solutions, the probability of uniquely constraining the starting orientation increases by 20 per cent on average (Fig. 9a).

In order to identify the starting orientation of the Ppv LPO, we found that a mixture of SK(K)S and ScS shear wave splitting measurements again provide the highest likelihood of constraining the orientation (Fig. 8a: bottom, middle panel). There is a clear dependence on angular dispersion (Fig. 9a: bottom, right-hand panel). Specifically, with low values of R (0–0.1) and middle values of R (0.5–0.7), there is a higher probability of constraining the orientation, while there is a decrease in probability between R = 0.1 and R = 0.4. In all other cases and starting models, we have not observed this pattern of dependence with R. While there is no explanation for this pattern, large values of R (0.8–1.0) resulting in low probabilities of finding the starting orientation is consistent with all other tests.

Returning to our consideration of single-crystal Ppv models, and as in Section 3.2, we considered the effect of adding a reflection measurement to shear wave splitting observations to constrain the orientation of the single-crystal Ppv starting model (Fig. 9b). For this test, we used a starting model that invokes an isotropic ppv layer over an anisotropic ppv layer with dominant [100](010) slip (Model A in Fig. 3). As in the previous test, we find that just adding one observation of reflection polarity measurements improves the probability of constraining the starting orientation (Fig. 9b), although the improvement was somewhat less dramatic. Again as with the previous tests, the relative improvement is greatest for data sets with number of measurements N roughly between 5 and 9.

Finally, in a test analogous to the Gaussian noise test discussed in Section 3.1 and illustrated in Fig. 7(c), we considered a single-crystal Ppv test in which we tried to retrieve the correct starting orientation using synthetic observations that included random, Gaussian distributed errors on the fast polarization predictions (Fig. 10). We found that adding Gaussian noise to the fast polarization directions, normally distributed between –20° and 20° with a mean of 0° and standard deviation of 9°, does not significantly hinder the probability of constraining the starting model's orientation (Fig. 10a). However, this test allowed us to explore the distinction between uniquely constraining the starting model's orientation and identifying a model with a minimum misfit value that corresponds to the correct starting orientation. For the error-free synthetic data sets, the minimum misfit value always corresponds to the correct orientation, even for cases in which other orientations are allowed by the data. For cases where Gaussian error is incorporated into the synthetic data set, the synthetic data set may result in an incorrect solution, where the minimum misfit value may be different from the correct solution. This observation led us to carry out a test (Fig. 10b) in which rather than attempting to uniquely constrain the correct starting orientation, we tested whether the best-fitting orientation (that is, the candidate orientation with the minimum misfit value) actually corresponded to the correct starting orientation. We further tested whether the best-fitting solution in terms of misfit was oriented within 10°–20° of the known, correct starting orientation. Encouragingly, we found that the probability that the minimum misfit solution was within 20° of the correct orientation exceeded 50 per cent for data sets with a relatively small number of shear wave splitting measurements (N ≈ 4).

4 DISCUSSION

4.1 Implications for the interpretation of real-world data sets

Understanding the scope of information about lowermost mantle anisotropy contained in shear wave splitting and reflection polarity observations is crucial for our ability to relate anisotropy observations to flow at the base of the mantle. While the mechanisms of lowermost mantle anisotropy remain imperfectly known, the results presented in this paper reveal observational strategies that can maximize the probability of constraining the mechanism and/or orientation, regardless of the actual anisotropic geometries present. This work shows that a diversity of shear wave splitting measurements and reflection polarity data is essential, and the modelling of single phases (e.g. ScS, SKS, SdS) is typically insufficient to constrain anisotropic geometry.

Specifically, this work demonstrates that because different seismic phases (ScS, SKS, SKKS, PdP and SdS) propagate through or reflect off the D" region at different angles from the horizontal, a combination of these phases is more useful for constraining anisotropy than data sets with wide azimuthal coverage. Consider, for example, a hypothetical case in which nine unique splitting measurements for ScS phases are used to probe an anisotropic structure consisting of horizontal, single crystal post-perovskite. In this case, post-perovskite can only be distinguished from other plausible anisotropic models less than 10 per cent of the time (Figs 8a and 9a). However, if SK(K)S phases and/or reflection polarities are incorporated into the analysis, then we can distinguish between the two possible mechanisms nearly 40 per cent of the time (Fig. 8a). In all cases of varying starting models and orientations, a combination of different types of data increases the probability of constraining the starting model by 10 per cent to 60 per cent. This pattern also generally holds true for finding the orientation of the Ppv elastic tensor. A diversity of data increases the likelihood of constraining the orientation of Ppv anisotropy anywhere from 10 to 50 per cent for 6 unique measurements. Interestingly, we observed an exception to this (Fig. 9a) for Ppv oriented at an azimuth of 90°, where only ScS splitting data (SKS number of 0) had the best chance to constrain the starting orientation.

Body wave data sets that probe seismic anisotropy in the lowermost mantle should combine both multiple data types and wide azimuthal coverage to maximize the probability that the anisotropic geometry can be tightly constrained. Fig. 11 illustrates regions in the mantle in which all of the body wave measurement methods could potentially be applied simultaneously. This map was generated by considering the actual distribution of high-magnitude (M > 6.5) seismicity on Earth, in combination with a database of long-running broadband seismic stations beneath which the upper mantle anisotropy pattern has been shown to be simple enough to correct for (Lynner & Long 2013, 2014b). While there are many regions of D’ with limited ray path coverage for the types of data considered in this study, we find that North America, the Arctic, northwestern Pacific and Australia are regions that represent ideal targets to collect a diverse set of observations to further constrain D’ anisotropy.

Our results inform our view of why previous studies that included crossing ray paths (e.g. Ford et al.2015; Creasy et al.2017) were unable to uniquely constrain a model for D’ anisotropy. Our study indicates that a relatively large number of shear wave splitting measurements (approximately 9 or more for most cases in Fig. 8a) are needed to have at least a 40–60 per cent chance of uniquely identifying the starting model. The observational data sets of Ford et al. (2015) and Creasy et al. (2017) included approximately four to eight shear wave splitting measurements over unique azimuths in the lowermost mantle (Table 3). The synthetic models presented in this paper help to provide context for why these studies have not been able to uniquely constrain a particular mechanism or orientation for anisotropy (e.g. Ford et al.2015; Creasy et al.2017). For example, each of these studies (Table 3) had relatively high angular dispersion values for their range of predicted fast splitting directions (greater than 0.4 in all cases). As discussed in Section 3.2, data sets with lower angular dispersion values are generally more successful at constraining a unique elastic tensor. Therefore, even though many of the studies listed in Table 3 used diverse data types with combinations of SKS, SKKS and ScS, they could not uniquely constrain an anisotropy mechanism or orientation when testing the elastic tensors considered in this study. The studies that used one type of observation (Nowacki et al.2010; Thomas et al.2011) did not consider all possible elastic tensors and orientations that we tested here, so we cannot directly compare them with the results of our synthetic tests. If the mechanisms for anisotropy and the associated elastic tensors can be reliably assumed, there is generally a higher chance of identifying the correct orientation and inferring the correct mantle flow geometry. With only nine measurements, there is a 40–80 per cent chance of uniquely constraining the orientation of post-perovskite (Fig. 9a), an improvement from the chance of uniquely identifying the elastic tensor itself (a 40–60 per cent chance). Consideration of these results in future studies of D’ anisotropy, as well as a more detailed consideration of how errors and uncertainties propagate in forward models, should enhance our ability to characterize anisotropy at the base of the mantle. A more detailed statistical analysis may be required similar to this study to fully explore the error and model space.

4.2 Practical considerations

Our tests that assumed Gaussian error on the predicted fast splitting directions (Fig. 10b) demonstrate that it does not significantly affect the probability of constraining the model, as compared to noise-free synthetic data. When including Gaussian error, we found that as few as four shear wave splitting measurements can identify the correct orientation within 20° more than 50 per cent of the time. Data sets of this size (roughly four unique measurements in the same region of D’) can likely be reasonably achieved in many regions of the lowermost mantle, based on the distribution of available ray paths (Fig. 11). This finding may help with the interpretation of modelling results for real splitting data sets, such as those considered by Ford et al. (2015) and Creasy et al. (2017), for which multiple possible anisotropic orientations were identified, but particular orientations had significantly lower misfit values than others.

Our synthetic modelling results also shed light on potential complications in interpretation caused by the different symmetry classes of some of the candidate elasticity scenarios that have been proposed to explain lowermost mantle anisotropy. To effectively differentiate these scenarios using shear wave splitting data alone, it is crucial for splitting data sets to probe the symmetry of the mineral such that no other elastic tensor simulates that pattern for a similar range of propagation directions. Of the candidate scenarios we tested in this study, Fp has the highest (cubic) symmetry with only three unique constants in the elastic tensor. SPO models have the next highest symmetry, since tubule and oblate SPO models are hexagonal (transversely isotropic) with five unique elastic constants. Ppv and Br are both orthorhombic, with the same order of symmetry and only nine unique elastic constants. In more complicated models, such as LPO calculations of single crystals, the symmetry is much lower than its single crystal counterpart with 21 unique elastic constants.

4.3 Limitations of our modelling approach

We caution that our synthetic tests must be interpreted in light of the still-considerable limitations in our understanding of the elasticity of anisotropic materials at lowermost mantle conditions. We have focused mainly on single-crystal elastic tensors, derived mainly from ab initio simulations, as a reasonable starting point in this study; however, predictions of single-crystal elasticity are likely imperfect and do not take into account effects such as variation in composition. Furthermore, single-crystal elasticity is an imperfect proxy for the likely anisotropic geometry of polycrystalline aggregates, particularly for minerals with high symmetry such as Fp (e.g. Yamazaki & Karato 2002). The further consideration of elasticity models that explicitly take into account texture development will be an important step, although texture models include a number of poorly known parameters (such as activation energies for different slip systems) and consensus on the dominant slip systems in different lowermost mantle minerals remains elusive (e.g. Nowacki et al.2011).

Another limitation of the work proposed here is that it is carried out in the context of ray theoretical predictions, assuming infinite frequency, rather than considering the full characteristics of the waveform at finite frequencies. With improvements on both observational and modelling techniques that model the full waveform (e.g. Kawai & Geller 2010; Nowacki & Wookey 2016; Parisi et al.2018), the interpretation of seismic anisotropy observations can very likely be improved. In particular, future work must investigate how the measurement techniques used influence the interpretation of finite frequency waveform effects and to what extent ray theoretical predictions are a useful approximation. Despite these limitations, we expect that future work that predicts body wave observations in the presence of lowermost mantle anisotropy in a finite-frequency framework will likely find similar results: a diversity of seismic phases and measurement yields the best probability of capturing the symmetry, orientation and properties of an elastic tensor. While this study is limited to a specific set of currently available elastic tensors from the mineral physics literature, our overall findings should be generally applicable and adaptable to future improvements of our knowledge of lowermost mantle elasticity.

5 SUMMARY

To summarize, the complete characterization and interpretation of seismic anisotropy at the base of the mantle would have profound effects on our understanding of lower mantle dynamics, potentially yielding insights into the pattern of mantle flow. Many recent studies have pointed to the difficulty of distinguishing different models of lowermost mantle anisotropy with body wave observations, given challenges with data coverage and uncertainties in the mechanism for anisotropy and the relationships between deformation and the resulting anisotropy at lower mantle conditions. In this study, we conducted a series of Monte Carlo simulations to determine what combination of body wave data sets (shear wave splitting and reflection polarities) are required to constrain D" anisotropy. We tested various starting models, orientations, and methods for the detection and identification of D’ anisotropy. The modelling approach in this study is applicable to a wide range of elasticity models, and can be extended as our knowledge of the physical properties of the lowermost mantle increases. This approach can be used in future work on D’ anisotropy to further explore how well a data set can discriminate among possible elastic tensors.

Our results show that a diversity of observational techniques, including different types of seismic phases propagating over a range of ray path directions, are necessary in order to maximize the chances of constraining anisotropy at the base of the mantle. A combination of shear wave splitting measurements and observations of PdP and SdS reflection polarities in the same regions may be particularly powerful. We have further shown that if the mineralogy and/or mechanism for anisotropy can be constrained from independent data, then the orientation of the elastic tensor (and thus information about patterns of mantle flow) can likely be retrieved from observational data sets that include a relatively modest number of measurements.

SUPPORTING INFORMATION

Figure S1. Schematic diagram showing the definition of angular dispersion (R), with arrows indicating a direction anywhere from 0° to 360°. Small values of R indicate a wide distribution of directions, while larger values indicate a tight configuration.

Figure S2. Angular dispersion—R—plots of all predicted fast axis directions for SKS (red), SKKS (orange), SKS and SKKS (violet) and ScS (blue). Eqs (2) and (3) from the main text were used to calculate R–; however, since fast axis directions can only vary from –90° to 90°, the fast axis orientations were adjusted from –90 to 90 to 0 to 360. Angular dispersion is plotted in terms of the clockwise rotation angle about the [100] axis of (a) post-perovskite (PPV), (b) bridgmanite (Br), (c) ferropericlase (Fp) and (d) LPO of PPV of the given models, respectively. Magenta triangles indicate the starting models’ orientations used in the modelling in the main text as in Figs 8 and 9.

Figure S3. Results of synthetic tests that aim to uniquely constrain the orientation of a Ppv starting model, as discussed in Section 3.3 in the main text, by selecting nine random the starting orientations of the Ppv tensor. Each figure illustrates the probability of uniquely identifying the given starting model's orientation using the synthetic data, for nine different orientations of Ppv. As in Fig. 7, each column represents the variable that was allowed to vary, while the other two were fixed. For (a), the SKS number was fixed (0.67). For (b), the number of measurements was fixed to six. For (c), the number of measurements was fixed to nine measurements and an SKS number of 0.67. The blue lines in (b) and (c) correlate to the blue line in (a), where the variation of SKS number and angular dispersion were only tested for one of the starting model orientations.

Figure S4. Example illustrating the identification of unstable solutions 49 in a test that aimed to identify the starting orientation of a horizontal textured Ppv LPO. The polar plots show all possible orientations that fit a given synthetic data set. The colours represent misfit values. The white circles mark the minimum misfit of each cluster of possible orientations. The magenta circles show the correct solution. The orange boxes highlight some of the unstable solutions. Each row represents a different misfit cut-off. The top row represents the cut-off used in this study (20°). The second row uses a cut-off of 15° and the third row 10°. With a lower cut-off, the unstable solutions are eliminated. The bottom elastic tensors on the left show and example of one of the resulting unstable solution from the orange boxes in the first row. The elastic tensor on the right shows the same tensor but rotated by 5°, which fails the misfit criterion. The magenta lines represent the measurements used in this simulation. Colours here represent %S-wave anisotropy.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

ACKNOWLEDGEMENTS

This work was supported by a National Science Foundation (NSF) Graduate Research Fellowship grant DGE-1122492 to N.C. and by NSF grant EAR-1547499 to M.D.L. Some figures were prepared using the Generic Mapping Tools (Wessel & Smith 1991). We thank the Yale Center for Research Computing for guidance and use of the research computing infrastructure, specifically Kaylea Nelson. AP received funding from the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No 642029–ITN CREEP. We are grateful to Sanne Cottaar, Jeroen Ritsema, and editor Ana Ferreira for thoughtful and constructive comments that helped us to improve the paper.

REFERENCES