Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Exact free oscillation spectra, splitting functions and the resolvability of Earth's density structure

Exact free oscillation spectra, splitting functions and the resolvability of Earth's density... Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 Geophysical Journal International Geophys. J. Int. (2018) 213, 58–76 doi: 10.1093/gji/ggx539 Advance Access publication 2017 December 18 GJI Seismology Exact free oscillation spectra, splitting functions and the resolvability of Earth’s density structure 1 2 1 1 3 F. Akbarashrafi, D. Al-Attar, A. Deuss, J. Trampert and A.P. Valentine Department of Earth Sciences, Universiteit Utrecht, Heidelberglaan 2, 3584CS Utrecht, The Netherlands Bullard Laboratories, Department of Earth Sciences, University of Cambridge, Cambridge CB3 0EZ,UK Research School of Earth Sciences, The Australian National University, Canberra, ACT 0200, Australia Accepted 2017 December 15. Received 2017 December 7; in original form 2017 August 23 SUMMARY Seismic free oscillations, or normal modes, provide a convenient tool to calculate low- frequency seismograms in heterogeneous Earth models. A procedure called ‘full mode cou- pling’ allows the seismic response of the Earth to be computed. However, in order to be theoretically exact, such calculations must involve an infinite set of modes. In practice, only a finite subset of modes can be used, introducing an error into the seismograms. By systematically increasing the number of modes beyond the highest frequency of interest in the seismograms, we investigate the convergence of full-coupling calculations. As a rule-of-thumb, it is nec- essary to couple modes 1–2 mHz above the highest frequency of interest, although results depend upon the details of the Earth model. This is significantly higher than has previously been assumed. Observations of free oscillations also provide important constraints on the het- erogeneous structure of the Earth. Historically, this inference problem has been addressed by the measurement and interpretation of splitting functions. These can be seen as secondary data extracted from low frequency seismograms. The measurement step necessitates the calculation of synthetic seismograms, but current implementations rely on approximations referred to as self- or group-coupling and do not use fully accurate seismograms. We therefore also investi- gate whether a systematic error might be present in currently published splitting functions. We find no evidence for any systematic bias, but published uncertainties must be doubled to prop- erly account for the errors due to theoretical omissions and regularization in the measurement process. Correspondingly, uncertainties in results derived from splitting functions must also be increased. As is well known, density has only a weak signal in low-frequency seismograms. Our results suggest this signal is of similar scale to the true uncertainties associated with currently published splitting functions. Thus, it seems that great care must be taken in any attempt to robustly infer details of Earth’s density structure using current splitting functions. Key words: Computational seismology; Seismic tomography; Surface waves and free oscil- lations; Theoretical seismology. et al. 2017). However, as this paper will show, computational limita- 1 INTRODUCTION tions have led most such studies to rely upon seismological approx- Our understanding of the Earth’s large-scale interior structure and imations that might be inadequate for imaging density variations. dynamics draws heavily on observations of seismic free oscillations As such, their conclusions must be approached with considerable (‘normal modes’). In particular, measurements of free-oscillation caution. spectra at long periods represent one of the few classes of geo- In order to derive information about Earth’s interior structure physical observable with meaningful sensitivity to lateral variations from seismic data, we must search for models which yield synthetic in density within the Earth. As such, free oscillation spectra have (predicted) observables which agree with those actually measured. played a central role in attempts to determine the nature of the Many strategies exist for performing this search, but all rely on the ‘Large Low Shear Velocity Provinces’ that have been identified presumption that synthetic observables are computed accurately. within the lowermost mantle, and in attempts to identify the relative The theoretical basis for computing synthetic normal mode spectra importance of thermal and chemical effects as driving forces for is well-developed, relying on a technique known as ‘normal mode mantle convection (e.g. Ishii & Tromp 1999; Trampert et al. 2004; coupling theory’ (e.g. Dahlen 1968, 1969; Woodhouse & Dahlen Lay 2007;Mosca et al. 2012; Moulik & Ekstrom 2016; Koelemeijer 1978; Woodhouse 1980; Woodhouse & Girnius 1982;Park 1986, The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 Resolvability of Earth’s 3-D Density 59 1987, 1990; Romanowicz 1987; Tromp & Dahlen 1990; Lognonne´ foregone conclusion—the validity of any approximation is entirely 1991; Hara et al. 1991, 1993; Um & Dahlen 1992; Deuss & context-specific—but the question must clearly be addressed. The Woodhouse 2001; Al-Attar et al. 2012;Yang&Tromp 2015;Lau first to address this question were Resovsky & Ritzwoller (1998). et al. 2015). Unlike some earlier studies, our calculations are based Being limited in computational resources, they estimated the the- on the work of Al-Attar et al. (2012) and take the full effect of 3-D oretical error by a transfer function. Having at our disposal an density variations into account, following the theory developed in arbitrarily precise forward solver, we can, for the first time, quantify Woodhouse (1980). Our calculations also allow us to take atten- this theoretical error properly. uation fully into account, although in the examples shown in this We therefore perform a simple, but definitive test: we com- paper, only the spherically symmetric attenuation profile of PREM pute accurate synthetic seismograms for the earth model S20RTS (Dziewonski & Anderson 1981) is used. Thus, our framework is (Ritsema et al. 1999), making use of full-coupling. We then mea- physically complete, and allows the seismic response of any Earth- sure splitting functions from this data set as if it were observed data. like body to be expressed exactly. However, the resulting expres- We also compute the ‘true’ splitting functions, by calculating the sions involve an infinite series expansion. This must be truncated appropriate radial integrals for the known input model. Comparing in any computational implementation of the theory, introducing a the two, we find noticeable discrepancies. Repeating the experiment truncation error within the synthetic spectra. for synthetic spectra using the self-coupling approximation allows Various common truncation schemes exist, motivated from physi- us to separate the contribution from theoretical approximations and cal principles. Computational costs scale unpleasantly with the num- regularization in the uncertainty of the measured splitting func- ber of terms retained in the series, so early studies had little choice tions. Finally repeating the experiment for models where density but to adopt simplifying strategies. The most limiting approximation perturbations are set to zero, we conclude that splitting functions is known as ‘self-coupling’, while ‘group-coupling’ incorporates are contaminated by a significant approximation error, a similar additional terms (we explain the distinction in more detail below). regularization error and contain a density signal of the level of the However, a number of studies, starting with Deuss & Woodhouse total measurement uncertainty. (2001), have demonstrated that neither can be relied upon within the frequency range of interest. In particular, Al-Attar et al. (2012) highlighted that use of either the self- or group-coupling approxi- 2 THEORETICAL BACKGROUND mations when computing synthetic spectra leads to errors of similar Imaging the Earth’s interior involves searching for the model param- magnitude to the signal likely to be attributable to lateral variations eters which can match observed seismic data. We therefore require in density within the Earth. Subsequently, a comprehensive study by a forward model, or a way to compute the observations that we Yang & Tromp (2015) concluded that self-coupling is ‘marginally would expect to see if the earth had any given structure. Ideally, this acceptable’ when considering spectra below around 1.5 mHz, and forward model needs to encapsulate the full physics of wave propa- ‘unacceptable’ above this point. Group-coupling was found to also gation. If it does not, the seismograms predicted for each candidate yield significant errors, albeit at a lower level than in self-coupling. structure will be inaccurate, and the results of the imaging process In principle, one could envisage improving group-coupling strate- may be biased by this systematic error (for a full discussion, see gies by computing coupling strengths between modes—perhaps Valentine & Trampert 2016). based on the work of Park (1987)—and ensuring that all signifi- cant interactions are accounted for. However, the notion of coupling strength is itself based on single-scattering approximations, and it 2.1 The forward problem is not clear how accurate the results would be. The optimal approach is therefore to adopt a strategy known as ‘Normal mode seismology’ represents one framework for con- ‘full-coupling’. Despite the name, this continues to involve a trun- structing a seismological forward model, particularly suited to low- cation of the infinite series, but whereas self- and group-coupling frequency, global-scale simulations. The theory incorporates all im- are based around the presumption that the series is dominated by portant physical effects, and can be used to produce seismograms only a few terms, full-coupling aims to retain all terms that have that faithfully reproduce physical reality. However, the computa- any potential to contribute significantly; again, a more precise def- tional resources required are significant, rendering full calculations inition will be given in due course. Building on the existing body intractable for early studies. As a result, a variety of simplifications of work, the first aim of the present paper is to identify more pre- were developed. These can greatly reduce the computational costs, cisely the conditions under which full-coupling can be regarded as at the expense of a degradation in accuracy. A brief, self-contained ‘sufficiently accurate’. account of the underlying theory is given in Appendix, emphasizing This information then permits us to perform high-quality syn- how and why these various approximations can be introduced. Here, thetic experiments to address the second question: can robust in- we sketch some key points. formation about earth structure be inferred from measurements of ‘splitting functions’ (Woodhouse & Giardini 1985; Giardini et al. 1987, 1988)? These are a particular representation of information 2.1.1 Reference eigenfunctions derived from seismic spectra, and they are an important ingredi- ent in many current models of Earth structure (e.g. Ritsema et al. For any spherically symmetric reference model (e.g. PREM, 1999; Ishii & Tromp 1999; Trampert et al. 2004; Ritsema et al. Dziewonski & Anderson 1981), it is relatively straightforward 2011;Mosca et al. 2012; Moulik & Ekstrom ¨ 2016; Koelemeijer to compute eigenfunctions S and their associated eigenfrequen- et al. 2017). However, the theory that connects splitting functions to cies ω ; each pair (S ,ω ) represents one ‘normal mode’ (e.g. k k k structural parameters relies upon precisely the same assumptions as Woodhouse 1988). Conventionally, four indices are used to describe are inherent to self- and group-coupling. Given that these have been the modes: k → (q, l, m, n), with q denoting the ‘class’ (principally demonstrated to introduce significant bias into synthetic spectra, ‘spheroidal’, q = S, or ‘toroidal’, q = T ); l, the angular degree, and are splitting function inferences meaningful? The answer is not a m, the azimuthal order, identifying the spherical harmonic Y (θ, φ) lm Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 60 F. Akbarashrafi et al. which describes the lateral form of S ;and n, the overtone number, Motivated by this, the ‘self-coupling approximation’ modifies the identifying its radial behaviour. definition of M, introducing an over-riding assertion that The eigenfrequency can be shown to depend only upon q, l and n. M (ω) = 0if S and S come from different multiplets. (3b) ij i j For each (q, l, n) combination, there are therefore 2l + 1 eigenfunc- tions that have a different spatial form, but share a common eigen- This results in M taking a block-diagonal form, with each multiplet frequency. This degenerate set is often referred to as a ‘multiplet’, q contributing a block of dimension (2l + 1) × (2l + 1). Solution n l and identified by a label in the form q (e.g. S ). The constituent n l 0 2 of eq. (2) is then much easier, as it decouples into separate calcula- eigenfunctions of each multiplet are then referred to as ‘singlets’. tions for each multiplet. Of course, the resulting solution is only an approximation to the one that would be obtained if the full-coupling matrix were used. 2.1.2 Computing synthetic seismograms The ‘group-coupling’ approximation allows for limited interac- tion between multiplets that are adjacent in frequency. The spectrum A central theme of normal mode seismology is that these eigenfunc- is divided into ‘groups’, with each group containing one or more tions, obtained for a spherically symmetric reference model, may be used as a convenient set of basis functions for representing vector multiplets. Studies may differ in the groupings used. Then, eq. (3b) fields within the (laterally heterogeneous) Earth. In particular, the is relaxed slightly, so that instead seismic waveforms s(x, t) can be expressed in the form (cf. eq. A5) M (ω) = 0 unless S and S belong to the same group. (3c) ij i j s(x, t) = u (t)S (x)(1) k k Again, this simplifies the solution of eq. (2), allowing separate computations for each group. Studies have shown that the resulting where S is a vector field representing an eigenfunction of the seismograms are more accurate than those of self-coupling, but reference model. errors remain significant (e.g. Deuss & Woodhouse 2001; Al-Attar There are an infinite number of eigenfunctions S ,and so in et al. 2012; Yang & Tromp 2015). principle this representation of the seismic wavefield requires an infinite summation. However, if seismograms are to be filtered so that they contain no component above a given frequency, ω ,it max 2.1.5 Normal modes of an aspherical earth model is justifiable to neglect modes with eigenfrequency ω  ω .To k max So far, we have outlined how to obtain seismograms in hetero- make this concrete, we define a cut-off frequency, ω ≥ ω ,and c max geneous earth models. It is also possible to compute the normal work with the finite set of modes having ω ≤ ω . Then, the Fourier- k c modes of a general (non-spherically symmetric) earth model (see transformed expansion coefficients u˜ (ω) may be found by solving Appendix A4) and again, group- or self-coupling approximations a system of linear equations (cf. eq. A9) may be adopted if required. We find that the modes are no longer M(ω)u˜ = q˜ (2) degenerate: each singlet within a multiplet has a distinct eigenfre- quency. This property is sometimes referred to as ‘splitting’ of the Here, M(ω) is a frequency-dependent ‘coupling matrix’ which can multiplet, relative to the spherically symmetric reference. Given that be computed for any earth model presenting 3-D variations in den- realistic earth models contain only weak lateral heterogeneity, the sity, elastic constants and attenuation, while the vector q˜ is a repre- typical frequency shifts involved are small, and due to finite-length sentation of the seismic source (typically an earthquake). time series, we do not generally expect to be able to resolve isolated singlets in spectral observations. Nevertheless, aspherical structure 2.1.3 Full-coupling manifests itself in the overall shape of each ‘peak’ in the spectrum of a seismogram. The term ‘full-coupling’ implies that we solve this linear system (i.e. eq. 2) directly. If we do this, we obtain seismograms that are essentially physically complete. They incorporate only one funda- 2.2 From spectra to structure: the splitting functions mental approximation: truncation of the infinite series at frequency Given the ability to compute synthetic spectra, one may contemplate ω . The first goal of this paper is to quantify the effect of this trunca- formulating an inverse problem that aims to infer earth structure tion upon synthetic seismic spectra, and to assess where ω should from observations. This could be implemented in a variety of ways, lie in order for the truncation to have no appreciable impact upon each with pros and cons. One popular approach involves the use of calculations for frequencies below ω . max ‘splitting functions’. These were introduced within Woodhouse & Giardini (1985), and developed in detail by Giardini et al. (1987, 2.1.4 The self- and group-coupling approximations 1988). These original papers are developed under the self-coupling assumption, where each multiplet is regarded as entirely isolated As ω increases, the dimension of the matrix M increases, and thus from every other. This was recognized to be a simplification of the cost of solving eq. (2) grows rapidly. The elements of M are the underlying physics, and therefore an approximation, but it was found by computing integrals of the form (cf. eq. A10) necessary if calculations were to be feasible using the computational ∗ 3 resources then available. M (ω) = S (x) · A(ω)S (x)d x (3a) ij j Making this approximation, the contribution of a given mul- where A(ω) is an integro-differential operator that depends upon tiplet κ = q within seismograms is entirely characterized by the n l frequency, and also upon the earth model. The elements tend to be (2l + 1) × (2l + 1) diagonal sub-block in the coupling matrix arising (κκ) relatively small unless S and S have similar eigenfrequencies. from eqs. (3a & 3b), which we denote M . Similarly, the multi- i j This leads to structure within M, linked to the degeneracy within plet can be directly identified with a specific portion of an observed each multiplet; it is dominated by the elements lying close to the seismic spectrum. Thus, given a suite of observations of the rele- diagonal. vant portion at different locations and from various seismic events, Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 Resolvability of Earth’s 3-D Density 61 one can formulate an inverse problem to estimate the elements of specifies P-wave speed (α), S-wave speed (β), and density (ρ)per- (κκ) M . Splitting functions provide a mechanism for parametrizing turbations at every point with respect to a spherically symmetric this inversion, designed to provide a straightforward link between reference model. Furthermore, we assume that the model is lat- (κκ) M and the underlying earth model. erally parametrized in terms of spherical harmonics to maximum degree L. Thus, the P-wave speed field may be expressed 2.2.1 Splitting functions for isolated groups of multiplets α(r,θ,φ) = α (r)Y (θ, φ)(8) lm lm l=0 The splitting function framework was extended by Resovsky & Ritzwoller (1995) and others, allowing for cross-coupling between with identical expressions for β and ρ. By deriving the exact ex- (κκ ) specific pairs (or small groups) of multiplets with similar frequen- pression for H (e.g. Woodhouse 1980), and then comparing this mm cies as in group-coupling. For simplicity, we describe here the case with eq. (5), it is possible to identify formulae for the kernels K where two multiplets, κ = q and κ = q , are allowed to interact n l n such that with one another, but are assumed to be isolated from the remainder (κκ ) (α;κκ ) (β;κκ ) of the spectrum. Extension to larger groups, or restriction to the c = α (r)K (r) + β (r)K (r) st st st st st self-coupling case, follows straightforwardly. (ρ;κκ ) + ρ (r)K (r)dr. (9) The isolated-group assumption implies that the only relevant part st st of M(ω) is the symmetric diagonal block If the aspherical model is specified in terms of a different set of (κκ) (κκ ) parameters, equivalent expressions can be found. Thus, an earth M M . (4) (κ κ) (κ κ ) model can be obtained by performing a linear inversion upon mea- M M sured splitting coefficients. A particular attraction of this procedure As discussed in Appendix A4, some additional assumptions allow us is the possibility of selecting multiplets that are presumed to have (κκ ) (κκ ) 2 to introduce the splitting matrix, H,suchthat M ≈ H − ω I. sensitivity at particular depths, and hence to produce models fo- We then parametrize this by introducing a set of complex-valued cused on a given region of the Earth (e.g. Koelemeijer et al. 2013). ‘splitting coefficients’ c such that st l+l s (κκ ) (κκ ) (κκ ) 2.2.2 Splitting functions in the real Earth 2 mm t H = ω δ + ω W + γ  c (5) κκ 0  st mm 0 mm ll s t =−s s=l−l All the foregoing discussion is self-consistent, subject to the fun- damental assumption that the (groups of) multiplets are spectrally where the γ are defined in terms of spherical harmonic triple prod- isolated. However, at least in the context of forward-modelling, it ucts, has been conclusively shown that this approach is insufficiently 2π π accurate for modern seismology, and instead high-quality seismic mm t ∗ γ = Y (θ, φ)Y  (θ, φ)Y (θ, φ)sin θ dθdφ. (6) l m st ll s lm calculations must account for interactions throughout the spectrum. 0 0 Various studies have shown that allowing for cross-coupling in split- This integral can be expressed straightforwardly in terms of the ting function analysis improves results (e.g. Resovsky & Ritzwoller Wigner-3j symbols, and vanishes for many combinations of spheri- 1995, 1998; Irving et al. 2008; Deuss et al. 2013). However, it is cal harmonics. Together, the coefficients describe a ‘splitting func- not possible to extend the splitting function approach to ever-wider tion’ coupling groups: the linearization described in Appendix A4 re- l+l s quires the coupling band-width to be small. If, instead, we choose (κκ ) (κκ ) η (θ, φ) = c Y (θ, φ)(7) st st to work with the full frequency-dependent coupling matrix M(ω), s=l−l t =−s it is unclear whether a splitting-function–style parametrization can which may be regarded as a map depicting spatial variations in spec- be usefully introduced. Where does this leave splitting-function tral properties. Given the symmetry of M, three splitting functions studies? are required to fully characterize our isolated pair of multiplets: In effect, splitting functions should be regarded as a derived (κκ) (κ κ ) (κκ ) η , η and η . data type or ‘secondary observables’ (e.g. Cara & Leveque 1987), ´ ˆ Following Woodhouse & Girnius (1982), the contribution of the intended to convey spectral information in a more manageable form. two multiplets within each seismogram can be computed via diag- The measurement procedure should therefore be seen as a data- onalization of eq. (4) (see eq. A13 and the associated discussion). processing operation, or transformation applied to raw seismic data. Thus, it is possible to pose the inverse problem whereby the c are The particular form of this transformation may be motivated by st determined from observed data; typically, the fiducial frequency ω an out-dated assumption, but it remains a valid transformation to and parameters describing attenuation (which we will not discuss adopt. As such, published measurements of splitting functions may further) are also refined within this process. The inversion is found to be taken at face value: they need not be regarded as somehow be weakly non-linear, and is typically tackled using a Bayesian for- ‘approximate’. The splitting function is defined to be the result of a mulation of the least-squares algorithm (as per Tarantola & Valette specified measurement procedure. 1982). Invariably, the data is insufficient to uniquely constrain a so- However, if we are to adopt the isolated-multiplet assumption in lution, and it is therefore necessary to introduce prior information the measurement procedure, eq. (9) provides only an approximate (i.e. regularization). Discussion of the extent to which this influ- relationship between the measured splitting functions and the un- ences results may be found in numerous sources (e.g. Resovsky & derlying earth structure. Alternatively, one may take the view that Ritzwoller 1999; Kuo & Romanowicz 2002; Pachhai et al. 2016). eq. (9) is the defining property of the splitting functions: in this case, In the isolated-group approximation, the splitting coefficients the measurement procedure is deficient, since it does not properly can be straightforwardly related to structural parameters. For the account for physical effects in wave propagation. Regardless of the purposes of illustration, we assume that the aspherical earth model viewpoint adopted, ‘measured’ and ‘predicted’ splitting functions Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 62 F. Akbarashrafi et al. Figure 1. Position of the stations (red stars) and the Bolivian event (beach ball) used in this study. Station GRFO is highlighted and shown by a blue star. differ, due to approximations within the measurement procedure. may vary considerably between different locations and experiments, Thus, the second key goal of this paper is to assess the impact of limiting the generality of this definition. In addition, the earthquake this approximation: can models produced using eq. (9) nevertheless source itself is also not perfectly well known, although it is appar- be interpreted usefully, or are they too severely contaminated by ent from inspection of eq. (2) that this uncertainty can be absorbed systematic errors? into the observational data uncertainty, by adding the corresponding If fully accurate earth models are to be produced from splitting variances. function data, it is necessary to replace eq. (9) and properly account Alternatively, given that we typically wish to fit seismograms to for multiplet interactions throughout the spectrum. Unfortunately, infer Earth structure, another definition of an ‘acceptable’ truncation it is not clear that this can be achieved by simple modifications of error is that the error should be much smaller than the signal we the modal depth functions K(r). It seems likely that the most fruitful wish to infer. This is perhaps more straightforward, since it does course would involve obtaining sensitivity kernels via adjoint tech- not depend upon the vagaries of seismic noise. Of course, the signal niques within a full-coupling framework, allowing the full range of itself must be above the noise level if we are to make meaningful physical effects to be accounted for. However, we do not pursue this inferences. idea further within the current paper. Indeed, given the capacity to We discuss both cases below and show that they can lead to compute exact sensitivity kernels for spectra using adjoint meth- different conclusions depending on the earth model. To quantify ods, one may question whether there is value in continuing to adopt these ideas, we will only consider synthetic data calculated for the two-stage inversion procedure inherent to splitting function various earth models. Of the three major seismic parameters (S- studies. wave speed, P-wave speed & density), density is known to have the least influence on seismograms, and is a parameter of considerable importance for understanding global geodynamics. We therefore 3 TRUNCATION ERRORS IN focus on the density signal when assessing whether truncation errors FULL-COUPLING SEISMOGRAMS are ‘negligible’. First, we consider the problem of computing seismograms using full-coupling. As we have described, a theoretically complete ex- pression for normal mode seismograms requires summation of an 3.2 Observations using S20RTS infinite set of seismograms. To make the computations tractable, We solved eq. (2) exactly, taking all 3-D effects of varying elas- we must truncate this modal expansion (eq. 1), and only include tic constants and density into account, using the method described those modes with eigenfrequencies below some cut-off, ω .This in Al-Attar et al. (2012), and computed the same set of seismo- truncation will inevitably introduce errors into the computed seis- grams several times varying ω . In all cases, we used the Global mograms. How do we choose ω to ensure that synthetic spectra are CMT Catalog source mechanism for the 1994 Bolivian event (event accurate below some frequency ω ? max code: 060994A) and a global distribution of 129 stations based on those within the Global Seismic Network and IRIS/IDA Seis- mic Network (see Fig. 1). As basis functions, we used eigen- 3.1 The scale of acceptable error functions calculated in in the spherically symmetric model PREM In order to address this question, we must somehow define the mag- (Dziewonski & Anderson 1981). The coupling matrix, explicitly de- nitude of error that we consider acceptable. One way to go about fined in eq. (A10a), is constructed for the shear wave speed model this is to consider the noise in the observed seismograms that we S20RTS (Ritsema et al. 1999), which provides values for δ ln V with ultimately wish to fit. To allow meaningful analysis, the trunca- respect to PREM, and where we added scaled compressional wave tion error should be well below that noise. However, noise levels speed perturbations using δ ln V = 0.5 δ ln V and scaled density p s Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 Resolvability of Earth’s 3-D Density 63 −3 x 10 Effect of source errors 7 Effect of omitting density Truncation error at 5mHz Truncation error at 4mHz Truncation error at 3mHz Reference spectrum (coupling up to 6mHz) 0.5 1 1.5 2 2.5 3 Frequency(mHz) Figure 2. Amplitude spectrum of vertical acceleration at station GRFO (black line) calculated by coupling all modes having frequencies less than ω = 6mHz. We refer to this as the ‘reference’ spectrum, s . Also shown are amplitude spectra of the difference between this reference spectrum and those produced using ref lower values of ω (blue lines), which we describe as the ‘truncation error’ for that cut-off. The red line shows the effect of removing any contribution from lateral density variations within the reference spectrum. Note that this is of broadly similar scale to the truncation error at ω = 3 mHz. The green line shows the effect of perturbing the seismic source with three standard deviations on the reported source uncertainties. perturbations using δ ln ρ = 0.3 δ ln V . These scaling relations are Rather than plotting similar figures for every station around the commonly used and are obtained from mineral physics considera- globe, we sought a method to assess the truncation errors more tions (e.g. Karato 1993). On top of the mantle model, we imple- systematically. We defined a measure of the relative spectral misfit mented the simple crustal model from Woodhouse & Dziewonski within any given frequency range (ω , ω ) 1 2 (1984) and included rotation and ellipticity. In the following, we 2 2 |s(ω) − s (ω)| dω ref will consider vertical component spectra. 1 ξ(s,ω ,ω ) =  . (10) 1 2 2 2 We chose ω = 3 mHz and used ω = 3, 4, 5, and 6 mHz, |s (ω)| dω ref max c increasing the number of coupled singlets (toroidals and spheroidals We then plot histograms of ξ evaluated for every station in our combined) from approximately 2000 to 13 000. We treated the data sets (Fig. 3). Again, we compare the reference spectra to those results from ω = 6 mHz as a set of ‘reference’ synthetic spectra, truncated at ω = 3, 4, and 5 mHz, and to those obtained by using s , and looked at amplitude spectra of the differences s − ref (ω =3mHz) ω = 6 mHz, but omitting lateral density. If we wish to ensure s , etc., to provide an indication of the convergence of the full- ref that the truncation error is consistently smaller than the density coupling calculations (Fig. 2). These difference spectra represent signal, we require the corresponding histograms to not overlap. We the truncation error, and are referred to as such in the figure. It is clearly see that a truncation at 3 mHz will generally lead to errors clear that a truncation at ω = ω can still have substantial errors, c max of a similar scale to the signal anticipated from density. Even for but they decrease rapidly as ω is increased. As far as we know, most a truncation at 4 mHz some overlap of the distributions occurs for published studies that adopt a ‘full-coupling’ strategy have used frequencies close to ω . It therefore seems necessary to ensure max ω ≈ ω . c max that the truncation frequency is significantly above the maximum To provide a sense of scale for this truncation error, we repeated frequency of interest (i.e. ω  ω ). These results are broadly in c max the calculation using ω = 6 mHz, but using a version of S20RTS accordance with what we might expect based on the knowledge of where all lateral density heterogeneity has been omitted (i.e. with the coupling-matrix elements. Generally, modes close together in density structure as given by PREM). The amplitude spectrum of frequency interact most strongly, so truncation errors ought to grow s − s is also shown in Fig. 2 and gives some indication of the ρ¯ ref as we approach the cut-off frequency. magnitude of seismic signal that might be attributable to density To quantify the truncation effects further, we defined a cumulative structure. Perhaps surprisingly, we observe that this is similar to misfit for the entire synthetic data set the magnitude of errors arising from truncation at ω = ω .To c max estimate the effect of source uncertainties, we make use of those (i) (i)   2 |s (ω ) − s (ω )| dω i 0 ref reported in the Global CMT Catalogue, although we note that these ψ(ω) = (11) ω (i) |s (ω )| dω i ref are likely to be an under-estimate (Valentine & Trampert 2012). 0 We therefore implement a large source variation by perturbing each where the summation is over all the stations in the data set. We then parameter by three times the corresponding reported standard devi- plot this quantity as a function of frequency for the various cut-off ation, with positive or negative perturbations chosen randomly, and frequencies (Fig. 4), and for the data set where density heterogeneity again computed spectra. Given the other effects seen in Fig. 2,the has been omitted. There are several ways of reading this figure. First, source effect is very modest. We will therefore neglect its contribu- focusing on the curves for the various choices of ω : suppose we tion in the discussion below. want a precision of (say) 0.1 per cent of the total energy in the signal Amplitude Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 64 F. Akbarashrafi et al. (0.2−3)mHz (0.2−2.5)mHz (2.5−3)mHz 60 60 60 Truncation at 3mHz Truncation at 4mHz Truncation at 5mHz 40 40 40 Omitting density 20 20 20 0 0 0 0.00 0.02 0.04 0.06 0.08 0.10 0.00 0.02 0.04 0.06 0.08 0.10 0.00 0.02 0.04 0.06 0.08 0.10 Spectral misfit Figure 3. Histograms showing relative spectral misfit, ξ (eq. 10), across all stations. Again, we compare reference spectra (ω = 6 mHz) to those truncated at lower frequencies (ω = 3, 4, 5 mHz), with all calculations being performed in the model S20RTS (where we incorporated V and density heterogeneity c p as described in the main text). Results from comparing the same reference spectra to those calculated in a version of S20RTS without any lateral density heterogeneity are also shown. Truncation at 3mHz Truncation at 4mHz Truncation at 5mHz Omitting density 0.01 0.001 0.0001 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 Frequency [mHz] Figure 4. Cumulative misfits defined by eq. (11), relative to reference spectra computed in S20RTS, for three different cut-off frequencies and for the omission of lateral density heterogeneity. The horizontal grey line denotes a misfit of 0.1 per cent, as discussed in the main text. (this choice should be guided by knowledge of the uncertainties 3.3 Model-dependence of the truncation error in observed seismic spectra and the seismic source). In this case, We expect results to vary depending on the model, and in partic- truncation at ω = 3 mHz will produces seismograms accurate up ular on the power spectrum of heterogeneity. Coupling effects are to ω = 2 mHz (as shown by a horizontal line on the figure). max somewhat analogous to scattering, and it is well-known that short- For a truncation at ω = 4 mHz, seismograms are accurate up to wavelength heterogeneity promotes multiple-scattering effects. We 3 mHz within the same precision. It appears that we roughly need might reasonably expect to observe similar effects in coupling cal- coupling of 1 mHz higher than the highest frequency of interest culations, with stronger short-wavelength structure necessitating in the seismogram. Alternatively, if we are interested in imaging relatively higher values of ω . S20RTS has modest amplitudes and a Earth’s density structure, we should require truncation errors to be ‘red’ spectrum, dominated by long-wavelength structure. The model at least an order of magnitude smaller than the density signal in the resolution is far below what the parametrization of the model would spectra, and we see roughly that the cut-off frequency should be at allow, ranging from 2000–4000 km horizontally to 250–750 km ver- 4 mHz for any frequency ω ≤ 3 mHz. These values are based on max tically (Ritsema et al. 2004). It is therefore important to investigate an assumption that Earth’s heterogeneity is similar to that present in how this is dictating the findings in the previous subsection. our implementation of S20RTS. The reader may choose her/his own To provide some insight into this question, we generated a ran- desired precision and adapt the conclusions accordingly, by direct dom model using the same parametrization as S20RTS—a spher- reading of Fig. 4. ical harmonic expansion to degree 20 laterally and 21 splines Number of stations Cumulative misfit Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 Resolvability of Earth’s 3-D Density 65 (0.2−3)mHz (0.2−2.5)mHz (2.5−3)mHz 60 60 60 Truncation at 3mHz Truncation at 4mHz Truncation at 5mHz 40 40 40 Omitting density 20 20 20 0 0 0 0.00 0.02 0.04 0.06 0.00 0.02 0.04 0.06 0.00 0.02 0.04 0.06 Spectral misfit Figure 5. Histograms of relative spectral misfit across all stations defined by eq. (10) for a random model. Truncation at 3mHz Truncation at 4mHz Truncation at 5mHz Omitting density 0.01 0.001 0.0001 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 Frequency [mHz] Figure 6. Cumulative misfits defined by eq. (11) for a randomly generated, spectrally white earth model. Compare Fig. 4. vertically. For δ ln V we drew random numbers uniformly in the Clearly the relation between ω , ω and the earth model is s c max range [−0.01, 0.01] giving an rms-amplitude of about 12 per cent not straightforward. Looking at cumulative misfits (Fig. 6), we see at all depths. We continue to assume that V and ρ are scaled ver- that white models appear to yield larger truncation errors than their sions of V , as with S20RTS. The rms-amplitude of this random red counterparts, as expected. Again, considering the data sets that model is between 5 and 38 times stronger than S20RTS depending include density heterogeneity, we can imagine seeking a precision on depth. This produces a model of true degree-20 lateral resolu- of 0.1 per cent of the total energy in the signal. A truncation at tion, a white spectrum, and much higher vertical resolution. We re- ω = 3 mHz now only produces seismograms accurate up to 1 mHz. peated all the full-coupling calculations already described, using this For a truncation at ω = 4 mHz, seismograms are accurate up to model. around 2 mHz within the same precision. It appears that we now need Again, we produce histograms of relative spectral misfit (Fig. 5) coupling to be at least 2 mHz higher than the highest frequency of and cumulative misfit (Fig. 6). From Fig. 5, we see that now the interest in the seismograms. If, however, we define the acceptable truncation effects are much stronger and overall we need to couple truncation error as being an order of magnitude smaller than the more modes than was necessary in S20RTS. In other words, more signal arising from density heterogeneity, we now see that coupling power in the structural model, and especially at short-wavelengths, up to 5 mHz is required for ω ≤ 3 mHz. Again it is straightforward max necessitates higher values of ω . It is also obvious that the contri- to adapt the conclusions for any other desired precision. bution of density to the seismograms is much lower in the random model than in S20RTS, despite the fact that the rms-amplitude of 3.4 Summary the model is much larger. This is because density imprints itself differently into the seismograms, compared to the model with a red The key conclusions we draw from these experiments may be sum- spectrum. marized: Number of stations Cumulative misfit Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 66 F. Akbarashrafi et al. S S S S S 2 6 2 12 2 13 5 3 0 2 CMB CMB CMB CMB CMB ICB ICB ICB ICB ICB S S S S S 0 5 0 6 0 7 1 4 1 7 CMB CMB CMB CMB CMB ICB ICB ICB ICB ICB S S S 1 8 1 9 1 10 Vs CMB CMB CMB Vp density ICB ICB ICB (κκ) Figure 7. Self-coupling sensitivity kernels K (as in eq. 9) corresponding to the angular order s = 0 splitting coefficient for our selected modes. The solid red line denotes V sensitivity, the solid black line V sensitivity and the dotted black line gives density sensitivity. The modes are ordered somewhat arbitrarily, p s with a general trend from mostly upper-mantle–sensitive to mostly lower-mantle–sensitive. (i) Accurate full-coupling spectra require the cut-off frequency 1999;Masters et al. 2000; Deuss et al. 2011, 2013) still represent ω to be significantly above the maximum frequency present in the the main ‘derived’ data set for long-period seismology. Measured seismograms, ω . To our knowledge, this criterion has not been centre frequencies (i.e. ω in eq. 5) are taken as representative of max 0 met in previously published studies based on full-coupling. the Earth’s spherical average structure (e.g. de Wit et al. 2014)and (ii) Accurate seismograms up to 3 mHz in earth-like models (hav- splitting function coefficients feature as a constraint in many current ing a red spectrum) require truncation at no less than 4 mHz. global tomography models (e.g. Ritsema et al. 2011; Moulik & (iii) Generally truncation errors are stronger at higher frequen- Ekstrom ¨ 2016). Splitting functions have also been used for targeted cies, as might be expected. studies investigating specific aspects of Earth structure, since one (iv) Effects seem more pronounced as models gain more power, can examine the kernels from eq. (9) and identify modes which especially at shorter wavelengths. In models with white spectra, have the desired sensitivity (see Fig. 7). This has particularly been accurate seismograms up to 3 mHz require a truncation level of at exploited to investigate deep-earth structure, including inner core least 5 mHz. studies (beginning with Woodhouse et al. 1986), targeted lower mantle studies (e.g. Koelemeijer et al. 2013, 2017) and density Overall, these observations are all broadly consistent with expec- inferences (e.g. Ishii & Tromp 1999; Resovsky & Trampert 2003). tations based on theoretical arguments. However, the effects are now It is common for splitting function catalogues to publish uncer- quantified, and we have obtained rules-of-thumb that can be applied tainty estimates on the recovered splitting coefficients and centre- to ensure acceptable accuracy when performing calculations. frequencies. In principle, these can be obtained directly from the Bayesian inversion framework (e.g. Tarantola & Valette 1982). However, it is not clear how to define a meaningful prior for the 4 THE INTERPRETABILITY splitting function inversion, and in any case such analysis would OF SPLITTING FUNCTION neglect the inherent non-linearity of the measurement process. This MEASUREMENTS point has previously been made by Resovsky & Ritzwoller (1998). After a careful analysis of all potential contributions to the uncer- We now turn to consider splitting functions, and the extent to which tainty, they estimated combined uncertainties using regressions per- they can provide useful constraints upon Earth structure. Catalogues formed on synthetic data with simulated theoretical errors and data of splitting function measurements (e.g. Giardini et al. 1987;He& uncertainties. More recently, it has been common to use some form Tromp 1996; Resovsky & Ritzwoller 1998; Durek & Romanowicz. Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 Resolvability of Earth’s 3-D Density 67 of boot-strapping procedure to investigate the degree of constraint group-coupling. We used the same event and station distributions upon results (e.g. Deuss et al. 2013). This uncertainty estimate there- and method (forward theory, cut-off frequencies and regularization) fore reflects any apparent inconsistency between different parts of as those in Deuss et al. (2013). The measured splitting coefficients (κκ) the data set: it may not reflect any systematic effects inherent to the c can then directly be compared to the expected coefficients for st measurement procedure or its implementation. model S20RTS (as obtained from eq. 9). Various studies have highlighted that considering interactions Two factors contribute to any difference between measured and between adjacent modes can alter splitting-function results sig- expected coefficients. First, differences may be the result of the ap- nificantly (e.g. Resovsky & Ritzwoller 1995; Irving et al. 2008) proximations in the measurement procedure: this is the effect we and group-coupling, where appropriate, is now the norm for var- wish to address. However, they may also arise as a result of the reg- ious frequency bands (e.g. Deuss et al. 2013). Studies have also ularization applied within the non-linear optimization procedure. investigated the extent to which splitting function inversions are To separate the two contributions, we also calculated S20RTS seis- well-constrained, typically focusing on aspects of the linearized mograms strictly using the self-coupling approximation, and again inversion scheme, such as regularization, parametrization and start- measured splitting functions from these. In this case, any difference ing model (e.g. Resovsky & Ritzwoller 1999; Kuo & Romanowicz between measured and observed coefficients can be attributed solely 2002; Pachhai et al. 2016). However, we are not aware of any previ- to regularization, as data and measurement procedure now contain ous studies that have taken a more holistic approach, to investigate the same physics. Assuming that all contributions are Gaussian dis- whether splitting functions provide a meaningful representation of tributed, the difference between the covariances corresponding to Earth structure that can be interpreted by eq. (9). As already dis- the two sets of splitting function then isolates the theoretical error. cussed, a potential issue arises because the measurement procedure We focused on a selection of modes, which were also measured is inherently based upon isolated-group approximations (either self- by Deuss et al. (2013), namely: S , S , S , S , S , S , S , 0 2 0 5 0 6 0 7 1 4 1 7 1 8 or group-coupling), although their interpretation (via eq. 9) is not. S , S , S , S , S , S . Fig. 7 shows some of the kernels K, 1 9 1 10 2 6 2 12 2 13 5 3 Our results in Section 3, and those of earlier work (e.g. Deuss & which we used in eq. (9) to predict the splitting functions we ought Woodhouse 2001; Al-Attar et al. 2012; Yang & Tromp 2015), all to retrieve if the measurements were perfect. Some of the chosen indicate that isolated-group seismograms provide a poor approx- modes have sensitivity mainly in the lower mantle, some mainly imation to the true physics of wave propagation in the earth. It in the upper mantle, and some are sensitive throughout the mantle. therefore seems distinctly possible that structural models estimated Our discussion will concentrate on coefficients of angular order 2 via isolated-group synthetics may inherit these inaccuracies. and 4, which carry most of the splitting function energy for the We address this by performing a straightforward synthetic test. modes we selected to analyse. This will also avoid the problem of For any earth model, we compute ‘predicted’ self-coupling split- neglected higher degree structure as the measurements are done to ting function coefficients using eq. (9). We also compute accurate higher orders (Resovsky & Ritzwoller 1998). All our synthetic ex- synthetic seismograms, using full-coupling theory and the results of periments are noise-free because we want to isolate the theoretical the previous section. From these seismograms, we then obtain ‘mea- error without a particular noise model (which would be difficult (κκ) sured’ splitting functions c , following the same procedure as is to choose objectively) propagating non-linearly into the results via st used for real data. We then assess the agreement between measured the regularization. Assuming that splitting function measurements and predicted splitting function values. Ideally, they should agree to are secondary data to be taken at face value, we will put our syn- within the stated uncertainty for the measurement procedure. If this thetic measurements into context by using actual uncertainties from is the case, we can conclude that the measured splitting functions are Deuss et al. (2013) obtained by boot-strapping real earthquake data interpretable using eq. (9), and the measurement uncertainties can and thought to represent only the data noise present in the Deuss be propagated through this expression to provide uncertainties on catalogue. structural parameters. If measured and predicted splitting functions differ by more than the measurement uncertainty—but show no ev- idence of any systematic bias—it may be reasonable to continue to 4.2 Observations from synthetic splitting functions apply eq. (9) for interpretation using increased uncertainty levels to reflect the theory error (as per Tarantola & Valette (1982)). How- The coefficients of the splitting function measurements for mode ever, if systematic biases are found, no meaningful interpretation S , a lower mantle mode, are depicted in Fig. 8.Whenself- 0 6 can be made. coupling seismograms are used, the measured coefficients are in- distinguishable from those calculated for S20RTS via the ker- nels, suggesting that the regularization used for this mode is not affecting the measurement. We therefore only show the cal- 4.1 Measurement of synthetic splitting functions culated coefficients in the figure, but not those measured from As stated in the detailed analysis of Resovsky & Ritzwoller (1998), self-coupling seismograms. There is however a difference with those splitting function measurements depend on theoretical errors, regu- measured using the full-coupling spectra as data. To calibrate this larization, noise in the data and event-station distributions. The final difference against the density signal, we also computed predicted splitting function depends non-linearly on all of these. We want to splitting coefficients for a version of S20RTS without lateral density quantify as precisely as possible the theoretical errors, and so we heterogeneity. While there are differences between the three sets of need to specify the other contributions as accurately as possible. splitting coefficients, these are small compared to the magnitude of To do this, all details matter. We therefore generated a synthetic the coefficients. It is therefore instructive to plot differences. full-coupling data set and measured splitting functions following In the bottom panels of Fig. 8 we plot the difference between the approach of Deuss et al. (2013). We considered these new spec- the splitting functions predicted for S20RTS, and those measured tra to be the observed data, referred to as above as s . We will from full-coupling synthetics (black line). Each coefficient is plot- ref (κκ) only focus on the self-coupling splitting function η (eq. 7), al- ted with error-bars, representing the uncertainty reported for this though measurements in certain parts of the spectrum employed mode by Deuss et al. (2013). We see that most coefficients touch Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 Im(C ) Re(C ) Im(C ) Re(C ) Im(C ) Re(C ) Im(C ) Re(C ) Im(C ) Re(C ) Im(C ) Re(C ) 68 F. Akbarashrafi et al. s=4 s=2 0.5 0.0 Measured from full−coupling synthetics Predicted (including density) −2 −0.5 Predicted (no density) 0.4 coupling effect 0.4 density effect 0.2 0.2 0.0 0.0 −0.2 −0.2 −0.4 −0.4 Figure 8. Observed splitting coefficients for mode S for s = 2and s = 4 (top panels). The full black circles represent measurements made on synthetic spectra 0 6 obtained using full-coupling. The red open circles are predicted splitting coefficients calculated using eq. (9), and correspond to those you would expect if the measurements were perfect. The grey diamonds represent splitting coefficients calculated for mantle model S20RTS without lateral density variations. Bottom panels: difference between those measured from full-coupling spectra and predicted coefficients (black circles) and difference between predicted coefficients without and with density (grey diamonds). The vertical bars correspond to error bars inferred by Deuss et al. (2013) for this mode. the zero-line when the error-bars are taken into account (certainly substantially greater than the stated noise level, it remains compara- when 2 standard deviations or the 95 per cent confidence level is ble in size to the theoretical error in splitting function observations. considered), implying that predicted and measured splitting func- Again, it seems unlikely that meaningful density inferences can be tions agree to within the measurement uncertainty. Thus, for this made from this mode using standard methods. These two examples mode, any systematic errors that arise from the use of self-coupling suggest that both density signal, and measurement error, become are negligible compared to the effect of noise within the data. more significant for modes with shallow sensitivity. To see if this We also plot the difference between splitting functions predicted holds in general, we must examine more modes. in S20RTS, and our version without density (grey line). Again, we Comparing the full suite of splitting function coefficients mea- plot error-bars representing the uncertainty reported by Deuss et al. sured from full-coupling synthetics to those predicted using eq. (9), (2013), and find that many of these also touch the zero-line and we find a very strong correlation (with an R-value of 0.99), with certainly within 2 standard deviations. This suggests that the signal no obvious evidence of any systematic trends. Splitting functions from density lies at, or below the noise level of measured splitting depend non-linearly on the spectra, and a more sophisticated error functions. As a result, it is not straightforward to extract meaningful analysis is required before the possibility of any systematic effects information about density from splitting functions for this mode. In can be fully eliminated. However, our results suggest that a reason- addition, we observe that the signal attributable to lateral density able starting point may be to assume that theoretical measurement heterogeneity appears to be similar in magnitude to the difference errors can be modelled by a zero-mean Gaussian distribution, al- between observed and predicted splitting coefficients. Thus, even lowing them to be straightforwardly incorporated into measurement if noise levels could be reduced significantly, it would appear that uncertainties. current splitting function measurement procedures are insufficiently To interpret our observations, it is then most instructive to show accurate to make use of the density signal. the reduced χ misfits for all measured coefficients of a given split- Similar results are seen in Fig. 9, which shows results for the ting function with respect to the coefficients predicted for S20RTS upper mantle mode S . In this example, we see differences which using eq. (9) 2 12 are more significant relative to uncertainties. Errors due to the use pred. meas. c − c of self-coupling in the measurement procedure are now markedly st 1 st χ = (12) greater than the quoted uncertainty, which accounts only for effects N σ st s,t due to noise in the data, even within 2 standard deviations for some coefficients. This is likely to affect any interpretation based on eq. 9. where σ are the uncertainties from Deuss et al. (2013) for the st Although lateral density heterogeneity has a signal that is now also given mode and N represents the number of coefficients in the Cst(μHz) Cst(μHz) Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 Im(C ) Re(C ) Im(C ) Re(C ) Im(C ) Re(C ) Im(C ) Re(C ) Im(C ) Re(C ) Im(C ) Re(C ) Resolvability of Earth’s 3-D Density 69 s=4 s=2 Measured from full−coupling synthetics Predicted (including density) Predicted (no density) −2 −4 1.5 coupling effect 1.0 1.0 density effect 0.5 0.5 0.0 0.0 −0.5 −0.5 −1.0 −1.0 Figure 9. As Fig. 8, but for mode S . 2 12 0.1 full−coupling seismograms self−coupling seismograms 0.01 S S S S S S S S S S S S S 2 62 12 2 13 5 30 20 50 60 71 41 71 81 91 10 Figure 10. Reduced χ misfit defined by eq. (12) for measurements on full- (black squares) and self-coupling (open squares) seismograms using all parameters as determined by Deuss et al. (2013). corresponding splitting function. This statistical measure allows re- that can be placed in models derived partly or fully from these sults to be summarized for easy interpretation; χ = 1 indicates splitting function data. that the difference between measured and predicted values is, on As has already been discussed, discrepancies between measured average, of similar size to the uncertainties. First, we consider the and predicted splitting coefficients may arise due to reliance upon splitting function measurements made upon synthetic seismograms approximations in the measurement procedure, or as a consequence computed using full-coupling (see Fig. 10). For most modes, the of regularization. To compare these effects, Fig. 10 also shows 2 2 value of χ is well above one, indicating that observations differ χ for splitting functions measured from synthetic seismograms from predicted values by an amount significantly greater than cur- generated using self-coupling. In this case, seismogram genera- rently published data uncertainties for splitting coefficients. In other tion and splitting function measurement employ a consistent the- words, the Deuss catalogue underestimates the true uncertainty on oretical framework. Thus, any differences here arise solely from splitting coefficients. This in turn may impact upon the confidence regularization. For several modes, we see that χ obtained from Cst( Hz) Reduced χ misfit Cst( Hz) Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 70 F. Akbarashrafi et al. 0.1 δ in V 0.01 δ in V δ in 0.001 S S S S S S S S S S S S S 2 62 12 2 13 5 30 20 50 60 71 41 71 81 91 10 Figure 11. Reduced χ misfit defined by eq. (12) for model S20RTS, except with δ ln V = 0 (black squares), δ ln V = 0 (open squares) and δ ln ρ = 0(grey s p diamonds). This provides an indication of the average strength of the respective fields within the splitting functions for the 13 modes analysed. self-coupling synthetics is similar to that for full-coupling syn- uncertainty in the Deuss catalogue by a factor of 2 to use a more thetics. This suggests that the discrepancy between prediction and representative number for the true measurement uncertainty. observation in these cases is mostly attributable to regularization. It While the quoted numbers reflect the contribution to uncertainty is possible that reducing the damping applied during the measure- in the Deuss catalogue, can they be used for other catalogues? ment procedure would improve the agreement between predictions The answer is no, because regularization errors and the data uncer- and observations. However, this might also be expected to increase tainty translate differently into coefficient uncertainty, depending the sensitivity to data noise, resulting in increased error-bars on the intimately on all the details on the measurement procedure (event- recovered coefficients (Resovsky & Ritzwoller 1998). It is therefore station distribution, regularization, number of iterations, cut-off fre- not obvious that a net benefit would be gained; re-analysis of real quencies, . . . ). Because we do not have this information available data to investigate this rigorously is beyond the scope of the present for other catalogues, we cannot give detailed results for those. We study. can however say Resovsky & Ritzwoller (1998) made a real effort to Assuming that all contributions to the uncertainties are Gaussian quantify the true uncertainties given the computational restrictions distributed, covariances will simply add up, in other words the to- of 30 years ago. tal posterior covariance in actually measured coefficients may be To put these results into perspective, it is instructive to look at po- expressed as C = C + C + C ,where C is the theoretical tential information on the Earth’s structure contained in the splitting tot th reg d th uncertainty we wish to quantify, C is due to regularization in the functions. We have already described how we generate a version of reg regression and C = σ I corresponds to data uncertainty including S20RTS omitting any lateral density heterogeneity, and use this to the influence of the event-station distribution. This decomposition estimate the magnitude of signal attributable to density. Similarly, is similar to that used in Resovsky & Ritzwoller (1998). The uncer- we can generate seismograms in versions of S20RTS that assume tainties quoted in Deuss et al. (2013) were determined using boot- there is no lateral V heterogeneity (having only lateral variations in strapping: randomly eliminating stations, and also whole events. V and ρ), or no lateral V structure (only V and ρ). In each case, we p p s This was done for a fixed regularization. Thus, these uncertainties can compute expected splitting coefficients using eq. (9), and obtain therefore represent the data covariance, C , including the particular χ values as in eq. (12). When we plot these in Fig. 11, we see that event-station distribution used to create the Deuss catalogue. the signal of V in the splitting functions is 2-3 orders of magnitude Our measurements based on self-coupling seismograms (see larger than the observational uncertainties for most modes, while Fig. 10) do not reflect the calculated splitting functions perfectly. that of V is on average one order of magnitude larger and ρ mostly Because these seismograms are noise free, the discrepancies may be of the same order as the uncertainties. This is interesting in itself, a measure of the regularization noise, which can be expressed as a as it shows that most modes do not contain any significant informa- multiple of the data noise, C = ασ I.Valuesof α can be directly tion on density, regardless of whether we can accurately measure reg read from Fig. 10. In the same figure we also display the result for their splitting function or not. The reason for this might be that measurements on full-coupling seismograms. These represent the wave speeds mostly affect the phase of a seismic spectrum, whereas discrepancies due to regularization and the theoretical error com- density mostly affects its amplitude. When we further consider our bined, which correspond then to C + C = βσ I. Hence β − α is findings that the measured data variances should be multiplied by a th reg a measure of C expressed in units of the covariances in the Deuss factor of 4.3 to represent the total measurement errors, we find that th catalogue. We see that C and C changes for every mode, but as on average the reduced χ will be around 100 for δ ln V , around 4 th reg s a quantitative indication: for the sample of modes we considered, for δ ln V and around 1 for δ ln ρ. on average β − α = 1.3 and α = 2.0. Thus for the Deuss cata- Most studies infer mantle structure using these splitting function logue, the theoretical measurement error, is on average 1.3 times the in a least-squares inversion (e.g. Ritsema et al. 1999; Ishii & Tromp quoted variance and the regularization error is on average 2 times the 1999; Ritsema et al. 2011; Moulik & Ekstrom ¨ 2016; Koelemeijer 2 2 quoted variance, resulting in C = (1 + 1.3 + 2.0)σ I = 4.3σ I. et al. 2017). Given that the V signal is well above the uncertainty tot s d d We therefore recommend users should multiply the published of the splitting functions, shear wave mantle models using such Reduced χ misfit Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 Resolvability of Earth’s 3-D Density 71 data should be fairly robust. However, because the density signal significantly beyond the highest frequency range of interest in the is of the level of the measurement uncertainty, its robustness was seismograms. For seismograms to be accurate up to 3 mHz in typical justifiably questioned in the past (e.g. Resovsky & Ritzwoller 1999; models, coupling to 5 mHz is preferable; 4 mHz may be sufficient Kuo & Romanowicz 2002). Only by moving away from standard for some applications if the Earth has a structure similar to S20RTS. techniques and using full sampling, can density be inferred with Our experimental evidence appears to corroborate theoretical ex- a meaningful uncertainty (Resovsky & Trampert 2003). The latter pectations that truncation errors become larger as short-wavelength study showed that the density model had on average an uncer- structure becomes more dominant. tainty of about 50 per cent of its magnitude, or 100 per cent within Self-coupling splitting functions, often used in the construction a 95 per cent confidence level. This generated lively debates, and of 3-D Earth models, are themselves inherently based on the self- for instance, the question of buoyant versus heavy large low-shear or group-coupling approximations, and may therefore contain sys- wave provinces is still waiting to be settled, although the heavy side tematic measurement errors. By measuring self-coupling splitting (Ishii & Tromp 1999; Resovsky & Trampert 2003; Trampert et al. functions on exact full-coupling seismograms, we showed that the 2004;Mosca et al. 2012; Moulik & Ekstrom ¨ 2016) got a recent recovered splitting function does not match theoretical predictions. boost with an independent study inverting tidal data (Lau et al. There is also a significant regularization contribution to the recov- 2017). Given that the Deuss catalogue uses far more seismograms ered splitting functions. There is no visible apparent bias from any to infer splitting functions than any previous catalogue, and we now of these contributions, but recently published (Deuss et al. 2013) have a good representations of the true measurement uncertainty, standard deviations of data uncertainties for splitting coefficients constructing a mantle model using a sampling algorithm is well need to be increased by a factor of around 2 to represent the total worth repeating. Alternatively, direct spectra inversions using full error within the measurement process. In that case, our results show coupling theory, where density has a higher signal-to-noise ratio, that the information on the 3-D density structure contained in the are becoming computationally feasible. measured splitting functions has the same magnitude as the total un- certainty. The information on the variation of compressional-wave speed is about 4 times that of the uncertainty and the information 4.3 Summary on the shear-wave speed has a signal-to-noise ratio of about 100. To go beyond a statistical treatment of the theoretical error, one Again, our results can be summarized as: route forward may be to simply avoid the use of splitting functions: (i) We confirm that there is a significant theoretical measurement with modern computational resources, it may be feasible to infer error in currently published splitting functions. structural information directly from spectra via full-coupling cal- (ii) There is no obvious bias, but this theoretical error can, to a culations. Alternatively, it may also be possible to continue to use first approximation, be assumed to be random noise. splitting functions, measured using self- or group-coupling approx- (iii) There is also a significant contribution from regularization imations, and instead replace the kernels within eq. (9) by those to the uncertainty in splitting functions. that account for these measuring approximations correctly (e.g. via (iv) Taking the total uncertainty into account, currently published adjoint theory). However, until such studies are undertaken, we splitting functions contain information on V with a signal-to-noise must assume markedly larger uncertainties than reported for mea- ratio of about 100, information on V with a signal-to-noise ratio sured splitting coefficients. Given this, our modelling suggests that of about 4, and information on ρ with a signal-to-noise ratio of Earth’s density signal lies at or below the presumed noise level. In about 1. consequence, there is no prospect of recovering robust estimates (v) Since the density signal is small in currently published split- of density structure. It therefore follows that currently published ting functions, advanced techniques should be employed to infer density models ought to be interpreted with great caution. density structure. (vi) The signal from V is very strong, and its inference is un- ACKNOWLEDGEMENTS likely to be significantly affected by the increased measurement uncertainty in currently published splitting functions. Thus, we do We are grateful for constructive reviews by Philippe Lognonnea ´ nd not believe our results should affect published models of V structure Joseph Resovsky. The research leading to these results has received that may be derived these splitting function measurements. funding from the European Research Council (ERC) under the Eu- ropean Union’s Seventh Framework Programme (FP/2007-2013) These results are specific to the Deuss catalogue (Deuss et al. 2013). grant agreement number 320639 (iGEO) and under the European Using other catalogues requires repeating the work in this study Union’s Horizon 2020 research and innovation programme grant if authors wish to infer density, the exception possibly being the agreement number 681535 (ATUNE). The code used to perform Colorado catalogue (Resovsky & Ritzwoller 1998). We also remark the full-coupling calculations makes extensive use of routines de- that the density signal may be more effectively isolated by working veloped by John Woodhouse, and we thank him for allowing us to directly upon seismic spectra, where the observational uncertainty use them. is expected to be smaller: the results of Section 3 suggest that this should be feasible. REFERENCES 5 CONCLUSIONS Al-Attar, D. & Crawford, O., 2016. Particle relabelling transformations in elastodynamics, Geophys. J. Int., 205, 575–593. Previous studies have shown that self- and group-coupling approxi- Al-Attar, D., Woodhouse, J. & Deuss, A., 2012. Calculation of normal mode mations ought to be regarded as obsolete in the context of computing spectra in laterally heterogeneous earth models using an iterative direct seismograms in 3-D Earth models. By systematically increasing the solution method, Geophys. J. Int., 189, 1038–1046. number of modes in full-coupling calculations, we have demon- Cara, M. & Lev ´ eque, ˆ J., 1987. Waveform inversion using secondary observ- strated that accurate full-coupling requires inclusion of modes ables, Geophys. Res. Lett., 14, 1046–1049. Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 72 F. Akbarashrafi et al. Dahlen, F., 1968. The normal modes of a rotating, elliptical Earth, Geophys. Lognonne, ´ P., 1991. Normal modes and seismograms in an anelastic rotating J. R. astr. Soc., 16, 329–367. Earth, J. geophys. Res., 96, 20 309–20 319. Dahlen, F., 1969. The normal modes of a rotating, elliptical Earth—II. Near- Masters, G., Laske, G. & Gilbert, F., 2000. Matrix autoregressive analysis resonance multiplet coupling, Geophys.J.R.astr.Soc., 18, 397–436. of free-oscillation coupling and splitting, Geophys. J. Int., 143, 478–489. Dahlen, T. & Tromp, J., 1998. Theoretical Global Seismology, Princeton Masters, G., Barmine, M. & Kientz, S., 2011. Mineos User Manual, Com- University Press. putational Infrastructure for Geodynamics, v. 1.0.2. de Wit, R., Kaufl, ¨ P., Valentine, A. & Trampert, J., 2014. Bayesian inversion Mosca, I., Cobden, L., Deuss, A., Ritsema, J. & Trampert, J., 2012. Seis- of free oscillations for Earth’s radial (an)elastic structure, Phys. Earth mic and mineralogical structures of the lower mantle from probabilistic planet. Inter., 237, 1–17. tomography, J. geophys. Res., 117, B06304, doi:10.1029/2011JB008851. Deuss, A. & Woodhouse, J., 2001. Theoretical free-oscillation spectra: the Moulik, P. & Ekstrom, ¨ G., 2016. The relationships between large-scale importance of wide band coupling, Geophys. J. Int., 146, 833–842. variations in shear velocity, density, and compressional velocity in the Deuss, A. & Woodhouse, J., 2004. Iteration method to determine the eigen- Earth’s mantle, J. geophys. Res., 121, 2737–2771. values and eigenvectors of a target multiplet including full mode coupling, Pachhai, S., Tkalci ˇ c, ´ H. & Masters, G., 2016. Estimation of splitting func- Geophys. J. Int., 159, 326–332. tions from Earth’s normal mode spectra using the neighbourhood algo- Deuss, A., Ritsema, J. & van Heijst, H., 2011. Splitting function measure- rithm, Geophys. J. Int., 204, 111–126. ments for Earth’s longest period normal modes using recent large earth- Park, J., 1986. Synthetic seismograms from coupled free oscillations: Effects quakes, Geophys. Res. Lett., 38, L04303, doi:10.1029/2010GL046115. of lateral structure and rotation, J. geophys. Res., 91, 6441–6464. Deuss, A., Ritsema, J. & van Heijst, H., 2013. A new catalogue of normal- Park, J., 1987. Asymptotic coupled-mode expressions for multiplet ampli- mode splitting function measurements up to 10 mHz, Geophys. J. Int., tude anomalies and frequency shifts on an aspherical earth, Geophys. J. 193, 920–937. R. astr. Soc., 90, 129–169. Durek, J.J. & Romanowicz, B., 1999. Inner core anisotropy inferred by direct Park, J., 1990. The subspace projection method for constructing coupled- inversion of normal mode spectra, Geophys. J. Int., 139, 599–622. mode synthetic seismograms, Geophys. J. Int., 101, 111–123. Dziewonski, A. & Anderson, D., 1981. Preliminary reference Earth model, Pekeris, C. & Jarosch, H., 1958. The free oscillations of the earth, in Phys. Earth planet. Inter., 25, 297–356. Contributions in Geophysics in Honor of Beno Butenberg, eds Benioff, Gharti, H. & Tromp, J., 2017. A spectral-infinite-element solution of Pois- H., Ewing, M., Howell, B. & Press, F., Pergamon Press. son’s equation: an application to self gravity,arXiv:1706.00855v1. Resovsky, J. & Ritzwoller, M., 1995. Constraining odd-degree Earth struc- Giardini, D., Li, X.-D. & Woodhouse, J., 1987. Three-dimensional structure ture with coupled free oscillations, Geophys. Res. Lett., 22, 2301–2304. of the Earth from splitting in free-oscillation spectra, Nature, 325, 405– Resovsky, J. & Ritzwoller, M., 1998. New and refined constraints on three- 411. dimensional Earth structure from normal modes below 3 mHz, J. geophys. Giardini, D., Li, X.-D. & Woodhouse, J., 1988. Splitting functions of long- Res., 103, 783–810. period normal modes of the Earth, J. geophys. Res., 93, 13 716–13 742. Resovsky, J. & Ritzwoller, M., 1999. Regularization uncertainty in density Gilbert, F., 1970. Excitation of the normal modes of the Earth by earthquake models estimated from normal mode data, Geophys. Res. Lett., 26, 2319– sources, Geophys. J. R. astr. Soc., 22, 223–226. 2322. Hara, T., Tsuboi, S. & Geller, R., 1991. Inversion for laterally heterogeneous Resovsky, J. & Trampert, J., 2003. Using probabilistic seismic tomogra- earth structure using a laterally heterogeneous starting model: preliminary phy to test mantle velocity-density relationships, Earth planet. Sci. Lett., results, Geophys. J. Int., 104, 523–540. 215,121–234. Hara, T., Tsuboi, S. & Geller, R., 1993. Inversion for laterally heteroge- Ritsema, J., van Heijst, H. & Woodhouse, J., 1999. Complex shear wave neous upper mantle S-wave velocity structure using iterative waveform velocity structure imaged beneath Africa and Iceland, Science, 286(5446), inversion, Geophys. J. Int., 115, 667–698. 1925–1931. He, X. & Tromp, J., 1996. Normal-mode constraints on the structure of the Ritsema, J., van Heijst, H. & Woodhouse, J., 2004. Global transition zone Earth, J. geophys. Res., 101, 20 053–20 082. tomography, J. geophys. Res., 109, doi:10.1029/2003JB002610. Irving, J., Deuss, A. & Andrews, J., 2008. Wide-band coupling of Earth’s Ritsema, J., Deuss, A., van Heijst, H. & Woodhouse, J., 2011. S40RTS: a normal modes due to anisotropic inner core structure, Geophys. J. Int., degree-40 shear-velocity model for the mantle from new Rayleigh wave 174, 919–929. dispersion, teleseismic traveltime and normal-mode splitting function Ishii, M. & Tromp, J., 1999. Normal-mode and free-air gravity constraints measurements, Geophys. J. Int., 184, 1223–1236. on lateral variations in velocity and density of Earth’s mantle, Science, Rogister, Y. & Valette, B., 2009. Influence of liquid core dynamics on 285, 1231–1236. rotational modes, Geophys. J. Int., 176, 368–388. Karato, S.-i., 1993. Importance of anelasticity in the interpretation of seismic Romanowicz, B., 1987. Multiplet-multiplet coupling due to lateral hetero- tomography, Geophys. Res. Lett., 15, 1623–1626. geneity: asymptotic effects on the amplitude and frequency of the Earth’s Koelemeijer, P., Deuss, A. & Ritsema, J., 2013. Observations of core-mantle normal modes, Geophys. J. R. astr. Soc., 90, 75–100. boundary Stoneley modes, Geophys. Res. Lett., 40, 2557–2561. Tarantola, A. & Valette, B., 1982. Generalized nonlinear inverse problems Koelemeijer, P., Deuss, A. & Ritsema, J., 2017. Density structure of Earth’s solved using the least squares criterion, Rev. Geophys. Space Phys., 20, lowermost mantle from Stoneley mode splitting observations, Nat. Com- 219–232. mun., 8, 15241. Trampert, J., Deschamps, F., Resovsky, J. & Yuen, D., 2004. Probabilistic Komatitsch, D. & Tromp, J., 2002. Spectral-element simulations of global tomography maps chemical heterogeneities throughout the lower mantle, seismic wave propagation—I. Validation, Geophys. J. Int., 149, 390–412. Science, 306, 853–856. Kuo, C. & Romanowicz, B., 2002. On the resolution of density anomalies in Tromp, J. & Dahlen, F., 1990. Free oscillations of a spherical anelastic earth, the Earth’s mantle using spectral fitting of normal-mode data, Geophys. Geophys. J. Int., 103, 707–723. J. Int., 150, 162–179. Um, J. & Dahlen, F., 1992. Normal mode multiplet coupling on an aspherical, Lau, H., Yang, H.-Y., Tromp, J., Mitrovica, J., Latychev, K. & Al-Attar, anelastic earth, Geophys. J. Int., 111, 11–34. D., 2015. A normal mode treatment of semi-diurnal body tides on an Valentine, A. & Trampert, J., 2012. Assessing the uncertainties on seismic aspherical, rotating and anelastic Earth, Geophys. J. Int., 202, 1392–1406. source parameters: towards realistic error estimates for centroid-moment- Lau, H., Mitrovica, J., Davies, J., Tromp, J., Yang, H.-Y. & Al-Attar, D., 2017. tensor determinations, Phys. Earth. planet. Inter., 210–211, 36–49. Tidal tomography constrains Earth’s deep-mantle buoyancy, Nature, 551, Valentine, A. & Trampert, J., 2016. The impact of approximations and 321–326. arbitrary choices on geophysical images, Geophys. J. Int., 204, 59–73. Lay, T., 2007. Deep earth structure—lower mantle and D ,in Treatise on Woodhouse, J., 1980. The coupling and attenuation of nearly resonant mul- Geophysics, Vol. 1, pp. 620–654, eds Romanowicz, B. & Dziewonski, ´ A., tiplets in the Earth’s free oscillation spectrum, Geophys.J.R.astr.Soc., Elsevier. 61, 261–283. Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 Resolvability of Earth’s 3-D Density 73 Woodhouse, J., 1988. The calculation of the eigenfrequencies and eigenfunc- algorithms. In a seismological context, one strategy is to adopt a tions of the free oscillations of the earth and the sun, in Seismological basis constructed on a finite-element–style mesh. This results in Algorithms, pp. 321–370, ed. Doornbos, D., Academic Press. computational tools such as SPECFEM3D_GLOBE (e.g Komatitsch & Woodhouse, J. & Dahlen, F., 1978. The effect of a general aspherical pertur- Tromp 2002), where simulation costs scale with the length of time- bation on the free oscillations of the Earth, Geophys. J. R. astr. Soc., 53, series that must be calculated. Normal mode seismology takes a 335–354. different approach: it relies upon basis functions that have global Woodhouse, J. & Deuss, A., 2007. Earth’s free oscillations, in Seismology support, and the resulting algorithms have computational costs that and structure of the Earth, vol. 1 of Treatise on Geophysics, chap. 2, depend largely on the desired frequency range of solution. Both pp. 31–65, eds Romanowicz, B. & Dziewonski, ´ A., Elsevier. approaches are formally complete and can provide faithful simula- Woodhouse, J. & Dziewonski, A., 1984. Mapping the upper mantle: Three- tion of physical reality. However, we note that SPECFEM3D_GLOBE dimensional modeling of Earth structure by inversion of seismic wave- forms, J. geophys. Res., 89, 5953–5986. currently omits self-gravitation effects, as these are challenging to Woodhouse, J. & Giardini, D., 1985. Inversion for the splitting function of incorporate in a spectral element framework, and is therefore not isolated low order normal mode multiplets, EOS, Trans. Am. geophys. suitable for modelling seismograms at the very low frequencies we Un., 66, 300. consider in this paper (although efforts are being made to overcome Woodhouse, J. & Girnius, T., 1982. Surface waves and free oscillations in a this difficulty: see Gharti & Tromp 2017). regionalized earth model, Geophys. J. R. astr. Soc., 68, 653–673. Woodhouse, J., Giardini, D. & Li, X.-D., 1986. Evidence for inner core anisotropy from free oscillations, Geophys. Res. Lett., 13, 1549–1552. Yang, H.-Y. & Tromp, J., 2015. Synthetic free oscillation spectra: an ap- A2 Spherical-earth eigenfunctions praisal of various mode-coupling methods, Geophys. J. Int., 203, 1179– The particular family of basis functions we use are the eigenfunc- tions, or ‘normal modes’, of a spherically symmetric earth model. In this case, the material properties vary only with depth, and eq. (A1) simplifies considerably. The Coriolis term can be discarded, and the APPENDIX: SYNTHETIC SEISMIC operator V takes a simpler form, which we denote by V. Similarly, WAVEFORMS AND SPECTRA—THE we use ρ¯ to highlight that the density distribution is also spherically NORMAL MODE FRAMEWORK symmetric. Neglecting the force term, and working in the frequency A complete description of the theoretical basis for normal mode domain, it can then be shown that the differential equations admit a seismology lies well beyond the scope of this paper. Here, we aim non-trivial solution, S (x), at a discrete set of frequencies ω . Thus, k k to provide sufficient background to allow the reader to appreciate we have the motivations for—and distinctions between—self-, group- and V − ρω ¯ S =0(A2) full-coupling. In particular, we aim to highlight precisely where and k why approximations need to be introduced within the various theo- where S represents the kth normal mode, and ω is its correspond- k k ries, providing insight into potential for systematic errors. Readers ing eigenfrequency. wishing for a more in-depth treatment are directed to a summary pa- This eigenvalue problem can be solved relatively straightfor- per by Woodhouse & Deuss (2007), or the comprehensive account wardly (Pekeris & Jarosch 1958). There are an infinite number contained within Dahlen & Tromp (1998). Of course, our discussion of (ω ,S ) pairs that satisfy eq. (A2), but only a finite set exist k k here is very general: specific implementations of these techniques within any given frequency interval ω ≤ ω ≤ ω .Thus,itisfea- 0 k 1 invariably adopt additional simplifications and restrictions, which sible to compute a complete set of eigenfunctions within the given may impact further upon accuracy. frequency range, and an algorithm described by Woodhouse (1988) provides a means for doing so. This method is implemented, for ex- ample, within the MINEOS software package described by Masters A1 The seismic wave equation et al. (2011). The seismic wavefield within the Earth, s(x, t), can be described by It turns out (see e.g. Pekeris & Jarosch 1958) that distinct families the seismic wave equation. This has the general form of normal modes may be identified. Principally, two classes exist, known as ‘spheroidal’ and ‘toroidal’ modes; in addition, ‘radial’ ∂ ∂ V + W + ρ s = f (A1) modes are sometimes discussed as a separate category, although they ∂t ∂t are essentially a particular subset of the spheroidals. Furthermore, where V represents a linear integro-differential operator that de- distinct sets of toroidal modes may be identified for each isolated pends upon the material properties of the Earth; W is an operator solid region within the structural model, thus two groups of toroidal representing effects due to the Coriolis force; ρ = ρ(x) is the equi- modes may be found for earth models such as PREM, where a solid librium density distribution within the Earth; and f = f(x, t)repre- inner core and mantle/crust region are separated by a fluid outer core. sents the force due to the seismic source. If viscoelastic effects are to The modes take the spatial form of spherical harmonics laterally; be properly accounted for, the operator V is itself a function of time, their radial functions must be found by numerical integration. Thus, and in this case Vs implies a convolution operation. Given specific the index k on S may be recast as a set of four indices, k → (q, realizations of the various quantities, we aim to solve eq. (A1) to l, m, n). Here, q defines the ‘family’ of mode being considered; l obtain synthetic seismic observables. and m are the angular degree and azimuthal order specifying the A general strategy for solving differential equations entails ex- appropriate spherical harmonic Y (θ, φ); and n is the ‘overtone lm panding the solution in terms of a known basis set, thereby reducing number’, which indexes the set of permissible radial functions. the equations to algebraic form. Although the solution remains the For present purposes, the main point of interest is that ω de- same irrespective of the basis chosen, the computational procedure pends only on q, l and n: in a spherically symmetric model, modes required to obtain that solution can vary considerably. Thus, an that differ only by their m-value share the same eigenfrequency. In appropriate choice of basis can be critical to developing efficient the literature of normal mode seismology, it is common to see an Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 74 F. Akbarashrafi et al. individual (q, l, m, n) state described as a ‘singlet’; the (2l + 1) symmetric reference structure. In addition, we implicitly assume degenerate singlets corresponding to a given (q, l, n) triplet are then that all interfaces within the aspherical model remain concentric referred to as a ‘multiplet’. This triplet is often written in the form spheres: in other words, no topography is present, either on the q ,with q being either S (in the case of spheroidal modes), or T free surface or on internal discontinuities (such as the core–mantle n l (toroidals). For brevity, we will continue to make use of the index k boundary). To incorporate such topography, studies have tradition- within this discussion, unless there is reason to do otherwise. ally relied on first-order boundary perturbation theory (Woodhouse & Dahlen 1978; Woodhouse 1980), which is inherently an approx- imation. Recently, an exact treatment of boundary perturbations A2.1 Completeness and orthogonality has been developed by Al-Attar & Crawford (2016), although this theory has yet to be fully implemented. In any case, the effects of The eigenfunctions of a spherically symmetric model can be shown realistic-scale topography are not likely to impact the conclusions to have some important properties, which stem from the self- of this paper, and need not be considered further here. adjointness of the operator V. In particular, the infinite set of S Given that the S form a complete basis, we can express the form a complete basis for vector fields within the earth model, and wavefield s(x, t)inthe form they may be shown to be orthogonal, in the sense that s(x, t) = u (t)S (x), (A5) k k ∗ 3 ρ¯ (x)S (x) · S (x)d x = δ (A3) k jk where the real part is again understood. u represent time-varying where δ represents the Kronecker delta. Here, the integration vol- jk expansion coefficients. Strictly, this is only valid if the summa- ume V is the complete interior of the earth model, and an asterisk tion includes the entire and infinite set of eigenfunctions for the is used to denote complex conjugation. spherically symmetric model. However, we note the similarity with Strictly, these properties are only proven for entirely solid models, eq. (A4), and the fact that realistic aspherical models contain only and we ignore any complications that might arise from the existence relatively weak perturbations to the spherically symmetric reference of ‘undertone modes’ in the Earth’s fluid outer core. These are not model. It therefore seems reasonable to assert that modes with eigen- thought to be problematic in a seismological context, as they exist frequency far in excess of ω are unlikely to contribute within the max at lower frequencies than we typically consider (e.g. Rogister & frequency band of interest. We introduce a cut-off frequency, ω ≥ Valette 2009). ω , and work with the finite set of S having ω ≤ ω . max k k c A2.2 Seismograms in a spherically symmetric model A3.1 Full-coupling Knowledge of the spherical-earth eigenfunctions allows straight- In order to make use of this expansion, eq. (A5), we need to ob- forward synthesis of the wavefield expected within a spherically tain expressions for the coefficients u . Substituting into the wave symmetric body. As shown by Gilbert (1970), and recounted in equation, eq. (A1), we find Woodhouse & Deuss (2007), the three-component wavefield can be reduced to the form (u VS + u˙ WS + ρu¨ S ) = f (A6) k k k k k k s¯(x, t) = E (1 − exp ı ω t)S (x)(A4) k k k ω and therefore, for any eigenfunction S , we obtain k j where the real part is understood and we only consider positive ∗ 3 ∗ 3 u S · VS d x + u˙ S · WS d x k k k k j j times. E represent the modal excitations, and may be obtained by V V calculating certain spatial integrals of the force term f(x, t). The ∗ 3 ∗ 3 process of computing seismograms by this formula is often referred + u¨ ρ S · S d x = S · f d x. (A7) k k j j V V to—for obvious reasons—as ‘mode summation’. In principle, the summation over k includes an infinite number Since the sum only contains a finite number of terms corresponding of modes. However, we see that the dominant contribution of each to frequencies lower than ω , this can equivalently be expressed as mode to the seismogram lies at its eigenfrequency. Thus, if the a matrix equation, seismogram is only required within a finite frequency band, only Vu + W˙u + P¨u = q (A8) modes lying within that band need to be summed, and this makes the summation finite. In practice, we always work with band-limited where the elements of P, W, V and q are obtained by evaluating the seismic data, which is a necessary consequence of digital signal appropriate integrals from eq. (A7). Taking the Fourier transform recording. If the maximum frequency present in our data set is ω , max allows us to eliminate the derivatives, and introduce the matrix it is straightforward—and sufficient—to compute a mode catalogue M(ω), such that containing all modes with ω ≤ ω . k max V(ω) + ı ωW − ω P u˜ = M(ω)u˜ = q˜ (A9) where the frequency-dependence of V arises from the inclusion of A3 Seismograms in an aspherical earth model viscoelasticity within V; strictly, eq. (A8) involves a convolution for We now return to consider general, aspherical earth models, and the first term. The elements of M and q˜ are readily evaluated, and the solution of eq. (A1). Invariably, the aspherical model will be thus obtaining u˜ is a straightforward exercise in linear algebra. specified as a 3-D structure superimposed upon a spherically sym- The ‘full-coupling’ approach involves solving eq. (A9) exactly. metric reference model (often based on the model PREM, as de- In reaching this point, we have introduced only one approximation: scribed in Dziewonski & Anderson 1981). In what follows, we as- we have truncated the modal expansion at frequency ω . To date, sume that the S have been obtained using the relevant spherically there has been little effort to investigate the errors introduced by k Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 Resolvability of Earth’s 3-D Density 75 this truncation, or assess how ω should be chosen. Most of the for a single group, G. Thus, some limited interactions are allowed existing studies that employ full-coupling assume it is sufficient to between multiplets. However, (assuming a sensible ordering of the use ω = ω , but there is little experimental or theoretical justi- modes), both self- and group-coupling cause the matrix M to take c max fication to support this. Thus, one goal of this paper is to explore on a block-diagonal form. Exploiting this structure allows a great how changes to ω impact the seismograms produced, and to iden- reduction in the computational effort required to solve eq. (A9), at tify where the cut-off must be placed to ensure computations are the expense, of course, of a deterioration in accuracy. essentially exact at frequencies below ω . max We note that a number of ‘full-coupling’ studies introduce addi- A4 Normal modes in aspherical models tional approximations within eq. (A9)—similar to those discussed in Appendix A4, below—designed to reduce computational com- Thus far, we have discussed the use of full-, group- and self-coupling plexity. It has been shown that these introduce only minor errors in to generate seismic spectra in aspherical earth models. These tech- results (Deuss & Woodhouse 2001). However, it is possible to avoid niques make use of the eigenfunctions and eigenfrequencies of these simplifications, and the full-coupling calculations performed spherically symmetric models. However, it is also possible to com- in this paper are based on solving eq. (A9) exactly. To achieve this, pute the normal modes of general, aspherical models. One route to we adopt an iterative, preconditioned method described in Al-Attar doing so is to adopt the same general strategy as we have already out- et al. (2012), building on earlier work by Hara et al. (1993)and lined: assume that the aspherical eigenfunction can be expressed in Deuss & Woodhouse (2004). terms of the modes of a spherically symmetric reference model, and then solve for the necessary coefficients. This results in essentially the same analysis as in Appendix A3.1, except with time-invariant A3.2 Self- and group-coupling u , and no force term. One may then regard eq. (A9) as a nonlinear eigenvalue problem, which can (in principle) be solved to yield the Although solving eq. (A9) is conceptually straightforward, com- necessary expansion coefficients. putational costs grow rapidly as ω increases. Broadly speaking, The procedure is computationally challenging, limiting its range the number of modes having ω ≤ ω (and thus, the dimension k c of application, especially in early studies. To simplify matters, of the system in eq. (A9)) scales proportionally with ω ;the cost we can make the same assumption as in self-coupling: that of obtaining u˜ then grows as ω . Computational costs rapidly be- ∗ 3 integrals of the form S (x) · A(ω)S (x)d x are only non-zero come prohibitive, even with modern computational resources; for i when the singlets S and S correspond to the same multiplet. This i j early studies, full-coupling was intractable. In consequence, it was allows the eigenvalue problem of eq. (A9) to be decomposed into a necessary to identify strategies for simplifying the computational sequence of independent calculations, one for each multiplet. task. Furthermore, a perturbation-theory approach may be used to Looking again at eq. (A7), we can consider how the equations approximate the nonlinear eigenvalue problem as one in standard behave in the absence of any aspherical structure. In this case, all form. Given the orthogonality relationship, eq. (A3), we can write the integrals vanish unless j = k and thus, in the limit of a purely P = I + P ,where I is an identity matrix and P is the part arising 1 1 spherical model, the matrix M has non-zero elements only along its only from the aspherical density structure (i.e. ρ − ρ¯ ). We then as- diagonal. This leads directly to our expression for seismograms in sume that M(ω) can be expanded around a reference eigenfrequency a spherically symmetric earth, eq. (A4), as is expected. However, it of the multiplet, ω , such that (cf. eq. A9) also motivates the assumption—borne out by empirical evidence— that for typical weakly aspherical models, M remains diagonally (κκ) (κκ) 2 (κκ) 2 M(ω) ≈ V (ω ) + ı ω W − ω P − ω I (A11) 0 0 1 dominant. Elements lying far from the diagonal are, in general, where the superscript (κκ) is used to emphasize that each multi- relatively small—and in group- and self-coupling, they are therefore plet is being treated independently. This approximation relies on ignored. the assumption that |ω − ω | is small, and that the earth model is In eq. (A9), a general element of M is defined dominated by spherical structure, allowing various additional terms ∗ 2 3 involving P , W and ∂V/∂ω can be neglected. It is common to in- M (ω) = S (x) · V + ı ωW − ω ρ(x) S (x)d x. (10a) jk k 1 V ˜ troduce a ‘splitting matrix’, H = V(ω ) + ı ω W − ω P , allowing 0 0 1 eq. (A9) to be written in the form In self-coupling, this is modified by introducing the assertion that (κκ) 2 (κκ) H − ω I u˜ = q˜ . (A12) q = q ⎨ j k M = M = 0 unless l = l (10b) jk (q ,l ,n )(q ,l ,n ) j k j j j k k k Thus, the contribution of each multiplet to spectra—within the n = n . j k self-coupling approximation—is then found by diagonalizing the (2l + 1) × (2l + 1) matrix H.The(2l + 1) constituent singlets within In other words, self-coupling computes integrals using only the the multiplet each contribute at a slightly different eigenfrequency, constituent singlets of each multiplet, and assumes that interactions removing the degeneracy seen in spherically symmetric models. between different multiplets can be entirely ignored. Thus, the spectral contribution of each multiplet becomes a set of In group-coupling, the multiplets are assigned to groups. We (2l + 1) closely spaced—but distinct—peaks. This phenomenon introduce κ to denote a given multiplet, κ = {(q, l, m, n): m ∈ [−l, l]}. We then introduce groups, G, containing one or more is known as ‘splitting’. In practice, the resolution of spectral mea- (i) (i) multiplets: G ={κ ,κ ,...}. In general, all multiplets within a surements may not allow individual singlets to be distinguished, but 1 2 group will have similar eigenfrequency, although individual authors the overall shape of the peak corresponding to each multiplet de- may differ over the particular groupings adopted. Group-coupling pends upon the splitting, and thus upon the detailed properties of the then replaces eq. (A10b) by the assertion that Earth model. In the time domain, the contribution of the multiplet to seismograms can be written in the form (q , l , m , n ) ∈ G j j j j M = 0 unless (10c) jk (κκ) (κκ) (q , l , m , n ) ∈ G s (x, t) = r(x) · exp ıH t · v (A13) k k k k Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 76 F. Akbarashrafi et al. where v and r are vectors that may be computed for the source and ing the reference frequency ω becomes more challenging, as receiver, respectively; this result is due to Woodhouse & Girnius the various multiplets within the group have different eigenfre- (1982). quencies. Typically, some sort of average of these eigenfrequen- The discussion in this section has focused on self-coupling in- cies is used (for a recent discussion on specific averages, see volving only individual isolated multiplets. Similar results can Deuss & Woodhouse 2001). We note that the requirement for be obtained for group-coupling–assumptions: in effect, group- |ω − ω | to remain small acts to restrict the size of the cou- coupling can be regarded as self-coupling of a ‘super-multiplet’, pling group: only multiplets with similar eigenfrequency may be allowing our analysis to translate straightforwardly. However, defin- used. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

Exact free oscillation spectra, splitting functions and the resolvability of Earth's density structure

Loading next page...
1
 
/lp/ou_press/exact-free-oscillation-spectra-splitting-functions-and-the-hN0u9zr3yj

References (69)

Publisher
Oxford University Press
Copyright
Copyright © 2022 The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
DOI
10.1093/gji/ggx539
Publisher site
See Article on Publisher Site

Abstract

Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 Geophysical Journal International Geophys. J. Int. (2018) 213, 58–76 doi: 10.1093/gji/ggx539 Advance Access publication 2017 December 18 GJI Seismology Exact free oscillation spectra, splitting functions and the resolvability of Earth’s density structure 1 2 1 1 3 F. Akbarashrafi, D. Al-Attar, A. Deuss, J. Trampert and A.P. Valentine Department of Earth Sciences, Universiteit Utrecht, Heidelberglaan 2, 3584CS Utrecht, The Netherlands Bullard Laboratories, Department of Earth Sciences, University of Cambridge, Cambridge CB3 0EZ,UK Research School of Earth Sciences, The Australian National University, Canberra, ACT 0200, Australia Accepted 2017 December 15. Received 2017 December 7; in original form 2017 August 23 SUMMARY Seismic free oscillations, or normal modes, provide a convenient tool to calculate low- frequency seismograms in heterogeneous Earth models. A procedure called ‘full mode cou- pling’ allows the seismic response of the Earth to be computed. However, in order to be theoretically exact, such calculations must involve an infinite set of modes. In practice, only a finite subset of modes can be used, introducing an error into the seismograms. By systematically increasing the number of modes beyond the highest frequency of interest in the seismograms, we investigate the convergence of full-coupling calculations. As a rule-of-thumb, it is nec- essary to couple modes 1–2 mHz above the highest frequency of interest, although results depend upon the details of the Earth model. This is significantly higher than has previously been assumed. Observations of free oscillations also provide important constraints on the het- erogeneous structure of the Earth. Historically, this inference problem has been addressed by the measurement and interpretation of splitting functions. These can be seen as secondary data extracted from low frequency seismograms. The measurement step necessitates the calculation of synthetic seismograms, but current implementations rely on approximations referred to as self- or group-coupling and do not use fully accurate seismograms. We therefore also investi- gate whether a systematic error might be present in currently published splitting functions. We find no evidence for any systematic bias, but published uncertainties must be doubled to prop- erly account for the errors due to theoretical omissions and regularization in the measurement process. Correspondingly, uncertainties in results derived from splitting functions must also be increased. As is well known, density has only a weak signal in low-frequency seismograms. Our results suggest this signal is of similar scale to the true uncertainties associated with currently published splitting functions. Thus, it seems that great care must be taken in any attempt to robustly infer details of Earth’s density structure using current splitting functions. Key words: Computational seismology; Seismic tomography; Surface waves and free oscil- lations; Theoretical seismology. et al. 2017). However, as this paper will show, computational limita- 1 INTRODUCTION tions have led most such studies to rely upon seismological approx- Our understanding of the Earth’s large-scale interior structure and imations that might be inadequate for imaging density variations. dynamics draws heavily on observations of seismic free oscillations As such, their conclusions must be approached with considerable (‘normal modes’). In particular, measurements of free-oscillation caution. spectra at long periods represent one of the few classes of geo- In order to derive information about Earth’s interior structure physical observable with meaningful sensitivity to lateral variations from seismic data, we must search for models which yield synthetic in density within the Earth. As such, free oscillation spectra have (predicted) observables which agree with those actually measured. played a central role in attempts to determine the nature of the Many strategies exist for performing this search, but all rely on the ‘Large Low Shear Velocity Provinces’ that have been identified presumption that synthetic observables are computed accurately. within the lowermost mantle, and in attempts to identify the relative The theoretical basis for computing synthetic normal mode spectra importance of thermal and chemical effects as driving forces for is well-developed, relying on a technique known as ‘normal mode mantle convection (e.g. Ishii & Tromp 1999; Trampert et al. 2004; coupling theory’ (e.g. Dahlen 1968, 1969; Woodhouse & Dahlen Lay 2007;Mosca et al. 2012; Moulik & Ekstrom 2016; Koelemeijer 1978; Woodhouse 1980; Woodhouse & Girnius 1982;Park 1986, The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 Resolvability of Earth’s 3-D Density 59 1987, 1990; Romanowicz 1987; Tromp & Dahlen 1990; Lognonne´ foregone conclusion—the validity of any approximation is entirely 1991; Hara et al. 1991, 1993; Um & Dahlen 1992; Deuss & context-specific—but the question must clearly be addressed. The Woodhouse 2001; Al-Attar et al. 2012;Yang&Tromp 2015;Lau first to address this question were Resovsky & Ritzwoller (1998). et al. 2015). Unlike some earlier studies, our calculations are based Being limited in computational resources, they estimated the the- on the work of Al-Attar et al. (2012) and take the full effect of 3-D oretical error by a transfer function. Having at our disposal an density variations into account, following the theory developed in arbitrarily precise forward solver, we can, for the first time, quantify Woodhouse (1980). Our calculations also allow us to take atten- this theoretical error properly. uation fully into account, although in the examples shown in this We therefore perform a simple, but definitive test: we com- paper, only the spherically symmetric attenuation profile of PREM pute accurate synthetic seismograms for the earth model S20RTS (Dziewonski & Anderson 1981) is used. Thus, our framework is (Ritsema et al. 1999), making use of full-coupling. We then mea- physically complete, and allows the seismic response of any Earth- sure splitting functions from this data set as if it were observed data. like body to be expressed exactly. However, the resulting expres- We also compute the ‘true’ splitting functions, by calculating the sions involve an infinite series expansion. This must be truncated appropriate radial integrals for the known input model. Comparing in any computational implementation of the theory, introducing a the two, we find noticeable discrepancies. Repeating the experiment truncation error within the synthetic spectra. for synthetic spectra using the self-coupling approximation allows Various common truncation schemes exist, motivated from physi- us to separate the contribution from theoretical approximations and cal principles. Computational costs scale unpleasantly with the num- regularization in the uncertainty of the measured splitting func- ber of terms retained in the series, so early studies had little choice tions. Finally repeating the experiment for models where density but to adopt simplifying strategies. The most limiting approximation perturbations are set to zero, we conclude that splitting functions is known as ‘self-coupling’, while ‘group-coupling’ incorporates are contaminated by a significant approximation error, a similar additional terms (we explain the distinction in more detail below). regularization error and contain a density signal of the level of the However, a number of studies, starting with Deuss & Woodhouse total measurement uncertainty. (2001), have demonstrated that neither can be relied upon within the frequency range of interest. In particular, Al-Attar et al. (2012) highlighted that use of either the self- or group-coupling approxi- 2 THEORETICAL BACKGROUND mations when computing synthetic spectra leads to errors of similar Imaging the Earth’s interior involves searching for the model param- magnitude to the signal likely to be attributable to lateral variations eters which can match observed seismic data. We therefore require in density within the Earth. Subsequently, a comprehensive study by a forward model, or a way to compute the observations that we Yang & Tromp (2015) concluded that self-coupling is ‘marginally would expect to see if the earth had any given structure. Ideally, this acceptable’ when considering spectra below around 1.5 mHz, and forward model needs to encapsulate the full physics of wave propa- ‘unacceptable’ above this point. Group-coupling was found to also gation. If it does not, the seismograms predicted for each candidate yield significant errors, albeit at a lower level than in self-coupling. structure will be inaccurate, and the results of the imaging process In principle, one could envisage improving group-coupling strate- may be biased by this systematic error (for a full discussion, see gies by computing coupling strengths between modes—perhaps Valentine & Trampert 2016). based on the work of Park (1987)—and ensuring that all signifi- cant interactions are accounted for. However, the notion of coupling strength is itself based on single-scattering approximations, and it 2.1 The forward problem is not clear how accurate the results would be. The optimal approach is therefore to adopt a strategy known as ‘Normal mode seismology’ represents one framework for con- ‘full-coupling’. Despite the name, this continues to involve a trun- structing a seismological forward model, particularly suited to low- cation of the infinite series, but whereas self- and group-coupling frequency, global-scale simulations. The theory incorporates all im- are based around the presumption that the series is dominated by portant physical effects, and can be used to produce seismograms only a few terms, full-coupling aims to retain all terms that have that faithfully reproduce physical reality. However, the computa- any potential to contribute significantly; again, a more precise def- tional resources required are significant, rendering full calculations inition will be given in due course. Building on the existing body intractable for early studies. As a result, a variety of simplifications of work, the first aim of the present paper is to identify more pre- were developed. These can greatly reduce the computational costs, cisely the conditions under which full-coupling can be regarded as at the expense of a degradation in accuracy. A brief, self-contained ‘sufficiently accurate’. account of the underlying theory is given in Appendix, emphasizing This information then permits us to perform high-quality syn- how and why these various approximations can be introduced. Here, thetic experiments to address the second question: can robust in- we sketch some key points. formation about earth structure be inferred from measurements of ‘splitting functions’ (Woodhouse & Giardini 1985; Giardini et al. 1987, 1988)? These are a particular representation of information 2.1.1 Reference eigenfunctions derived from seismic spectra, and they are an important ingredi- ent in many current models of Earth structure (e.g. Ritsema et al. For any spherically symmetric reference model (e.g. PREM, 1999; Ishii & Tromp 1999; Trampert et al. 2004; Ritsema et al. Dziewonski & Anderson 1981), it is relatively straightforward 2011;Mosca et al. 2012; Moulik & Ekstrom ¨ 2016; Koelemeijer to compute eigenfunctions S and their associated eigenfrequen- et al. 2017). However, the theory that connects splitting functions to cies ω ; each pair (S ,ω ) represents one ‘normal mode’ (e.g. k k k structural parameters relies upon precisely the same assumptions as Woodhouse 1988). Conventionally, four indices are used to describe are inherent to self- and group-coupling. Given that these have been the modes: k → (q, l, m, n), with q denoting the ‘class’ (principally demonstrated to introduce significant bias into synthetic spectra, ‘spheroidal’, q = S, or ‘toroidal’, q = T ); l, the angular degree, and are splitting function inferences meaningful? The answer is not a m, the azimuthal order, identifying the spherical harmonic Y (θ, φ) lm Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 60 F. Akbarashrafi et al. which describes the lateral form of S ;and n, the overtone number, Motivated by this, the ‘self-coupling approximation’ modifies the identifying its radial behaviour. definition of M, introducing an over-riding assertion that The eigenfrequency can be shown to depend only upon q, l and n. M (ω) = 0if S and S come from different multiplets. (3b) ij i j For each (q, l, n) combination, there are therefore 2l + 1 eigenfunc- tions that have a different spatial form, but share a common eigen- This results in M taking a block-diagonal form, with each multiplet frequency. This degenerate set is often referred to as a ‘multiplet’, q contributing a block of dimension (2l + 1) × (2l + 1). Solution n l and identified by a label in the form q (e.g. S ). The constituent n l 0 2 of eq. (2) is then much easier, as it decouples into separate calcula- eigenfunctions of each multiplet are then referred to as ‘singlets’. tions for each multiplet. Of course, the resulting solution is only an approximation to the one that would be obtained if the full-coupling matrix were used. 2.1.2 Computing synthetic seismograms The ‘group-coupling’ approximation allows for limited interac- tion between multiplets that are adjacent in frequency. The spectrum A central theme of normal mode seismology is that these eigenfunc- is divided into ‘groups’, with each group containing one or more tions, obtained for a spherically symmetric reference model, may be used as a convenient set of basis functions for representing vector multiplets. Studies may differ in the groupings used. Then, eq. (3b) fields within the (laterally heterogeneous) Earth. In particular, the is relaxed slightly, so that instead seismic waveforms s(x, t) can be expressed in the form (cf. eq. A5) M (ω) = 0 unless S and S belong to the same group. (3c) ij i j s(x, t) = u (t)S (x)(1) k k Again, this simplifies the solution of eq. (2), allowing separate computations for each group. Studies have shown that the resulting where S is a vector field representing an eigenfunction of the seismograms are more accurate than those of self-coupling, but reference model. errors remain significant (e.g. Deuss & Woodhouse 2001; Al-Attar There are an infinite number of eigenfunctions S ,and so in et al. 2012; Yang & Tromp 2015). principle this representation of the seismic wavefield requires an infinite summation. However, if seismograms are to be filtered so that they contain no component above a given frequency, ω ,it max 2.1.5 Normal modes of an aspherical earth model is justifiable to neglect modes with eigenfrequency ω  ω .To k max So far, we have outlined how to obtain seismograms in hetero- make this concrete, we define a cut-off frequency, ω ≥ ω ,and c max geneous earth models. It is also possible to compute the normal work with the finite set of modes having ω ≤ ω . Then, the Fourier- k c modes of a general (non-spherically symmetric) earth model (see transformed expansion coefficients u˜ (ω) may be found by solving Appendix A4) and again, group- or self-coupling approximations a system of linear equations (cf. eq. A9) may be adopted if required. We find that the modes are no longer M(ω)u˜ = q˜ (2) degenerate: each singlet within a multiplet has a distinct eigenfre- quency. This property is sometimes referred to as ‘splitting’ of the Here, M(ω) is a frequency-dependent ‘coupling matrix’ which can multiplet, relative to the spherically symmetric reference. Given that be computed for any earth model presenting 3-D variations in den- realistic earth models contain only weak lateral heterogeneity, the sity, elastic constants and attenuation, while the vector q˜ is a repre- typical frequency shifts involved are small, and due to finite-length sentation of the seismic source (typically an earthquake). time series, we do not generally expect to be able to resolve isolated singlets in spectral observations. Nevertheless, aspherical structure 2.1.3 Full-coupling manifests itself in the overall shape of each ‘peak’ in the spectrum of a seismogram. The term ‘full-coupling’ implies that we solve this linear system (i.e. eq. 2) directly. If we do this, we obtain seismograms that are essentially physically complete. They incorporate only one funda- 2.2 From spectra to structure: the splitting functions mental approximation: truncation of the infinite series at frequency Given the ability to compute synthetic spectra, one may contemplate ω . The first goal of this paper is to quantify the effect of this trunca- formulating an inverse problem that aims to infer earth structure tion upon synthetic seismic spectra, and to assess where ω should from observations. This could be implemented in a variety of ways, lie in order for the truncation to have no appreciable impact upon each with pros and cons. One popular approach involves the use of calculations for frequencies below ω . max ‘splitting functions’. These were introduced within Woodhouse & Giardini (1985), and developed in detail by Giardini et al. (1987, 2.1.4 The self- and group-coupling approximations 1988). These original papers are developed under the self-coupling assumption, where each multiplet is regarded as entirely isolated As ω increases, the dimension of the matrix M increases, and thus from every other. This was recognized to be a simplification of the cost of solving eq. (2) grows rapidly. The elements of M are the underlying physics, and therefore an approximation, but it was found by computing integrals of the form (cf. eq. A10) necessary if calculations were to be feasible using the computational ∗ 3 resources then available. M (ω) = S (x) · A(ω)S (x)d x (3a) ij j Making this approximation, the contribution of a given mul- where A(ω) is an integro-differential operator that depends upon tiplet κ = q within seismograms is entirely characterized by the n l frequency, and also upon the earth model. The elements tend to be (2l + 1) × (2l + 1) diagonal sub-block in the coupling matrix arising (κκ) relatively small unless S and S have similar eigenfrequencies. from eqs. (3a & 3b), which we denote M . Similarly, the multi- i j This leads to structure within M, linked to the degeneracy within plet can be directly identified with a specific portion of an observed each multiplet; it is dominated by the elements lying close to the seismic spectrum. Thus, given a suite of observations of the rele- diagonal. vant portion at different locations and from various seismic events, Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 Resolvability of Earth’s 3-D Density 61 one can formulate an inverse problem to estimate the elements of specifies P-wave speed (α), S-wave speed (β), and density (ρ)per- (κκ) M . Splitting functions provide a mechanism for parametrizing turbations at every point with respect to a spherically symmetric this inversion, designed to provide a straightforward link between reference model. Furthermore, we assume that the model is lat- (κκ) M and the underlying earth model. erally parametrized in terms of spherical harmonics to maximum degree L. Thus, the P-wave speed field may be expressed 2.2.1 Splitting functions for isolated groups of multiplets α(r,θ,φ) = α (r)Y (θ, φ)(8) lm lm l=0 The splitting function framework was extended by Resovsky & Ritzwoller (1995) and others, allowing for cross-coupling between with identical expressions for β and ρ. By deriving the exact ex- (κκ ) specific pairs (or small groups) of multiplets with similar frequen- pression for H (e.g. Woodhouse 1980), and then comparing this mm cies as in group-coupling. For simplicity, we describe here the case with eq. (5), it is possible to identify formulae for the kernels K where two multiplets, κ = q and κ = q , are allowed to interact n l n such that with one another, but are assumed to be isolated from the remainder (κκ ) (α;κκ ) (β;κκ ) of the spectrum. Extension to larger groups, or restriction to the c = α (r)K (r) + β (r)K (r) st st st st st self-coupling case, follows straightforwardly. (ρ;κκ ) + ρ (r)K (r)dr. (9) The isolated-group assumption implies that the only relevant part st st of M(ω) is the symmetric diagonal block If the aspherical model is specified in terms of a different set of (κκ) (κκ ) parameters, equivalent expressions can be found. Thus, an earth M M . (4) (κ κ) (κ κ ) model can be obtained by performing a linear inversion upon mea- M M sured splitting coefficients. A particular attraction of this procedure As discussed in Appendix A4, some additional assumptions allow us is the possibility of selecting multiplets that are presumed to have (κκ ) (κκ ) 2 to introduce the splitting matrix, H,suchthat M ≈ H − ω I. sensitivity at particular depths, and hence to produce models fo- We then parametrize this by introducing a set of complex-valued cused on a given region of the Earth (e.g. Koelemeijer et al. 2013). ‘splitting coefficients’ c such that st l+l s (κκ ) (κκ ) (κκ ) 2.2.2 Splitting functions in the real Earth 2 mm t H = ω δ + ω W + γ  c (5) κκ 0  st mm 0 mm ll s t =−s s=l−l All the foregoing discussion is self-consistent, subject to the fun- damental assumption that the (groups of) multiplets are spectrally where the γ are defined in terms of spherical harmonic triple prod- isolated. However, at least in the context of forward-modelling, it ucts, has been conclusively shown that this approach is insufficiently 2π π accurate for modern seismology, and instead high-quality seismic mm t ∗ γ = Y (θ, φ)Y  (θ, φ)Y (θ, φ)sin θ dθdφ. (6) l m st ll s lm calculations must account for interactions throughout the spectrum. 0 0 Various studies have shown that allowing for cross-coupling in split- This integral can be expressed straightforwardly in terms of the ting function analysis improves results (e.g. Resovsky & Ritzwoller Wigner-3j symbols, and vanishes for many combinations of spheri- 1995, 1998; Irving et al. 2008; Deuss et al. 2013). However, it is cal harmonics. Together, the coefficients describe a ‘splitting func- not possible to extend the splitting function approach to ever-wider tion’ coupling groups: the linearization described in Appendix A4 re- l+l s quires the coupling band-width to be small. If, instead, we choose (κκ ) (κκ ) η (θ, φ) = c Y (θ, φ)(7) st st to work with the full frequency-dependent coupling matrix M(ω), s=l−l t =−s it is unclear whether a splitting-function–style parametrization can which may be regarded as a map depicting spatial variations in spec- be usefully introduced. Where does this leave splitting-function tral properties. Given the symmetry of M, three splitting functions studies? are required to fully characterize our isolated pair of multiplets: In effect, splitting functions should be regarded as a derived (κκ) (κ κ ) (κκ ) η , η and η . data type or ‘secondary observables’ (e.g. Cara & Leveque 1987), ´ ˆ Following Woodhouse & Girnius (1982), the contribution of the intended to convey spectral information in a more manageable form. two multiplets within each seismogram can be computed via diag- The measurement procedure should therefore be seen as a data- onalization of eq. (4) (see eq. A13 and the associated discussion). processing operation, or transformation applied to raw seismic data. Thus, it is possible to pose the inverse problem whereby the c are The particular form of this transformation may be motivated by st determined from observed data; typically, the fiducial frequency ω an out-dated assumption, but it remains a valid transformation to and parameters describing attenuation (which we will not discuss adopt. As such, published measurements of splitting functions may further) are also refined within this process. The inversion is found to be taken at face value: they need not be regarded as somehow be weakly non-linear, and is typically tackled using a Bayesian for- ‘approximate’. The splitting function is defined to be the result of a mulation of the least-squares algorithm (as per Tarantola & Valette specified measurement procedure. 1982). Invariably, the data is insufficient to uniquely constrain a so- However, if we are to adopt the isolated-multiplet assumption in lution, and it is therefore necessary to introduce prior information the measurement procedure, eq. (9) provides only an approximate (i.e. regularization). Discussion of the extent to which this influ- relationship between the measured splitting functions and the un- ences results may be found in numerous sources (e.g. Resovsky & derlying earth structure. Alternatively, one may take the view that Ritzwoller 1999; Kuo & Romanowicz 2002; Pachhai et al. 2016). eq. (9) is the defining property of the splitting functions: in this case, In the isolated-group approximation, the splitting coefficients the measurement procedure is deficient, since it does not properly can be straightforwardly related to structural parameters. For the account for physical effects in wave propagation. Regardless of the purposes of illustration, we assume that the aspherical earth model viewpoint adopted, ‘measured’ and ‘predicted’ splitting functions Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 62 F. Akbarashrafi et al. Figure 1. Position of the stations (red stars) and the Bolivian event (beach ball) used in this study. Station GRFO is highlighted and shown by a blue star. differ, due to approximations within the measurement procedure. may vary considerably between different locations and experiments, Thus, the second key goal of this paper is to assess the impact of limiting the generality of this definition. In addition, the earthquake this approximation: can models produced using eq. (9) nevertheless source itself is also not perfectly well known, although it is appar- be interpreted usefully, or are they too severely contaminated by ent from inspection of eq. (2) that this uncertainty can be absorbed systematic errors? into the observational data uncertainty, by adding the corresponding If fully accurate earth models are to be produced from splitting variances. function data, it is necessary to replace eq. (9) and properly account Alternatively, given that we typically wish to fit seismograms to for multiplet interactions throughout the spectrum. Unfortunately, infer Earth structure, another definition of an ‘acceptable’ truncation it is not clear that this can be achieved by simple modifications of error is that the error should be much smaller than the signal we the modal depth functions K(r). It seems likely that the most fruitful wish to infer. This is perhaps more straightforward, since it does course would involve obtaining sensitivity kernels via adjoint tech- not depend upon the vagaries of seismic noise. Of course, the signal niques within a full-coupling framework, allowing the full range of itself must be above the noise level if we are to make meaningful physical effects to be accounted for. However, we do not pursue this inferences. idea further within the current paper. Indeed, given the capacity to We discuss both cases below and show that they can lead to compute exact sensitivity kernels for spectra using adjoint meth- different conclusions depending on the earth model. To quantify ods, one may question whether there is value in continuing to adopt these ideas, we will only consider synthetic data calculated for the two-stage inversion procedure inherent to splitting function various earth models. Of the three major seismic parameters (S- studies. wave speed, P-wave speed & density), density is known to have the least influence on seismograms, and is a parameter of considerable importance for understanding global geodynamics. We therefore 3 TRUNCATION ERRORS IN focus on the density signal when assessing whether truncation errors FULL-COUPLING SEISMOGRAMS are ‘negligible’. First, we consider the problem of computing seismograms using full-coupling. As we have described, a theoretically complete ex- pression for normal mode seismograms requires summation of an 3.2 Observations using S20RTS infinite set of seismograms. To make the computations tractable, We solved eq. (2) exactly, taking all 3-D effects of varying elas- we must truncate this modal expansion (eq. 1), and only include tic constants and density into account, using the method described those modes with eigenfrequencies below some cut-off, ω .This in Al-Attar et al. (2012), and computed the same set of seismo- truncation will inevitably introduce errors into the computed seis- grams several times varying ω . In all cases, we used the Global mograms. How do we choose ω to ensure that synthetic spectra are CMT Catalog source mechanism for the 1994 Bolivian event (event accurate below some frequency ω ? max code: 060994A) and a global distribution of 129 stations based on those within the Global Seismic Network and IRIS/IDA Seis- mic Network (see Fig. 1). As basis functions, we used eigen- 3.1 The scale of acceptable error functions calculated in in the spherically symmetric model PREM In order to address this question, we must somehow define the mag- (Dziewonski & Anderson 1981). The coupling matrix, explicitly de- nitude of error that we consider acceptable. One way to go about fined in eq. (A10a), is constructed for the shear wave speed model this is to consider the noise in the observed seismograms that we S20RTS (Ritsema et al. 1999), which provides values for δ ln V with ultimately wish to fit. To allow meaningful analysis, the trunca- respect to PREM, and where we added scaled compressional wave tion error should be well below that noise. However, noise levels speed perturbations using δ ln V = 0.5 δ ln V and scaled density p s Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 Resolvability of Earth’s 3-D Density 63 −3 x 10 Effect of source errors 7 Effect of omitting density Truncation error at 5mHz Truncation error at 4mHz Truncation error at 3mHz Reference spectrum (coupling up to 6mHz) 0.5 1 1.5 2 2.5 3 Frequency(mHz) Figure 2. Amplitude spectrum of vertical acceleration at station GRFO (black line) calculated by coupling all modes having frequencies less than ω = 6mHz. We refer to this as the ‘reference’ spectrum, s . Also shown are amplitude spectra of the difference between this reference spectrum and those produced using ref lower values of ω (blue lines), which we describe as the ‘truncation error’ for that cut-off. The red line shows the effect of removing any contribution from lateral density variations within the reference spectrum. Note that this is of broadly similar scale to the truncation error at ω = 3 mHz. The green line shows the effect of perturbing the seismic source with three standard deviations on the reported source uncertainties. perturbations using δ ln ρ = 0.3 δ ln V . These scaling relations are Rather than plotting similar figures for every station around the commonly used and are obtained from mineral physics considera- globe, we sought a method to assess the truncation errors more tions (e.g. Karato 1993). On top of the mantle model, we imple- systematically. We defined a measure of the relative spectral misfit mented the simple crustal model from Woodhouse & Dziewonski within any given frequency range (ω , ω ) 1 2 (1984) and included rotation and ellipticity. In the following, we 2 2 |s(ω) − s (ω)| dω ref will consider vertical component spectra. 1 ξ(s,ω ,ω ) =  . (10) 1 2 2 2 We chose ω = 3 mHz and used ω = 3, 4, 5, and 6 mHz, |s (ω)| dω ref max c increasing the number of coupled singlets (toroidals and spheroidals We then plot histograms of ξ evaluated for every station in our combined) from approximately 2000 to 13 000. We treated the data sets (Fig. 3). Again, we compare the reference spectra to those results from ω = 6 mHz as a set of ‘reference’ synthetic spectra, truncated at ω = 3, 4, and 5 mHz, and to those obtained by using s , and looked at amplitude spectra of the differences s − ref (ω =3mHz) ω = 6 mHz, but omitting lateral density. If we wish to ensure s , etc., to provide an indication of the convergence of the full- ref that the truncation error is consistently smaller than the density coupling calculations (Fig. 2). These difference spectra represent signal, we require the corresponding histograms to not overlap. We the truncation error, and are referred to as such in the figure. It is clearly see that a truncation at 3 mHz will generally lead to errors clear that a truncation at ω = ω can still have substantial errors, c max of a similar scale to the signal anticipated from density. Even for but they decrease rapidly as ω is increased. As far as we know, most a truncation at 4 mHz some overlap of the distributions occurs for published studies that adopt a ‘full-coupling’ strategy have used frequencies close to ω . It therefore seems necessary to ensure max ω ≈ ω . c max that the truncation frequency is significantly above the maximum To provide a sense of scale for this truncation error, we repeated frequency of interest (i.e. ω  ω ). These results are broadly in c max the calculation using ω = 6 mHz, but using a version of S20RTS accordance with what we might expect based on the knowledge of where all lateral density heterogeneity has been omitted (i.e. with the coupling-matrix elements. Generally, modes close together in density structure as given by PREM). The amplitude spectrum of frequency interact most strongly, so truncation errors ought to grow s − s is also shown in Fig. 2 and gives some indication of the ρ¯ ref as we approach the cut-off frequency. magnitude of seismic signal that might be attributable to density To quantify the truncation effects further, we defined a cumulative structure. Perhaps surprisingly, we observe that this is similar to misfit for the entire synthetic data set the magnitude of errors arising from truncation at ω = ω .To c max estimate the effect of source uncertainties, we make use of those (i) (i)   2 |s (ω ) − s (ω )| dω i 0 ref reported in the Global CMT Catalogue, although we note that these ψ(ω) = (11) ω (i) |s (ω )| dω i ref are likely to be an under-estimate (Valentine & Trampert 2012). 0 We therefore implement a large source variation by perturbing each where the summation is over all the stations in the data set. We then parameter by three times the corresponding reported standard devi- plot this quantity as a function of frequency for the various cut-off ation, with positive or negative perturbations chosen randomly, and frequencies (Fig. 4), and for the data set where density heterogeneity again computed spectra. Given the other effects seen in Fig. 2,the has been omitted. There are several ways of reading this figure. First, source effect is very modest. We will therefore neglect its contribu- focusing on the curves for the various choices of ω : suppose we tion in the discussion below. want a precision of (say) 0.1 per cent of the total energy in the signal Amplitude Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 64 F. Akbarashrafi et al. (0.2−3)mHz (0.2−2.5)mHz (2.5−3)mHz 60 60 60 Truncation at 3mHz Truncation at 4mHz Truncation at 5mHz 40 40 40 Omitting density 20 20 20 0 0 0 0.00 0.02 0.04 0.06 0.08 0.10 0.00 0.02 0.04 0.06 0.08 0.10 0.00 0.02 0.04 0.06 0.08 0.10 Spectral misfit Figure 3. Histograms showing relative spectral misfit, ξ (eq. 10), across all stations. Again, we compare reference spectra (ω = 6 mHz) to those truncated at lower frequencies (ω = 3, 4, 5 mHz), with all calculations being performed in the model S20RTS (where we incorporated V and density heterogeneity c p as described in the main text). Results from comparing the same reference spectra to those calculated in a version of S20RTS without any lateral density heterogeneity are also shown. Truncation at 3mHz Truncation at 4mHz Truncation at 5mHz Omitting density 0.01 0.001 0.0001 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 Frequency [mHz] Figure 4. Cumulative misfits defined by eq. (11), relative to reference spectra computed in S20RTS, for three different cut-off frequencies and for the omission of lateral density heterogeneity. The horizontal grey line denotes a misfit of 0.1 per cent, as discussed in the main text. (this choice should be guided by knowledge of the uncertainties 3.3 Model-dependence of the truncation error in observed seismic spectra and the seismic source). In this case, We expect results to vary depending on the model, and in partic- truncation at ω = 3 mHz will produces seismograms accurate up ular on the power spectrum of heterogeneity. Coupling effects are to ω = 2 mHz (as shown by a horizontal line on the figure). max somewhat analogous to scattering, and it is well-known that short- For a truncation at ω = 4 mHz, seismograms are accurate up to wavelength heterogeneity promotes multiple-scattering effects. We 3 mHz within the same precision. It appears that we roughly need might reasonably expect to observe similar effects in coupling cal- coupling of 1 mHz higher than the highest frequency of interest culations, with stronger short-wavelength structure necessitating in the seismogram. Alternatively, if we are interested in imaging relatively higher values of ω . S20RTS has modest amplitudes and a Earth’s density structure, we should require truncation errors to be ‘red’ spectrum, dominated by long-wavelength structure. The model at least an order of magnitude smaller than the density signal in the resolution is far below what the parametrization of the model would spectra, and we see roughly that the cut-off frequency should be at allow, ranging from 2000–4000 km horizontally to 250–750 km ver- 4 mHz for any frequency ω ≤ 3 mHz. These values are based on max tically (Ritsema et al. 2004). It is therefore important to investigate an assumption that Earth’s heterogeneity is similar to that present in how this is dictating the findings in the previous subsection. our implementation of S20RTS. The reader may choose her/his own To provide some insight into this question, we generated a ran- desired precision and adapt the conclusions accordingly, by direct dom model using the same parametrization as S20RTS—a spher- reading of Fig. 4. ical harmonic expansion to degree 20 laterally and 21 splines Number of stations Cumulative misfit Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 Resolvability of Earth’s 3-D Density 65 (0.2−3)mHz (0.2−2.5)mHz (2.5−3)mHz 60 60 60 Truncation at 3mHz Truncation at 4mHz Truncation at 5mHz 40 40 40 Omitting density 20 20 20 0 0 0 0.00 0.02 0.04 0.06 0.00 0.02 0.04 0.06 0.00 0.02 0.04 0.06 Spectral misfit Figure 5. Histograms of relative spectral misfit across all stations defined by eq. (10) for a random model. Truncation at 3mHz Truncation at 4mHz Truncation at 5mHz Omitting density 0.01 0.001 0.0001 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 Frequency [mHz] Figure 6. Cumulative misfits defined by eq. (11) for a randomly generated, spectrally white earth model. Compare Fig. 4. vertically. For δ ln V we drew random numbers uniformly in the Clearly the relation between ω , ω and the earth model is s c max range [−0.01, 0.01] giving an rms-amplitude of about 12 per cent not straightforward. Looking at cumulative misfits (Fig. 6), we see at all depths. We continue to assume that V and ρ are scaled ver- that white models appear to yield larger truncation errors than their sions of V , as with S20RTS. The rms-amplitude of this random red counterparts, as expected. Again, considering the data sets that model is between 5 and 38 times stronger than S20RTS depending include density heterogeneity, we can imagine seeking a precision on depth. This produces a model of true degree-20 lateral resolu- of 0.1 per cent of the total energy in the signal. A truncation at tion, a white spectrum, and much higher vertical resolution. We re- ω = 3 mHz now only produces seismograms accurate up to 1 mHz. peated all the full-coupling calculations already described, using this For a truncation at ω = 4 mHz, seismograms are accurate up to model. around 2 mHz within the same precision. It appears that we now need Again, we produce histograms of relative spectral misfit (Fig. 5) coupling to be at least 2 mHz higher than the highest frequency of and cumulative misfit (Fig. 6). From Fig. 5, we see that now the interest in the seismograms. If, however, we define the acceptable truncation effects are much stronger and overall we need to couple truncation error as being an order of magnitude smaller than the more modes than was necessary in S20RTS. In other words, more signal arising from density heterogeneity, we now see that coupling power in the structural model, and especially at short-wavelengths, up to 5 mHz is required for ω ≤ 3 mHz. Again it is straightforward max necessitates higher values of ω . It is also obvious that the contri- to adapt the conclusions for any other desired precision. bution of density to the seismograms is much lower in the random model than in S20RTS, despite the fact that the rms-amplitude of 3.4 Summary the model is much larger. This is because density imprints itself differently into the seismograms, compared to the model with a red The key conclusions we draw from these experiments may be sum- spectrum. marized: Number of stations Cumulative misfit Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 66 F. Akbarashrafi et al. S S S S S 2 6 2 12 2 13 5 3 0 2 CMB CMB CMB CMB CMB ICB ICB ICB ICB ICB S S S S S 0 5 0 6 0 7 1 4 1 7 CMB CMB CMB CMB CMB ICB ICB ICB ICB ICB S S S 1 8 1 9 1 10 Vs CMB CMB CMB Vp density ICB ICB ICB (κκ) Figure 7. Self-coupling sensitivity kernels K (as in eq. 9) corresponding to the angular order s = 0 splitting coefficient for our selected modes. The solid red line denotes V sensitivity, the solid black line V sensitivity and the dotted black line gives density sensitivity. The modes are ordered somewhat arbitrarily, p s with a general trend from mostly upper-mantle–sensitive to mostly lower-mantle–sensitive. (i) Accurate full-coupling spectra require the cut-off frequency 1999;Masters et al. 2000; Deuss et al. 2011, 2013) still represent ω to be significantly above the maximum frequency present in the the main ‘derived’ data set for long-period seismology. Measured seismograms, ω . To our knowledge, this criterion has not been centre frequencies (i.e. ω in eq. 5) are taken as representative of max 0 met in previously published studies based on full-coupling. the Earth’s spherical average structure (e.g. de Wit et al. 2014)and (ii) Accurate seismograms up to 3 mHz in earth-like models (hav- splitting function coefficients feature as a constraint in many current ing a red spectrum) require truncation at no less than 4 mHz. global tomography models (e.g. Ritsema et al. 2011; Moulik & (iii) Generally truncation errors are stronger at higher frequen- Ekstrom ¨ 2016). Splitting functions have also been used for targeted cies, as might be expected. studies investigating specific aspects of Earth structure, since one (iv) Effects seem more pronounced as models gain more power, can examine the kernels from eq. (9) and identify modes which especially at shorter wavelengths. In models with white spectra, have the desired sensitivity (see Fig. 7). This has particularly been accurate seismograms up to 3 mHz require a truncation level of at exploited to investigate deep-earth structure, including inner core least 5 mHz. studies (beginning with Woodhouse et al. 1986), targeted lower mantle studies (e.g. Koelemeijer et al. 2013, 2017) and density Overall, these observations are all broadly consistent with expec- inferences (e.g. Ishii & Tromp 1999; Resovsky & Trampert 2003). tations based on theoretical arguments. However, the effects are now It is common for splitting function catalogues to publish uncer- quantified, and we have obtained rules-of-thumb that can be applied tainty estimates on the recovered splitting coefficients and centre- to ensure acceptable accuracy when performing calculations. frequencies. In principle, these can be obtained directly from the Bayesian inversion framework (e.g. Tarantola & Valette 1982). However, it is not clear how to define a meaningful prior for the 4 THE INTERPRETABILITY splitting function inversion, and in any case such analysis would OF SPLITTING FUNCTION neglect the inherent non-linearity of the measurement process. This MEASUREMENTS point has previously been made by Resovsky & Ritzwoller (1998). After a careful analysis of all potential contributions to the uncer- We now turn to consider splitting functions, and the extent to which tainty, they estimated combined uncertainties using regressions per- they can provide useful constraints upon Earth structure. Catalogues formed on synthetic data with simulated theoretical errors and data of splitting function measurements (e.g. Giardini et al. 1987;He& uncertainties. More recently, it has been common to use some form Tromp 1996; Resovsky & Ritzwoller 1998; Durek & Romanowicz. Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 Resolvability of Earth’s 3-D Density 67 of boot-strapping procedure to investigate the degree of constraint group-coupling. We used the same event and station distributions upon results (e.g. Deuss et al. 2013). This uncertainty estimate there- and method (forward theory, cut-off frequencies and regularization) fore reflects any apparent inconsistency between different parts of as those in Deuss et al. (2013). The measured splitting coefficients (κκ) the data set: it may not reflect any systematic effects inherent to the c can then directly be compared to the expected coefficients for st measurement procedure or its implementation. model S20RTS (as obtained from eq. 9). Various studies have highlighted that considering interactions Two factors contribute to any difference between measured and between adjacent modes can alter splitting-function results sig- expected coefficients. First, differences may be the result of the ap- nificantly (e.g. Resovsky & Ritzwoller 1995; Irving et al. 2008) proximations in the measurement procedure: this is the effect we and group-coupling, where appropriate, is now the norm for var- wish to address. However, they may also arise as a result of the reg- ious frequency bands (e.g. Deuss et al. 2013). Studies have also ularization applied within the non-linear optimization procedure. investigated the extent to which splitting function inversions are To separate the two contributions, we also calculated S20RTS seis- well-constrained, typically focusing on aspects of the linearized mograms strictly using the self-coupling approximation, and again inversion scheme, such as regularization, parametrization and start- measured splitting functions from these. In this case, any difference ing model (e.g. Resovsky & Ritzwoller 1999; Kuo & Romanowicz between measured and observed coefficients can be attributed solely 2002; Pachhai et al. 2016). However, we are not aware of any previ- to regularization, as data and measurement procedure now contain ous studies that have taken a more holistic approach, to investigate the same physics. Assuming that all contributions are Gaussian dis- whether splitting functions provide a meaningful representation of tributed, the difference between the covariances corresponding to Earth structure that can be interpreted by eq. (9). As already dis- the two sets of splitting function then isolates the theoretical error. cussed, a potential issue arises because the measurement procedure We focused on a selection of modes, which were also measured is inherently based upon isolated-group approximations (either self- by Deuss et al. (2013), namely: S , S , S , S , S , S , S , 0 2 0 5 0 6 0 7 1 4 1 7 1 8 or group-coupling), although their interpretation (via eq. 9) is not. S , S , S , S , S , S . Fig. 7 shows some of the kernels K, 1 9 1 10 2 6 2 12 2 13 5 3 Our results in Section 3, and those of earlier work (e.g. Deuss & which we used in eq. (9) to predict the splitting functions we ought Woodhouse 2001; Al-Attar et al. 2012; Yang & Tromp 2015), all to retrieve if the measurements were perfect. Some of the chosen indicate that isolated-group seismograms provide a poor approx- modes have sensitivity mainly in the lower mantle, some mainly imation to the true physics of wave propagation in the earth. It in the upper mantle, and some are sensitive throughout the mantle. therefore seems distinctly possible that structural models estimated Our discussion will concentrate on coefficients of angular order 2 via isolated-group synthetics may inherit these inaccuracies. and 4, which carry most of the splitting function energy for the We address this by performing a straightforward synthetic test. modes we selected to analyse. This will also avoid the problem of For any earth model, we compute ‘predicted’ self-coupling split- neglected higher degree structure as the measurements are done to ting function coefficients using eq. (9). We also compute accurate higher orders (Resovsky & Ritzwoller 1998). All our synthetic ex- synthetic seismograms, using full-coupling theory and the results of periments are noise-free because we want to isolate the theoretical the previous section. From these seismograms, we then obtain ‘mea- error without a particular noise model (which would be difficult (κκ) sured’ splitting functions c , following the same procedure as is to choose objectively) propagating non-linearly into the results via st used for real data. We then assess the agreement between measured the regularization. Assuming that splitting function measurements and predicted splitting function values. Ideally, they should agree to are secondary data to be taken at face value, we will put our syn- within the stated uncertainty for the measurement procedure. If this thetic measurements into context by using actual uncertainties from is the case, we can conclude that the measured splitting functions are Deuss et al. (2013) obtained by boot-strapping real earthquake data interpretable using eq. (9), and the measurement uncertainties can and thought to represent only the data noise present in the Deuss be propagated through this expression to provide uncertainties on catalogue. structural parameters. If measured and predicted splitting functions differ by more than the measurement uncertainty—but show no ev- idence of any systematic bias—it may be reasonable to continue to 4.2 Observations from synthetic splitting functions apply eq. (9) for interpretation using increased uncertainty levels to reflect the theory error (as per Tarantola & Valette (1982)). How- The coefficients of the splitting function measurements for mode ever, if systematic biases are found, no meaningful interpretation S , a lower mantle mode, are depicted in Fig. 8.Whenself- 0 6 can be made. coupling seismograms are used, the measured coefficients are in- distinguishable from those calculated for S20RTS via the ker- nels, suggesting that the regularization used for this mode is not affecting the measurement. We therefore only show the cal- 4.1 Measurement of synthetic splitting functions culated coefficients in the figure, but not those measured from As stated in the detailed analysis of Resovsky & Ritzwoller (1998), self-coupling seismograms. There is however a difference with those splitting function measurements depend on theoretical errors, regu- measured using the full-coupling spectra as data. To calibrate this larization, noise in the data and event-station distributions. The final difference against the density signal, we also computed predicted splitting function depends non-linearly on all of these. We want to splitting coefficients for a version of S20RTS without lateral density quantify as precisely as possible the theoretical errors, and so we heterogeneity. While there are differences between the three sets of need to specify the other contributions as accurately as possible. splitting coefficients, these are small compared to the magnitude of To do this, all details matter. We therefore generated a synthetic the coefficients. It is therefore instructive to plot differences. full-coupling data set and measured splitting functions following In the bottom panels of Fig. 8 we plot the difference between the approach of Deuss et al. (2013). We considered these new spec- the splitting functions predicted for S20RTS, and those measured tra to be the observed data, referred to as above as s . We will from full-coupling synthetics (black line). Each coefficient is plot- ref (κκ) only focus on the self-coupling splitting function η (eq. 7), al- ted with error-bars, representing the uncertainty reported for this though measurements in certain parts of the spectrum employed mode by Deuss et al. (2013). We see that most coefficients touch Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 Im(C ) Re(C ) Im(C ) Re(C ) Im(C ) Re(C ) Im(C ) Re(C ) Im(C ) Re(C ) Im(C ) Re(C ) 68 F. Akbarashrafi et al. s=4 s=2 0.5 0.0 Measured from full−coupling synthetics Predicted (including density) −2 −0.5 Predicted (no density) 0.4 coupling effect 0.4 density effect 0.2 0.2 0.0 0.0 −0.2 −0.2 −0.4 −0.4 Figure 8. Observed splitting coefficients for mode S for s = 2and s = 4 (top panels). The full black circles represent measurements made on synthetic spectra 0 6 obtained using full-coupling. The red open circles are predicted splitting coefficients calculated using eq. (9), and correspond to those you would expect if the measurements were perfect. The grey diamonds represent splitting coefficients calculated for mantle model S20RTS without lateral density variations. Bottom panels: difference between those measured from full-coupling spectra and predicted coefficients (black circles) and difference between predicted coefficients without and with density (grey diamonds). The vertical bars correspond to error bars inferred by Deuss et al. (2013) for this mode. the zero-line when the error-bars are taken into account (certainly substantially greater than the stated noise level, it remains compara- when 2 standard deviations or the 95 per cent confidence level is ble in size to the theoretical error in splitting function observations. considered), implying that predicted and measured splitting func- Again, it seems unlikely that meaningful density inferences can be tions agree to within the measurement uncertainty. Thus, for this made from this mode using standard methods. These two examples mode, any systematic errors that arise from the use of self-coupling suggest that both density signal, and measurement error, become are negligible compared to the effect of noise within the data. more significant for modes with shallow sensitivity. To see if this We also plot the difference between splitting functions predicted holds in general, we must examine more modes. in S20RTS, and our version without density (grey line). Again, we Comparing the full suite of splitting function coefficients mea- plot error-bars representing the uncertainty reported by Deuss et al. sured from full-coupling synthetics to those predicted using eq. (9), (2013), and find that many of these also touch the zero-line and we find a very strong correlation (with an R-value of 0.99), with certainly within 2 standard deviations. This suggests that the signal no obvious evidence of any systematic trends. Splitting functions from density lies at, or below the noise level of measured splitting depend non-linearly on the spectra, and a more sophisticated error functions. As a result, it is not straightforward to extract meaningful analysis is required before the possibility of any systematic effects information about density from splitting functions for this mode. In can be fully eliminated. However, our results suggest that a reason- addition, we observe that the signal attributable to lateral density able starting point may be to assume that theoretical measurement heterogeneity appears to be similar in magnitude to the difference errors can be modelled by a zero-mean Gaussian distribution, al- between observed and predicted splitting coefficients. Thus, even lowing them to be straightforwardly incorporated into measurement if noise levels could be reduced significantly, it would appear that uncertainties. current splitting function measurement procedures are insufficiently To interpret our observations, it is then most instructive to show accurate to make use of the density signal. the reduced χ misfits for all measured coefficients of a given split- Similar results are seen in Fig. 9, which shows results for the ting function with respect to the coefficients predicted for S20RTS upper mantle mode S . In this example, we see differences which using eq. (9) 2 12 are more significant relative to uncertainties. Errors due to the use pred. meas. c − c of self-coupling in the measurement procedure are now markedly st 1 st χ = (12) greater than the quoted uncertainty, which accounts only for effects N σ st s,t due to noise in the data, even within 2 standard deviations for some coefficients. This is likely to affect any interpretation based on eq. 9. where σ are the uncertainties from Deuss et al. (2013) for the st Although lateral density heterogeneity has a signal that is now also given mode and N represents the number of coefficients in the Cst(μHz) Cst(μHz) Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 Im(C ) Re(C ) Im(C ) Re(C ) Im(C ) Re(C ) Im(C ) Re(C ) Im(C ) Re(C ) Im(C ) Re(C ) Resolvability of Earth’s 3-D Density 69 s=4 s=2 Measured from full−coupling synthetics Predicted (including density) Predicted (no density) −2 −4 1.5 coupling effect 1.0 1.0 density effect 0.5 0.5 0.0 0.0 −0.5 −0.5 −1.0 −1.0 Figure 9. As Fig. 8, but for mode S . 2 12 0.1 full−coupling seismograms self−coupling seismograms 0.01 S S S S S S S S S S S S S 2 62 12 2 13 5 30 20 50 60 71 41 71 81 91 10 Figure 10. Reduced χ misfit defined by eq. (12) for measurements on full- (black squares) and self-coupling (open squares) seismograms using all parameters as determined by Deuss et al. (2013). corresponding splitting function. This statistical measure allows re- that can be placed in models derived partly or fully from these sults to be summarized for easy interpretation; χ = 1 indicates splitting function data. that the difference between measured and predicted values is, on As has already been discussed, discrepancies between measured average, of similar size to the uncertainties. First, we consider the and predicted splitting coefficients may arise due to reliance upon splitting function measurements made upon synthetic seismograms approximations in the measurement procedure, or as a consequence computed using full-coupling (see Fig. 10). For most modes, the of regularization. To compare these effects, Fig. 10 also shows 2 2 value of χ is well above one, indicating that observations differ χ for splitting functions measured from synthetic seismograms from predicted values by an amount significantly greater than cur- generated using self-coupling. In this case, seismogram genera- rently published data uncertainties for splitting coefficients. In other tion and splitting function measurement employ a consistent the- words, the Deuss catalogue underestimates the true uncertainty on oretical framework. Thus, any differences here arise solely from splitting coefficients. This in turn may impact upon the confidence regularization. For several modes, we see that χ obtained from Cst( Hz) Reduced χ misfit Cst( Hz) Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 70 F. Akbarashrafi et al. 0.1 δ in V 0.01 δ in V δ in 0.001 S S S S S S S S S S S S S 2 62 12 2 13 5 30 20 50 60 71 41 71 81 91 10 Figure 11. Reduced χ misfit defined by eq. (12) for model S20RTS, except with δ ln V = 0 (black squares), δ ln V = 0 (open squares) and δ ln ρ = 0(grey s p diamonds). This provides an indication of the average strength of the respective fields within the splitting functions for the 13 modes analysed. self-coupling synthetics is similar to that for full-coupling syn- uncertainty in the Deuss catalogue by a factor of 2 to use a more thetics. This suggests that the discrepancy between prediction and representative number for the true measurement uncertainty. observation in these cases is mostly attributable to regularization. It While the quoted numbers reflect the contribution to uncertainty is possible that reducing the damping applied during the measure- in the Deuss catalogue, can they be used for other catalogues? ment procedure would improve the agreement between predictions The answer is no, because regularization errors and the data uncer- and observations. However, this might also be expected to increase tainty translate differently into coefficient uncertainty, depending the sensitivity to data noise, resulting in increased error-bars on the intimately on all the details on the measurement procedure (event- recovered coefficients (Resovsky & Ritzwoller 1998). It is therefore station distribution, regularization, number of iterations, cut-off fre- not obvious that a net benefit would be gained; re-analysis of real quencies, . . . ). Because we do not have this information available data to investigate this rigorously is beyond the scope of the present for other catalogues, we cannot give detailed results for those. We study. can however say Resovsky & Ritzwoller (1998) made a real effort to Assuming that all contributions to the uncertainties are Gaussian quantify the true uncertainties given the computational restrictions distributed, covariances will simply add up, in other words the to- of 30 years ago. tal posterior covariance in actually measured coefficients may be To put these results into perspective, it is instructive to look at po- expressed as C = C + C + C ,where C is the theoretical tential information on the Earth’s structure contained in the splitting tot th reg d th uncertainty we wish to quantify, C is due to regularization in the functions. We have already described how we generate a version of reg regression and C = σ I corresponds to data uncertainty including S20RTS omitting any lateral density heterogeneity, and use this to the influence of the event-station distribution. This decomposition estimate the magnitude of signal attributable to density. Similarly, is similar to that used in Resovsky & Ritzwoller (1998). The uncer- we can generate seismograms in versions of S20RTS that assume tainties quoted in Deuss et al. (2013) were determined using boot- there is no lateral V heterogeneity (having only lateral variations in strapping: randomly eliminating stations, and also whole events. V and ρ), or no lateral V structure (only V and ρ). In each case, we p p s This was done for a fixed regularization. Thus, these uncertainties can compute expected splitting coefficients using eq. (9), and obtain therefore represent the data covariance, C , including the particular χ values as in eq. (12). When we plot these in Fig. 11, we see that event-station distribution used to create the Deuss catalogue. the signal of V in the splitting functions is 2-3 orders of magnitude Our measurements based on self-coupling seismograms (see larger than the observational uncertainties for most modes, while Fig. 10) do not reflect the calculated splitting functions perfectly. that of V is on average one order of magnitude larger and ρ mostly Because these seismograms are noise free, the discrepancies may be of the same order as the uncertainties. This is interesting in itself, a measure of the regularization noise, which can be expressed as a as it shows that most modes do not contain any significant informa- multiple of the data noise, C = ασ I.Valuesof α can be directly tion on density, regardless of whether we can accurately measure reg read from Fig. 10. In the same figure we also display the result for their splitting function or not. The reason for this might be that measurements on full-coupling seismograms. These represent the wave speeds mostly affect the phase of a seismic spectrum, whereas discrepancies due to regularization and the theoretical error com- density mostly affects its amplitude. When we further consider our bined, which correspond then to C + C = βσ I. Hence β − α is findings that the measured data variances should be multiplied by a th reg a measure of C expressed in units of the covariances in the Deuss factor of 4.3 to represent the total measurement errors, we find that th catalogue. We see that C and C changes for every mode, but as on average the reduced χ will be around 100 for δ ln V , around 4 th reg s a quantitative indication: for the sample of modes we considered, for δ ln V and around 1 for δ ln ρ. on average β − α = 1.3 and α = 2.0. Thus for the Deuss cata- Most studies infer mantle structure using these splitting function logue, the theoretical measurement error, is on average 1.3 times the in a least-squares inversion (e.g. Ritsema et al. 1999; Ishii & Tromp quoted variance and the regularization error is on average 2 times the 1999; Ritsema et al. 2011; Moulik & Ekstrom ¨ 2016; Koelemeijer 2 2 quoted variance, resulting in C = (1 + 1.3 + 2.0)σ I = 4.3σ I. et al. 2017). Given that the V signal is well above the uncertainty tot s d d We therefore recommend users should multiply the published of the splitting functions, shear wave mantle models using such Reduced χ misfit Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 Resolvability of Earth’s 3-D Density 71 data should be fairly robust. However, because the density signal significantly beyond the highest frequency range of interest in the is of the level of the measurement uncertainty, its robustness was seismograms. For seismograms to be accurate up to 3 mHz in typical justifiably questioned in the past (e.g. Resovsky & Ritzwoller 1999; models, coupling to 5 mHz is preferable; 4 mHz may be sufficient Kuo & Romanowicz 2002). Only by moving away from standard for some applications if the Earth has a structure similar to S20RTS. techniques and using full sampling, can density be inferred with Our experimental evidence appears to corroborate theoretical ex- a meaningful uncertainty (Resovsky & Trampert 2003). The latter pectations that truncation errors become larger as short-wavelength study showed that the density model had on average an uncer- structure becomes more dominant. tainty of about 50 per cent of its magnitude, or 100 per cent within Self-coupling splitting functions, often used in the construction a 95 per cent confidence level. This generated lively debates, and of 3-D Earth models, are themselves inherently based on the self- for instance, the question of buoyant versus heavy large low-shear or group-coupling approximations, and may therefore contain sys- wave provinces is still waiting to be settled, although the heavy side tematic measurement errors. By measuring self-coupling splitting (Ishii & Tromp 1999; Resovsky & Trampert 2003; Trampert et al. functions on exact full-coupling seismograms, we showed that the 2004;Mosca et al. 2012; Moulik & Ekstrom ¨ 2016) got a recent recovered splitting function does not match theoretical predictions. boost with an independent study inverting tidal data (Lau et al. There is also a significant regularization contribution to the recov- 2017). Given that the Deuss catalogue uses far more seismograms ered splitting functions. There is no visible apparent bias from any to infer splitting functions than any previous catalogue, and we now of these contributions, but recently published (Deuss et al. 2013) have a good representations of the true measurement uncertainty, standard deviations of data uncertainties for splitting coefficients constructing a mantle model using a sampling algorithm is well need to be increased by a factor of around 2 to represent the total worth repeating. Alternatively, direct spectra inversions using full error within the measurement process. In that case, our results show coupling theory, where density has a higher signal-to-noise ratio, that the information on the 3-D density structure contained in the are becoming computationally feasible. measured splitting functions has the same magnitude as the total un- certainty. The information on the variation of compressional-wave speed is about 4 times that of the uncertainty and the information 4.3 Summary on the shear-wave speed has a signal-to-noise ratio of about 100. To go beyond a statistical treatment of the theoretical error, one Again, our results can be summarized as: route forward may be to simply avoid the use of splitting functions: (i) We confirm that there is a significant theoretical measurement with modern computational resources, it may be feasible to infer error in currently published splitting functions. structural information directly from spectra via full-coupling cal- (ii) There is no obvious bias, but this theoretical error can, to a culations. Alternatively, it may also be possible to continue to use first approximation, be assumed to be random noise. splitting functions, measured using self- or group-coupling approx- (iii) There is also a significant contribution from regularization imations, and instead replace the kernels within eq. (9) by those to the uncertainty in splitting functions. that account for these measuring approximations correctly (e.g. via (iv) Taking the total uncertainty into account, currently published adjoint theory). However, until such studies are undertaken, we splitting functions contain information on V with a signal-to-noise must assume markedly larger uncertainties than reported for mea- ratio of about 100, information on V with a signal-to-noise ratio sured splitting coefficients. Given this, our modelling suggests that of about 4, and information on ρ with a signal-to-noise ratio of Earth’s density signal lies at or below the presumed noise level. In about 1. consequence, there is no prospect of recovering robust estimates (v) Since the density signal is small in currently published split- of density structure. It therefore follows that currently published ting functions, advanced techniques should be employed to infer density models ought to be interpreted with great caution. density structure. (vi) The signal from V is very strong, and its inference is un- ACKNOWLEDGEMENTS likely to be significantly affected by the increased measurement uncertainty in currently published splitting functions. Thus, we do We are grateful for constructive reviews by Philippe Lognonnea ´ nd not believe our results should affect published models of V structure Joseph Resovsky. The research leading to these results has received that may be derived these splitting function measurements. funding from the European Research Council (ERC) under the Eu- ropean Union’s Seventh Framework Programme (FP/2007-2013) These results are specific to the Deuss catalogue (Deuss et al. 2013). grant agreement number 320639 (iGEO) and under the European Using other catalogues requires repeating the work in this study Union’s Horizon 2020 research and innovation programme grant if authors wish to infer density, the exception possibly being the agreement number 681535 (ATUNE). The code used to perform Colorado catalogue (Resovsky & Ritzwoller 1998). We also remark the full-coupling calculations makes extensive use of routines de- that the density signal may be more effectively isolated by working veloped by John Woodhouse, and we thank him for allowing us to directly upon seismic spectra, where the observational uncertainty use them. is expected to be smaller: the results of Section 3 suggest that this should be feasible. REFERENCES 5 CONCLUSIONS Al-Attar, D. & Crawford, O., 2016. Particle relabelling transformations in elastodynamics, Geophys. J. Int., 205, 575–593. Previous studies have shown that self- and group-coupling approxi- Al-Attar, D., Woodhouse, J. & Deuss, A., 2012. Calculation of normal mode mations ought to be regarded as obsolete in the context of computing spectra in laterally heterogeneous earth models using an iterative direct seismograms in 3-D Earth models. By systematically increasing the solution method, Geophys. J. Int., 189, 1038–1046. number of modes in full-coupling calculations, we have demon- Cara, M. & Lev ´ eque, ˆ J., 1987. Waveform inversion using secondary observ- strated that accurate full-coupling requires inclusion of modes ables, Geophys. Res. Lett., 14, 1046–1049. Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 72 F. Akbarashrafi et al. Dahlen, F., 1968. The normal modes of a rotating, elliptical Earth, Geophys. Lognonne, ´ P., 1991. Normal modes and seismograms in an anelastic rotating J. R. astr. Soc., 16, 329–367. Earth, J. geophys. Res., 96, 20 309–20 319. Dahlen, F., 1969. The normal modes of a rotating, elliptical Earth—II. Near- Masters, G., Laske, G. & Gilbert, F., 2000. Matrix autoregressive analysis resonance multiplet coupling, Geophys.J.R.astr.Soc., 18, 397–436. of free-oscillation coupling and splitting, Geophys. J. Int., 143, 478–489. Dahlen, T. & Tromp, J., 1998. Theoretical Global Seismology, Princeton Masters, G., Barmine, M. & Kientz, S., 2011. Mineos User Manual, Com- University Press. putational Infrastructure for Geodynamics, v. 1.0.2. de Wit, R., Kaufl, ¨ P., Valentine, A. & Trampert, J., 2014. Bayesian inversion Mosca, I., Cobden, L., Deuss, A., Ritsema, J. & Trampert, J., 2012. Seis- of free oscillations for Earth’s radial (an)elastic structure, Phys. Earth mic and mineralogical structures of the lower mantle from probabilistic planet. Inter., 237, 1–17. tomography, J. geophys. Res., 117, B06304, doi:10.1029/2011JB008851. Deuss, A. & Woodhouse, J., 2001. Theoretical free-oscillation spectra: the Moulik, P. & Ekstrom, ¨ G., 2016. The relationships between large-scale importance of wide band coupling, Geophys. J. Int., 146, 833–842. variations in shear velocity, density, and compressional velocity in the Deuss, A. & Woodhouse, J., 2004. Iteration method to determine the eigen- Earth’s mantle, J. geophys. Res., 121, 2737–2771. values and eigenvectors of a target multiplet including full mode coupling, Pachhai, S., Tkalci ˇ c, ´ H. & Masters, G., 2016. Estimation of splitting func- Geophys. J. Int., 159, 326–332. tions from Earth’s normal mode spectra using the neighbourhood algo- Deuss, A., Ritsema, J. & van Heijst, H., 2011. Splitting function measure- rithm, Geophys. J. Int., 204, 111–126. ments for Earth’s longest period normal modes using recent large earth- Park, J., 1986. Synthetic seismograms from coupled free oscillations: Effects quakes, Geophys. Res. Lett., 38, L04303, doi:10.1029/2010GL046115. of lateral structure and rotation, J. geophys. Res., 91, 6441–6464. Deuss, A., Ritsema, J. & van Heijst, H., 2013. A new catalogue of normal- Park, J., 1987. Asymptotic coupled-mode expressions for multiplet ampli- mode splitting function measurements up to 10 mHz, Geophys. J. Int., tude anomalies and frequency shifts on an aspherical earth, Geophys. J. 193, 920–937. R. astr. Soc., 90, 129–169. Durek, J.J. & Romanowicz, B., 1999. Inner core anisotropy inferred by direct Park, J., 1990. The subspace projection method for constructing coupled- inversion of normal mode spectra, Geophys. J. Int., 139, 599–622. mode synthetic seismograms, Geophys. J. Int., 101, 111–123. Dziewonski, A. & Anderson, D., 1981. Preliminary reference Earth model, Pekeris, C. & Jarosch, H., 1958. The free oscillations of the earth, in Phys. Earth planet. Inter., 25, 297–356. Contributions in Geophysics in Honor of Beno Butenberg, eds Benioff, Gharti, H. & Tromp, J., 2017. A spectral-infinite-element solution of Pois- H., Ewing, M., Howell, B. & Press, F., Pergamon Press. son’s equation: an application to self gravity,arXiv:1706.00855v1. Resovsky, J. & Ritzwoller, M., 1995. Constraining odd-degree Earth struc- Giardini, D., Li, X.-D. & Woodhouse, J., 1987. Three-dimensional structure ture with coupled free oscillations, Geophys. Res. Lett., 22, 2301–2304. of the Earth from splitting in free-oscillation spectra, Nature, 325, 405– Resovsky, J. & Ritzwoller, M., 1998. New and refined constraints on three- 411. dimensional Earth structure from normal modes below 3 mHz, J. geophys. Giardini, D., Li, X.-D. & Woodhouse, J., 1988. Splitting functions of long- Res., 103, 783–810. period normal modes of the Earth, J. geophys. Res., 93, 13 716–13 742. Resovsky, J. & Ritzwoller, M., 1999. Regularization uncertainty in density Gilbert, F., 1970. Excitation of the normal modes of the Earth by earthquake models estimated from normal mode data, Geophys. Res. Lett., 26, 2319– sources, Geophys. J. R. astr. Soc., 22, 223–226. 2322. Hara, T., Tsuboi, S. & Geller, R., 1991. Inversion for laterally heterogeneous Resovsky, J. & Trampert, J., 2003. Using probabilistic seismic tomogra- earth structure using a laterally heterogeneous starting model: preliminary phy to test mantle velocity-density relationships, Earth planet. Sci. Lett., results, Geophys. J. Int., 104, 523–540. 215,121–234. Hara, T., Tsuboi, S. & Geller, R., 1993. Inversion for laterally heteroge- Ritsema, J., van Heijst, H. & Woodhouse, J., 1999. Complex shear wave neous upper mantle S-wave velocity structure using iterative waveform velocity structure imaged beneath Africa and Iceland, Science, 286(5446), inversion, Geophys. J. Int., 115, 667–698. 1925–1931. He, X. & Tromp, J., 1996. Normal-mode constraints on the structure of the Ritsema, J., van Heijst, H. & Woodhouse, J., 2004. Global transition zone Earth, J. geophys. Res., 101, 20 053–20 082. tomography, J. geophys. Res., 109, doi:10.1029/2003JB002610. Irving, J., Deuss, A. & Andrews, J., 2008. Wide-band coupling of Earth’s Ritsema, J., Deuss, A., van Heijst, H. & Woodhouse, J., 2011. S40RTS: a normal modes due to anisotropic inner core structure, Geophys. J. Int., degree-40 shear-velocity model for the mantle from new Rayleigh wave 174, 919–929. dispersion, teleseismic traveltime and normal-mode splitting function Ishii, M. & Tromp, J., 1999. Normal-mode and free-air gravity constraints measurements, Geophys. J. Int., 184, 1223–1236. on lateral variations in velocity and density of Earth’s mantle, Science, Rogister, Y. & Valette, B., 2009. Influence of liquid core dynamics on 285, 1231–1236. rotational modes, Geophys. J. Int., 176, 368–388. Karato, S.-i., 1993. Importance of anelasticity in the interpretation of seismic Romanowicz, B., 1987. Multiplet-multiplet coupling due to lateral hetero- tomography, Geophys. Res. Lett., 15, 1623–1626. geneity: asymptotic effects on the amplitude and frequency of the Earth’s Koelemeijer, P., Deuss, A. & Ritsema, J., 2013. Observations of core-mantle normal modes, Geophys. J. R. astr. Soc., 90, 75–100. boundary Stoneley modes, Geophys. Res. Lett., 40, 2557–2561. Tarantola, A. & Valette, B., 1982. Generalized nonlinear inverse problems Koelemeijer, P., Deuss, A. & Ritsema, J., 2017. Density structure of Earth’s solved using the least squares criterion, Rev. Geophys. Space Phys., 20, lowermost mantle from Stoneley mode splitting observations, Nat. Com- 219–232. mun., 8, 15241. Trampert, J., Deschamps, F., Resovsky, J. & Yuen, D., 2004. Probabilistic Komatitsch, D. & Tromp, J., 2002. Spectral-element simulations of global tomography maps chemical heterogeneities throughout the lower mantle, seismic wave propagation—I. Validation, Geophys. J. Int., 149, 390–412. Science, 306, 853–856. Kuo, C. & Romanowicz, B., 2002. On the resolution of density anomalies in Tromp, J. & Dahlen, F., 1990. Free oscillations of a spherical anelastic earth, the Earth’s mantle using spectral fitting of normal-mode data, Geophys. Geophys. J. Int., 103, 707–723. J. Int., 150, 162–179. Um, J. & Dahlen, F., 1992. Normal mode multiplet coupling on an aspherical, Lau, H., Yang, H.-Y., Tromp, J., Mitrovica, J., Latychev, K. & Al-Attar, anelastic earth, Geophys. J. Int., 111, 11–34. D., 2015. A normal mode treatment of semi-diurnal body tides on an Valentine, A. & Trampert, J., 2012. Assessing the uncertainties on seismic aspherical, rotating and anelastic Earth, Geophys. J. Int., 202, 1392–1406. source parameters: towards realistic error estimates for centroid-moment- Lau, H., Mitrovica, J., Davies, J., Tromp, J., Yang, H.-Y. & Al-Attar, D., 2017. tensor determinations, Phys. Earth. planet. Inter., 210–211, 36–49. Tidal tomography constrains Earth’s deep-mantle buoyancy, Nature, 551, Valentine, A. & Trampert, J., 2016. The impact of approximations and 321–326. arbitrary choices on geophysical images, Geophys. J. Int., 204, 59–73. Lay, T., 2007. Deep earth structure—lower mantle and D ,in Treatise on Woodhouse, J., 1980. The coupling and attenuation of nearly resonant mul- Geophysics, Vol. 1, pp. 620–654, eds Romanowicz, B. & Dziewonski, ´ A., tiplets in the Earth’s free oscillation spectrum, Geophys.J.R.astr.Soc., Elsevier. 61, 261–283. Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 Resolvability of Earth’s 3-D Density 73 Woodhouse, J., 1988. The calculation of the eigenfrequencies and eigenfunc- algorithms. In a seismological context, one strategy is to adopt a tions of the free oscillations of the earth and the sun, in Seismological basis constructed on a finite-element–style mesh. This results in Algorithms, pp. 321–370, ed. Doornbos, D., Academic Press. computational tools such as SPECFEM3D_GLOBE (e.g Komatitsch & Woodhouse, J. & Dahlen, F., 1978. The effect of a general aspherical pertur- Tromp 2002), where simulation costs scale with the length of time- bation on the free oscillations of the Earth, Geophys. J. R. astr. Soc., 53, series that must be calculated. Normal mode seismology takes a 335–354. different approach: it relies upon basis functions that have global Woodhouse, J. & Deuss, A., 2007. Earth’s free oscillations, in Seismology support, and the resulting algorithms have computational costs that and structure of the Earth, vol. 1 of Treatise on Geophysics, chap. 2, depend largely on the desired frequency range of solution. Both pp. 31–65, eds Romanowicz, B. & Dziewonski, ´ A., Elsevier. approaches are formally complete and can provide faithful simula- Woodhouse, J. & Dziewonski, A., 1984. Mapping the upper mantle: Three- tion of physical reality. However, we note that SPECFEM3D_GLOBE dimensional modeling of Earth structure by inversion of seismic wave- forms, J. geophys. Res., 89, 5953–5986. currently omits self-gravitation effects, as these are challenging to Woodhouse, J. & Giardini, D., 1985. Inversion for the splitting function of incorporate in a spectral element framework, and is therefore not isolated low order normal mode multiplets, EOS, Trans. Am. geophys. suitable for modelling seismograms at the very low frequencies we Un., 66, 300. consider in this paper (although efforts are being made to overcome Woodhouse, J. & Girnius, T., 1982. Surface waves and free oscillations in a this difficulty: see Gharti & Tromp 2017). regionalized earth model, Geophys. J. R. astr. Soc., 68, 653–673. Woodhouse, J., Giardini, D. & Li, X.-D., 1986. Evidence for inner core anisotropy from free oscillations, Geophys. Res. Lett., 13, 1549–1552. Yang, H.-Y. & Tromp, J., 2015. Synthetic free oscillation spectra: an ap- A2 Spherical-earth eigenfunctions praisal of various mode-coupling methods, Geophys. J. Int., 203, 1179– The particular family of basis functions we use are the eigenfunc- tions, or ‘normal modes’, of a spherically symmetric earth model. In this case, the material properties vary only with depth, and eq. (A1) simplifies considerably. The Coriolis term can be discarded, and the APPENDIX: SYNTHETIC SEISMIC operator V takes a simpler form, which we denote by V. Similarly, WAVEFORMS AND SPECTRA—THE we use ρ¯ to highlight that the density distribution is also spherically NORMAL MODE FRAMEWORK symmetric. Neglecting the force term, and working in the frequency A complete description of the theoretical basis for normal mode domain, it can then be shown that the differential equations admit a seismology lies well beyond the scope of this paper. Here, we aim non-trivial solution, S (x), at a discrete set of frequencies ω . Thus, k k to provide sufficient background to allow the reader to appreciate we have the motivations for—and distinctions between—self-, group- and V − ρω ¯ S =0(A2) full-coupling. In particular, we aim to highlight precisely where and k why approximations need to be introduced within the various theo- where S represents the kth normal mode, and ω is its correspond- k k ries, providing insight into potential for systematic errors. Readers ing eigenfrequency. wishing for a more in-depth treatment are directed to a summary pa- This eigenvalue problem can be solved relatively straightfor- per by Woodhouse & Deuss (2007), or the comprehensive account wardly (Pekeris & Jarosch 1958). There are an infinite number contained within Dahlen & Tromp (1998). Of course, our discussion of (ω ,S ) pairs that satisfy eq. (A2), but only a finite set exist k k here is very general: specific implementations of these techniques within any given frequency interval ω ≤ ω ≤ ω .Thus,itisfea- 0 k 1 invariably adopt additional simplifications and restrictions, which sible to compute a complete set of eigenfunctions within the given may impact further upon accuracy. frequency range, and an algorithm described by Woodhouse (1988) provides a means for doing so. This method is implemented, for ex- ample, within the MINEOS software package described by Masters A1 The seismic wave equation et al. (2011). The seismic wavefield within the Earth, s(x, t), can be described by It turns out (see e.g. Pekeris & Jarosch 1958) that distinct families the seismic wave equation. This has the general form of normal modes may be identified. Principally, two classes exist, known as ‘spheroidal’ and ‘toroidal’ modes; in addition, ‘radial’ ∂ ∂ V + W + ρ s = f (A1) modes are sometimes discussed as a separate category, although they ∂t ∂t are essentially a particular subset of the spheroidals. Furthermore, where V represents a linear integro-differential operator that de- distinct sets of toroidal modes may be identified for each isolated pends upon the material properties of the Earth; W is an operator solid region within the structural model, thus two groups of toroidal representing effects due to the Coriolis force; ρ = ρ(x) is the equi- modes may be found for earth models such as PREM, where a solid librium density distribution within the Earth; and f = f(x, t)repre- inner core and mantle/crust region are separated by a fluid outer core. sents the force due to the seismic source. If viscoelastic effects are to The modes take the spatial form of spherical harmonics laterally; be properly accounted for, the operator V is itself a function of time, their radial functions must be found by numerical integration. Thus, and in this case Vs implies a convolution operation. Given specific the index k on S may be recast as a set of four indices, k → (q, realizations of the various quantities, we aim to solve eq. (A1) to l, m, n). Here, q defines the ‘family’ of mode being considered; l obtain synthetic seismic observables. and m are the angular degree and azimuthal order specifying the A general strategy for solving differential equations entails ex- appropriate spherical harmonic Y (θ, φ); and n is the ‘overtone lm panding the solution in terms of a known basis set, thereby reducing number’, which indexes the set of permissible radial functions. the equations to algebraic form. Although the solution remains the For present purposes, the main point of interest is that ω de- same irrespective of the basis chosen, the computational procedure pends only on q, l and n: in a spherically symmetric model, modes required to obtain that solution can vary considerably. Thus, an that differ only by their m-value share the same eigenfrequency. In appropriate choice of basis can be critical to developing efficient the literature of normal mode seismology, it is common to see an Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 74 F. Akbarashrafi et al. individual (q, l, m, n) state described as a ‘singlet’; the (2l + 1) symmetric reference structure. In addition, we implicitly assume degenerate singlets corresponding to a given (q, l, n) triplet are then that all interfaces within the aspherical model remain concentric referred to as a ‘multiplet’. This triplet is often written in the form spheres: in other words, no topography is present, either on the q ,with q being either S (in the case of spheroidal modes), or T free surface or on internal discontinuities (such as the core–mantle n l (toroidals). For brevity, we will continue to make use of the index k boundary). To incorporate such topography, studies have tradition- within this discussion, unless there is reason to do otherwise. ally relied on first-order boundary perturbation theory (Woodhouse & Dahlen 1978; Woodhouse 1980), which is inherently an approx- imation. Recently, an exact treatment of boundary perturbations A2.1 Completeness and orthogonality has been developed by Al-Attar & Crawford (2016), although this theory has yet to be fully implemented. In any case, the effects of The eigenfunctions of a spherically symmetric model can be shown realistic-scale topography are not likely to impact the conclusions to have some important properties, which stem from the self- of this paper, and need not be considered further here. adjointness of the operator V. In particular, the infinite set of S Given that the S form a complete basis, we can express the form a complete basis for vector fields within the earth model, and wavefield s(x, t)inthe form they may be shown to be orthogonal, in the sense that s(x, t) = u (t)S (x), (A5) k k ∗ 3 ρ¯ (x)S (x) · S (x)d x = δ (A3) k jk where the real part is again understood. u represent time-varying where δ represents the Kronecker delta. Here, the integration vol- jk expansion coefficients. Strictly, this is only valid if the summa- ume V is the complete interior of the earth model, and an asterisk tion includes the entire and infinite set of eigenfunctions for the is used to denote complex conjugation. spherically symmetric model. However, we note the similarity with Strictly, these properties are only proven for entirely solid models, eq. (A4), and the fact that realistic aspherical models contain only and we ignore any complications that might arise from the existence relatively weak perturbations to the spherically symmetric reference of ‘undertone modes’ in the Earth’s fluid outer core. These are not model. It therefore seems reasonable to assert that modes with eigen- thought to be problematic in a seismological context, as they exist frequency far in excess of ω are unlikely to contribute within the max at lower frequencies than we typically consider (e.g. Rogister & frequency band of interest. We introduce a cut-off frequency, ω ≥ Valette 2009). ω , and work with the finite set of S having ω ≤ ω . max k k c A2.2 Seismograms in a spherically symmetric model A3.1 Full-coupling Knowledge of the spherical-earth eigenfunctions allows straight- In order to make use of this expansion, eq. (A5), we need to ob- forward synthesis of the wavefield expected within a spherically tain expressions for the coefficients u . Substituting into the wave symmetric body. As shown by Gilbert (1970), and recounted in equation, eq. (A1), we find Woodhouse & Deuss (2007), the three-component wavefield can be reduced to the form (u VS + u˙ WS + ρu¨ S ) = f (A6) k k k k k k s¯(x, t) = E (1 − exp ı ω t)S (x)(A4) k k k ω and therefore, for any eigenfunction S , we obtain k j where the real part is understood and we only consider positive ∗ 3 ∗ 3 u S · VS d x + u˙ S · WS d x k k k k j j times. E represent the modal excitations, and may be obtained by V V calculating certain spatial integrals of the force term f(x, t). The ∗ 3 ∗ 3 process of computing seismograms by this formula is often referred + u¨ ρ S · S d x = S · f d x. (A7) k k j j V V to—for obvious reasons—as ‘mode summation’. In principle, the summation over k includes an infinite number Since the sum only contains a finite number of terms corresponding of modes. However, we see that the dominant contribution of each to frequencies lower than ω , this can equivalently be expressed as mode to the seismogram lies at its eigenfrequency. Thus, if the a matrix equation, seismogram is only required within a finite frequency band, only Vu + W˙u + P¨u = q (A8) modes lying within that band need to be summed, and this makes the summation finite. In practice, we always work with band-limited where the elements of P, W, V and q are obtained by evaluating the seismic data, which is a necessary consequence of digital signal appropriate integrals from eq. (A7). Taking the Fourier transform recording. If the maximum frequency present in our data set is ω , max allows us to eliminate the derivatives, and introduce the matrix it is straightforward—and sufficient—to compute a mode catalogue M(ω), such that containing all modes with ω ≤ ω . k max V(ω) + ı ωW − ω P u˜ = M(ω)u˜ = q˜ (A9) where the frequency-dependence of V arises from the inclusion of A3 Seismograms in an aspherical earth model viscoelasticity within V; strictly, eq. (A8) involves a convolution for We now return to consider general, aspherical earth models, and the first term. The elements of M and q˜ are readily evaluated, and the solution of eq. (A1). Invariably, the aspherical model will be thus obtaining u˜ is a straightforward exercise in linear algebra. specified as a 3-D structure superimposed upon a spherically sym- The ‘full-coupling’ approach involves solving eq. (A9) exactly. metric reference model (often based on the model PREM, as de- In reaching this point, we have introduced only one approximation: scribed in Dziewonski & Anderson 1981). In what follows, we as- we have truncated the modal expansion at frequency ω . To date, sume that the S have been obtained using the relevant spherically there has been little effort to investigate the errors introduced by k Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 Resolvability of Earth’s 3-D Density 75 this truncation, or assess how ω should be chosen. Most of the for a single group, G. Thus, some limited interactions are allowed existing studies that employ full-coupling assume it is sufficient to between multiplets. However, (assuming a sensible ordering of the use ω = ω , but there is little experimental or theoretical justi- modes), both self- and group-coupling cause the matrix M to take c max fication to support this. Thus, one goal of this paper is to explore on a block-diagonal form. Exploiting this structure allows a great how changes to ω impact the seismograms produced, and to iden- reduction in the computational effort required to solve eq. (A9), at tify where the cut-off must be placed to ensure computations are the expense, of course, of a deterioration in accuracy. essentially exact at frequencies below ω . max We note that a number of ‘full-coupling’ studies introduce addi- A4 Normal modes in aspherical models tional approximations within eq. (A9)—similar to those discussed in Appendix A4, below—designed to reduce computational com- Thus far, we have discussed the use of full-, group- and self-coupling plexity. It has been shown that these introduce only minor errors in to generate seismic spectra in aspherical earth models. These tech- results (Deuss & Woodhouse 2001). However, it is possible to avoid niques make use of the eigenfunctions and eigenfrequencies of these simplifications, and the full-coupling calculations performed spherically symmetric models. However, it is also possible to com- in this paper are based on solving eq. (A9) exactly. To achieve this, pute the normal modes of general, aspherical models. One route to we adopt an iterative, preconditioned method described in Al-Attar doing so is to adopt the same general strategy as we have already out- et al. (2012), building on earlier work by Hara et al. (1993)and lined: assume that the aspherical eigenfunction can be expressed in Deuss & Woodhouse (2004). terms of the modes of a spherically symmetric reference model, and then solve for the necessary coefficients. This results in essentially the same analysis as in Appendix A3.1, except with time-invariant A3.2 Self- and group-coupling u , and no force term. One may then regard eq. (A9) as a nonlinear eigenvalue problem, which can (in principle) be solved to yield the Although solving eq. (A9) is conceptually straightforward, com- necessary expansion coefficients. putational costs grow rapidly as ω increases. Broadly speaking, The procedure is computationally challenging, limiting its range the number of modes having ω ≤ ω (and thus, the dimension k c of application, especially in early studies. To simplify matters, of the system in eq. (A9)) scales proportionally with ω ;the cost we can make the same assumption as in self-coupling: that of obtaining u˜ then grows as ω . Computational costs rapidly be- ∗ 3 integrals of the form S (x) · A(ω)S (x)d x are only non-zero come prohibitive, even with modern computational resources; for i when the singlets S and S correspond to the same multiplet. This i j early studies, full-coupling was intractable. In consequence, it was allows the eigenvalue problem of eq. (A9) to be decomposed into a necessary to identify strategies for simplifying the computational sequence of independent calculations, one for each multiplet. task. Furthermore, a perturbation-theory approach may be used to Looking again at eq. (A7), we can consider how the equations approximate the nonlinear eigenvalue problem as one in standard behave in the absence of any aspherical structure. In this case, all form. Given the orthogonality relationship, eq. (A3), we can write the integrals vanish unless j = k and thus, in the limit of a purely P = I + P ,where I is an identity matrix and P is the part arising 1 1 spherical model, the matrix M has non-zero elements only along its only from the aspherical density structure (i.e. ρ − ρ¯ ). We then as- diagonal. This leads directly to our expression for seismograms in sume that M(ω) can be expanded around a reference eigenfrequency a spherically symmetric earth, eq. (A4), as is expected. However, it of the multiplet, ω , such that (cf. eq. A9) also motivates the assumption—borne out by empirical evidence— that for typical weakly aspherical models, M remains diagonally (κκ) (κκ) 2 (κκ) 2 M(ω) ≈ V (ω ) + ı ω W − ω P − ω I (A11) 0 0 1 dominant. Elements lying far from the diagonal are, in general, where the superscript (κκ) is used to emphasize that each multi- relatively small—and in group- and self-coupling, they are therefore plet is being treated independently. This approximation relies on ignored. the assumption that |ω − ω | is small, and that the earth model is In eq. (A9), a general element of M is defined dominated by spherical structure, allowing various additional terms ∗ 2 3 involving P , W and ∂V/∂ω can be neglected. It is common to in- M (ω) = S (x) · V + ı ωW − ω ρ(x) S (x)d x. (10a) jk k 1 V ˜ troduce a ‘splitting matrix’, H = V(ω ) + ı ω W − ω P , allowing 0 0 1 eq. (A9) to be written in the form In self-coupling, this is modified by introducing the assertion that (κκ) 2 (κκ) H − ω I u˜ = q˜ . (A12) q = q ⎨ j k M = M = 0 unless l = l (10b) jk (q ,l ,n )(q ,l ,n ) j k j j j k k k Thus, the contribution of each multiplet to spectra—within the n = n . j k self-coupling approximation—is then found by diagonalizing the (2l + 1) × (2l + 1) matrix H.The(2l + 1) constituent singlets within In other words, self-coupling computes integrals using only the the multiplet each contribute at a slightly different eigenfrequency, constituent singlets of each multiplet, and assumes that interactions removing the degeneracy seen in spherically symmetric models. between different multiplets can be entirely ignored. Thus, the spectral contribution of each multiplet becomes a set of In group-coupling, the multiplets are assigned to groups. We (2l + 1) closely spaced—but distinct—peaks. This phenomenon introduce κ to denote a given multiplet, κ = {(q, l, m, n): m ∈ [−l, l]}. We then introduce groups, G, containing one or more is known as ‘splitting’. In practice, the resolution of spectral mea- (i) (i) multiplets: G ={κ ,κ ,...}. In general, all multiplets within a surements may not allow individual singlets to be distinguished, but 1 2 group will have similar eigenfrequency, although individual authors the overall shape of the peak corresponding to each multiplet de- may differ over the particular groupings adopted. Group-coupling pends upon the splitting, and thus upon the detailed properties of the then replaces eq. (A10b) by the assertion that Earth model. In the time domain, the contribution of the multiplet to seismograms can be written in the form (q , l , m , n ) ∈ G j j j j M = 0 unless (10c) jk (κκ) (κκ) (q , l , m , n ) ∈ G s (x, t) = r(x) · exp ıH t · v (A13) k k k k Downloaded from https://academic.oup.com/gji/article/213/1/58/4757069 by DeepDyve user on 13 July 2022 76 F. Akbarashrafi et al. where v and r are vectors that may be computed for the source and ing the reference frequency ω becomes more challenging, as receiver, respectively; this result is due to Woodhouse & Girnius the various multiplets within the group have different eigenfre- (1982). quencies. Typically, some sort of average of these eigenfrequen- The discussion in this section has focused on self-coupling in- cies is used (for a recent discussion on specific averages, see volving only individual isolated multiplets. Similar results can Deuss & Woodhouse 2001). We note that the requirement for be obtained for group-coupling–assumptions: in effect, group- |ω − ω | to remain small acts to restrict the size of the cou- coupling can be regarded as self-coupling of a ‘super-multiplet’, pling group: only multiplets with similar eigenfrequency may be allowing our analysis to translate straightforwardly. However, defin- used.

Journal

Geophysical Journal InternationalOxford University Press

Published: Apr 1, 2018

There are no references for this article.