Seismic event coda-correlation’s formation: implications for global seismology

The seismic-event-coda correlograms are characterized by many prominent features, which, if understood thoroughly, could supply valuable information on the internal structure of the Earth. To further reﬁne our knowledge and be able to utilize that information, all-embracing comprehension of coda-correlation’s formation apart from a conjecture, is a pre-requisite. Here, we conduct a comprehensive analysis that aims at a quantitative ‘dissection’ of the formation mechanism of coda correlation. Our analysis presents relevant implications for global seismology. We demonstrate that coda correlation is dominated by a few contributions, most of which arise from the late-coda time window, 3 hr after the earthquake origin time. Our identiﬁcation analysis conﬁrms that the contributions are cross-terms between body waves. That represents an observational proof of the conjecture that coda-correlation features are formed due to body waves arriving at a pair of receivers with the same slowness. We further quantify the relationship between body-wave cross-terms and event-receiver geometries and Earth structure, which has signiﬁcant practical implications. Our analysis demonstrates that body-wave cross-terms that contribute to the same coda-correlation feature sample the Earth along fundamentally different paths. They are signiﬁcantly different depending on event locations, although the resulting time variation is quite small if the late coda (e.g. 3–9 hr after event origin time) is used. That explains why the late coda is more effective than an earlier time window in producing relatively stable features, as empirically suggested by previous studies. Our study enables quantitative and practical understanding of coda-correlation features in terms of their formation progress, and this opens a way to distill valuable


I N T RO D U C T I O N A N D M O T I VAT I O N
Half a century has passed since pioneering investigations of the properties of coda waves-the late and long-lasting wave trains of earthquakes (Aki 1969;Aki & Chouet 1975). The discovery of the long-range correlation in seismic wavefield allows for significant advancement in near-surface imaging (Campillo & Paul 2003;Shapiro et al. 2005) and inferences on the Earth's deep interior (Poli et al. 2012(Poli et al. , 2015. In ambient-noise seismology, cross-correlation stacks resemble the Earth's response between two receivers as if the seismic wave arrivals are measured at one of the receivers from a virtual source at the other. This provides a scope for recovery of subsurface architecture through a dense volumetric coverage by means of many receiver pairs on Earth's surface and with less uncertainty than in conventional tomographic imaging due to the ability to circumvent ubiquitous source effects. Theoretically, cross-correlation of a diffuse wavefield converges to an elastic response of a medium between a pair of receivers (Green's function). Ambient-noise wavefield or seismic coda can be diffuse if the medium is heterogeneous (Lobkis & Weaver 2001) or there are many uncorrelated and equally distributed sources (Wapenaar 2004;Wapenaar et al. 2010). Those advances led to the reconstruction of surface waves from ambient-noise wavefield given the configurations of sources and receivers on the Earth's free surface (Shapiro & Campillo 2004) and numerous tomographic applications on regional and global scales (e.g. Yao et al. 2006;Moschetti et al. 2007;Bensen et al. 2008;Lin et al. 2008).
It is appealing to interpret those features in correlograms as reconstructed body waves (e.g. Huang et al. 2015;Wang et al. 2015): this interpretation is seemingly on a par with surface wave reconstruction from ambient noise. However, Boué et al. (2013) reported abnormal amplitudes of 'deep-Earth phases' in coda correlation. These features present much larger amplitudes in correlograms than in the direct seismic wavefield. Furthermore, regarding the particlepolarization pattern, Boué et al. (2014) reported an unusual ScSlike feature in vertical-to-vertical correlograms. Boué et al. (2014) described time inaccuracy of some features in cross-correlation functions when selecting different event-receiver geometries. More significantly, there are many stable and prominent features in correlograms that do not exist in direct seismic wavefield (Boué et al. 2013;Pha . m et al. 2018;Sager et al. 2018). Apart from that, the seismic-event coda on global scale is made of strong reverberations with energy confined in the great-circle plane (Sens-Schönfelder et al. 2015). This contradicts the necessary condition for accurate reconstruction of seismic waves in which the energy distribution should be equal in different directions for the wavefield that is cross-correlated. The unequal distribution of energy was also suggested by the dominance of high-quality-factor normal modes in earthquake coda (Maeda et al. 2006;Poli et al. 2017).
For the above reasons, it is essential to distinguish between the true nature of the features in coda correlation and the notion of 'reconstructed body waves' and to obtain an accurate understanding of coda-correlation's formation.  suggested that coda-correlation features are related to the late coda of earthquakes (3-9 hr since event origin time). Boué et al. (2014) sketched the involvement of coherent reverberations after large earthquakes in the emergence of coda-correlation features. Poli et al. (2017) demonstrated that the body-wave-like features in global correlograms are formed by the normal modes that correspond to seismic reverberations. Pedersen & Colombi (2018) further suggested the influences by cross-terms of seismic body waves on cross-correlation functions.
All features in global correlograms can be explained by means of a mathematical conjecture and simple laws of physics that employs ray theory: they arise due to the similarity between body waves of the same slowness for a given pair of receivers at the Earth's surface . In other words, both causal (where the correlation feature follow the timing of arrivals in the seismic wavefield) and non-causal features (where the correlation feature do not correspond to any arrivals in the seismic wavefield) in correlograms are produced by interferences between seismic waves, and the term correlation wavefield was coined. There can be many efficient interferences between seismic waves (Sager et al. 2018), particularly for steeply travelling waves (Kennett & Pha . m 2018a).
Motivated by the recent developments described above, here, we present a comprehensive and quantitative analysis that provides practical insights into coda-correlation's formation. In doing so, we quantitatively examine the essence of contributions to codacorrelation's formation based on a temporal variation analysis. The separation and identification analyses confirm that dominant contributions are cross-terms between body waves. That represents the observational proof to the mathematical conjecture of codacorrelation's formation Kennett & Pha . m 2018a) and provides deterministic and quantitative separation of reverberated body waves, whose involvement in codacorrelation's formation was suggested in previous studies (Boué et al. 2014;Poli et al. 2017).
On the one hand, separation of contributions explains the effectiveness of late earthquake-coda in forming correlation  and the development of correlation wavefield when different earthquake-coda segments were used (Poli et al. 2017;Kennett & Pha . m 2018b). On the other hand, we quantify the dependency of body-wave cross-terms on Earth structure and event-receiver geometries, and we determine the variations between different bodywave cross-terms. We demonstrate the mechanism of seemingly stable, but time-varying coda-correlation features, especially in terms of the time inaccuracy that was observed in previous studies (e.g. Boué et al. 2014). Our analysis enables quantitative interpretation and utilization of coda correlation to constrain Earth structure.

C O DA C O R R E L O G R A M S
We analyse earthquake-coda wavefield recorded by the networks and individual receivers in the Australasian region (Fig. 1a). This region has considerably low background noise due to the remote location of many receivers, away from the coastline and humanrelated noise (Kennett & Furumura 2008;Kennett & Furumura 2016). The source-receiver geometry is characterized by diverse azimuths and epicentral distances for globally distributed earthquakes (Figs 1a and b). We select earthquakes with magnitudes greater than M w = 6.8 for sufficient energy-radiation strength and depth larger than 20 km to suppress the excitation of surface waves (Poli et al. 2017;Pha . m et al. 2018). We experiment with different time windows, including 3-10 hr after the event origin time , and converged on the time windows of 3-9 hr after the event origin time as optimal to extract coda waveforms for the Australasian region. We perform temporal normalization to suppress surface waves or aftershocks within the coda, and spectral whitening to balance the energy in different frequency bands (Bensen et al. 2007;Pha . m et al. 2018). We filter cross-correlation functions with the 15-50 s bandpass filter for feature prominence. The cross-correlation functions are binned to inter-receiver distance bins of 0.5 • and stacked linearly with respect to the distance bins to produce the correlograms.
As shown in Fig. 2(a), the correlograms produced for the recordings in the Australasian region reveal similar patterns to those that utilize global-scale networks (e.g. Boué et al. 2014;Pha . m et al. 2018). The most prominent features are reminiscent of widelyobserved seismic phases: PcP, ScS and PKIKPPKIKP, and are termed PcP * , ScS * and I2 * following the convention for naming and abbreviation of . Other features are reminiscent of 'exotic seismic phases' such as PKPPcPPKP and PcPPKPPcPPKP that are termed KcK * , and cKcK * . Additionally, there are features such as cS-cP that have no correspondence in seismic wavefield. Sens-Schönfelder et al. (2015) and Poli et al. (2017) reported that earthquake coda presents unequal energy distribution in different directions and the energy is dominant along great-circle path from an event to a receiver. The lack of equipartitioning hints that coda correlation is dominated by the energy that propagates along the greater circle path. To investigate that, we first construct a global correlogram consisting of all events greater than M w = 6.8 without any further selection (Fig. 2a). We then experiment with selecting subsets of those events chosen with respect to the geometry between the event location and receivers. We first calculate the angle, ζ , between two great-circle paths defined by an event and two different receivers. If ζ is 0 • , the event and receivers are perfectly aligned and in the same great-circle plane. ζ can approach 90 • for some events given a pair of receivers, and we expect to see a strong dependence of coda-correlation features on it. Figs 2(b)-(f) shows the coda

Figure 2.
Global correlograms made up by stacking multiple contributions from inter-receiver distance bins of 1 • . The most prominent features are labelled following the notation defined in . (a) All events are used in the stacking. (b) Events are selected for the angle, ζ , smaller than 5 • . ζ is the angle between two great-circle paths that are from the same event to two receivers, respectively. (c)-(f) Events are selected for the angle, ζ , in ranges of correlograms for varying ranges of ζ . Coda-correlation features are clearly prominent when ζ is close to 0 • . The subset of events that satisfy ζ < 5 • (Fig. 2b) is sufficient to reproduce all coda-correlation features (e.g. ScS * , I2 * and cS-cP) seen in the global correlation stack from all events (Fig. 2a). The same is somewhat true for ζ in the range 5 • -15 • (Fig. 2c), although some features are significantly diminished in comparison with those in Figs 2(a) and (b). In contrast, using the events with larger azimuth difference (away from the great-circle plane) results in the correlograms with deteriorating correlation features (Figs 2d and f). In short, when using 3-D global stacks, only those in the great-circle plane constructively contribute to the formation of coda-correlation features.
The main result obtained above and demonstrated in Fig. 2 narrows the search range to determine the constructive contributions to coda-correlation's formation. That is the first stage in resolving coda-correlation's formation mechanism. Then, in the second stage, we further determine the exact constructive contributions. The twostage analysis might seem like reducing a 3-D problem to a 2-D problem. However, it is simply a way to determine the formation of coda-correlation features.
To further determine the exact constructive contributions to codacorrelation's formation, we therefore focus on earthquake-coda energy in the great-circle plane. We avoid non-contributing energy that is off the great-circle plane by setting an upper limit of 5 • for the angle ζ .

S E PA R AT I O N O F C O N T R I B U T I O N S T O C O DA C O R R E L AT I O N
We aim to separate and determine seismic waves in earthquake coda that contribute to coda-correlation features. In order to achieve this goal, we conduct a temporal variation analysis to determine the timings and energy of individual contributions to each coda-correlation feature and an identification analysis to determine seismic waves for each contribution.

Temporal variation analysis
In coda correlation, interactions between different parts of the earthquake-coda wavefield is manifold. Given the earthquake-coda time-interval, wave propagation through Earth's interior evolves significantly. Moving time windows can explicitly separate different parts of earthquake coda. We use a relatively short, moving time window of 100 s with 1 s overlap to calculate coda correlograms. 100 s is twice as long as 50 s, the lowest cut-off frequency of the bandpass filter, as required by the Nyquist-Shannon sampling theorem. That is significantly different from previous studies that use hourly time windows to investigate the optimal time window that produces prominent coda-correlation features (e.g.  or temporal development of correlation wavefield (e.g. Poli et al. 2017;Kennett & Pha . m 2018b).
Although earthquake coda from 3 hr after event origin time is used (e.g. Pha . m et al. 2018), we select earthquake coda between 1 and 9 hr after event's origin in order to answer why the emergence of coda-correlation features is favoured by the late earthquake coda starting 3 hr after event's origin.
We select the correlogram feature I2 * at inter-receiver distance of 5 • for its robustness and signal prominence. This feature is reminiscent of PKIKPPKIKP phase in the seismic wavefield (Fig. 3a). I2 * exists in the small inter-receiver distances and is well formed and pronounced for local-scale network settings. Additionally, I2 * has been used as a probe of the inner core structure (e.g. Huang et al. 2015;Wang et al. 2015). Fig. 4(a) shows the temporal variations of I2 * using the short time window. To quantify the temporal variations, we calculate the coherence of I2 * as a function of time. Following Rawlinson & Kennett (2004), the coherence can be measured by taking the inverse of the sum of the squared residual between normalized I2 * using a 100 s time window and I2 * using 3-9 hr: where the τ denotes the time from the event's origin. The coherence denotes the similarity between I2 * feature using a moving 100 s length time window (I2 * 100 s ) and the prominent I2 * in global coda correlograms that are produced by using 3-9 hr after event's origin (I2 * 3-9 hr ). A relatively large coherence indicates that our chosen time window centred on τ captures a constructive contribution to I2 * , whereas low coherence means that corresponding contributions are minimal. Fig. 4(b, orange line) presents the coherence variation with the time τ . The spikes indicate discrete contributions to the formation of I2 * . There are fewer coherence spikes in the time interval of 1-3 hr than in the time interval of 3-9 hr.
Apart from the measurement for coherence, energy is equally important in investigating contributions. The energy strength of a single contribution can be measured following Sager et al. (2018): As shown in Fig. 4(b, blue line), the energy is relatively strong and oscillatory in 1-3 hr and gradually decreases with the time τ and becomes stable in later time windows. Only the contributions with large coherence and relatively strong energy can be constructive contributions to correlation features, while those with small coherence but strong energy will 'contaminate' correlation features. Finally, those with weak energy and either small or large coherence are negligible.
We investigate temporal variations of various contributions by combining the measure of coherence and energy. For the time interval of 1-3 hr, most of the contributions present small coherence and oscillatory energy, and there are a few spikes in the coherence plot however corresponding energy is weak. The former would contribute destructively and overwhelm the latter in correlationfeatures' formation. The observed pattern makes the coda wavefield within 1-3 hr ineffective in contributing toward I2 * . For the time interval of 3-9 hr, more contributions emerge with relatively high coherence and strong energy. They dominate the formation of I2 * and suppress phase deformations by incoherent contributions. This time-dependent variation of different contributions explains the effectiveness of the late coda in forming prominent features in correlograms, as suggested in previous studies Poli et al. 2017;Kennett & Pha . m 2018b).

Identification of coda-correlation constituents
We identify constructive contributions (constituents) in codacorrelation features, by matching the time τ of coherence spikes (Fig. 4b) with the arrival times of transient seismic waves, which can be predicted by the ray theory (e.g. Chapman 1983). We attribute a constituent represented by a coherence spike to the cross-term between two seismic waves. Specifically, a constituent at time τ is the cross-term between two seismic waves that arrive at τ and τ + t after event origin time, respectively, where t is the correlation time of a correlogram feature. Arrival-time predictions depend on event locations, although the dependency is moderate for I2 * that presents near-zero slowness features (Fig. 2). In the next section, we further analyse the influence of event locations on the correlation time of correlogram features.
As shown in Fig. 4(c), almost all of constructive contributions are identified as cross-terms between seismic waves that reverberate between the Earth's free surface and the core-mantle boundary (CMB). Contributions are insignificant for seismic waves related to the Moho and upper-mantle discontinuities. As for smaller 3-D heterogeneity, associated scattering effects would be less likely to contribute to the formation of coda correlation. Therefore, coda correlation consists of limited types of constituents that are cross-terms between the deep-Earth phases. This predominance of deep-Earth phases agrees with the lack of equipartitioning in the late coda (Sens-Schönfelder et al. 2015;Poli et al. 2017). These discrete constituents of coda correlation could be thought of as an observational proof of the conjecture that the correlation wavefield is formed by interferences between body waves .  Further to analysing I2 * , we conduct a similar temporal analysis of another coda-correlation feature: cS-cP (Fig. 3b) and its constituents. We select cS-cP at inter-receiver distance of 2 • . This relatively small inter-receiver distance allows a sufficient number of global events when applying the same great-circle plane requirement (ζ <5 • ). cS-cP displays much faster deterioration than I2 * if the same great-circle plane condition is not satisfied (Fig. 2). As shown in Fig. 5, there are discrete contributions to cS-cP, identified as cross-terms between reverberated body waves from the Earth's free surface and CMB, and most of the contributions are in the time interval 3-9 hr after event's origin.

C H A R A C T E R I S T I C O F C O DA C O R R E L AT I O N F RO M I T S C O N S T I T U E N T S
We conduct a quantitative analysis in order to enhance practical understanding of coda correlation. Once coda-correlation constituents are fully isolated and identified, we can quantify their relationship with Earth structure and dependency on event-receiver geometry.

Attributes of coda-correlation constituents
A notable feature is the temporal fluctuation of separate constituents of I2 * or cS-cP (Figs 4c or 5c). That variation is on the order of 10 s, which is equivalent to a P-wave isotropic velocity anomaly of ∼2 per cent in the inner core that is sampled by I2 * . The variation among waveforms is apparent both in amplitude and phase, which means that the contributing constituents have different sensitivities to Earth structure. This can be expected given that each pair of seismic waves, for example those listed in Figs 4(c) and 5(c), sample the Earth's interior along different paths. Varied focal mechanism and source-time functions are other reasons for the waveform differences among separated constituents, although they generate less significant time variations than Earth structure.
The correlation time of coda-correlation constituents and their variations contain information directly dependent on Earth structure. Each correlation time corresponds to the time difference between two seismic waves that can be measured by the cross-term between them. That can be quantified as: where the correlation time of a coda-correlation constituent t cc constituent can be predicted from the time difference between two seismic wave arrival times t wave 1 (x s , x r 1 ; m) and t wave 2 (x s , x r 2 ; m) that depend on Earth structure m and locations of an event x s and two receivers x r 1 and x r 2 . We assess that quantification through a synthetic experiment as follows. We use the axisymmetric spectral element method, AxiSEM (Nissen-Meyer et al. 2014), to simulate synthetic coda wavefield for the actual event-receiver geometry settings. In our simulations, we employ a spherically symmetric Earth model ak135 (Kennett et al. 1995). We use explosive mechanism and Gaussian source time function with a dominant period of 10 s following Pha . m et al. (2018). We consider two fixed receivers with a 5 • separation and events distributed globally along the great-circle path defined by the receivers (Fig. 6a). We choose the cross-term between I6 and I4 phases for our analysis, and we calculate cross-terms for events at different longitudes. We do not consider varied focal mechanisms or source time functions in the simulation to keep the synthetic experiment simple.
As shown in Fig. 6(c), the correlation time of cross-term I6-I4 equals the arrival-time difference between I6 and I4 that are predicted using ray theory. Significantly, the cross-term I6-I4 exhibits variation in correlation time when events are in different locations (hyperbolic in shape as a function of event longitude). The time variation curve has a saddle-point minimum around which the curvature is nearly flat (Fig. 6c). At the minimum, I6 and I4 have the same slowness and a portion of the ray paths in common (Figs 6b  and d). Only at the point, the differential ray path equals to PKIKPP-KIKP between two receivers. This equivalence does not exist for other event locations, and hence there are time variations, waveform distortion, and even the absence of I2 * in extreme cases (Fig. 6c). Similar time variations hold for other constituents, for example I8-I6, I10-I8, etc. (Fig. S1, Supporting Information). Eq.
(3) denotes a quantitative relationship between codacorrelation and Earth structure. Given an Earth structure model, the correlation time of a coda-correlation constituent can be predicted. Through decreasing the difference between the predictions and observations, an inverse problem can be built towards inverting Earth structure.

Formation of a coda-correlation feature by many constituents
The time variations introduced due to event locations have a significant bearing on the formation of the coda-correlation features. When the time variations around the same slowness point, as shown in Fig. 6(c), are nearly zero, cross-terms contribute constructively to cross-correlation stacks to produce a stable correlation feature. In contrast, large time variations result in non-constructive contributions whose cross-terms dilute the stacks. We further group events according to their distance to the same slowness point and stack I6-I4 cross-terms to produce I2 * feature. As shown in Fig. 7, the stacked I2 * feature has a relatively stable timing if we limit events to only those within 10 • around the same slowness point. There is a severe time shift if all events on global scale are used: this results in a time shift of I2 * up to ∼4 s. Such a time shift is equivalent to ∼1 per cent P-wave isotropic-velocity perturbation of the entire inner core sampled by I2 * . Furthermore, the time variations due to event locations are significantly different among different I2 * constituents, which is illustrated in Fig. 8. The cross-terms between seismic waves in the late earthquake-coda (3 hr after the event origin time) present much smaller time variation ranges than those formed by the seismic waves in early time windows (1-3 hr from the event origin time). For example, I6-I4 presents time variation ranging up to ∼61 s due to event locations, whereas I24-I22 has much smaller time variation of ∼6 s. That pronounced difference explains why the published late-coda correlograms favour the emergence of codacorrelation features (e.g. Poli et al. 2017). Namely, only for the late coda, body-wave cross-terms exhibit less location-dependent time variations, and hence the cross-terms can effectively produce stable correlation features in the process of stacking.
In terms of comparing the 2-D and 3-D effects, the crossterms for 3-D Earth include both those in the great-circle plane and those outside the great-circle plane. As demonstrated in Section 2, the former mainly contribute to the formation of codacorrelation features, while the latter do not contribute significantly. Therefore, stacking of all cross-terms for 3-D Earth becomes well represented by a synthetic experiment reduced to the great-circle plane configurations. The only difference is that the global-scale stacks include more non-contributing cross-terms of seismic waves that propagate outside the great-circle plane in three dimensions.
The relationship between cross-terms of seismic waves and Earth structure always varies significantly with respect to event locations, although cross-terms presents small time variations when late earthquake coda is used. As shown in Fig. 9, from different events, ray paths are completely different for cross-terms between seismic waves that are either in earlycoda such as I6-I4 or in the late coda such as I22-I20 and I28-I26, although the latter present much smaller time variations.
Therefore, to produce stable coda-correlation features, we can either use the late coda, as suggested in previous studies (e.g. Poli et al. 2017) or carefully select events and crossterms in early time windows. If using the late coda without event selection, it is important to recognize the variation of cross-term ray paths, although the associated correlation time can be nearly invariant. A good example is one of the I2 * constituents, I28-I26, whose time variation is ∼5 s in maximum (Fig. 8b). Fig. 9 shows the ray paths for some of the I2 * constituents for the events located in the stationary or the point of the same slowness for the constitutive phases (left-hand column) and randomly selected event locations (right-hand columns). For example, I6-I4 is shown in Fig. 9(a), and I28-I26 is shown in Fig. 9(d). Even when the event is located in the corresponding stationary point, the sampled volume of the Earth increases considerably for I28-I26 in comparison to I6-I4, and we lose the lateral resolution of the deep earth.

D I S C U S S I O N S A N D C O N C L U S I O N S
We have shown that coda-correlation consists of discrete and timedependent contributions. Notably, most of these contributions reside in the late coda, 3 hr after the origin time and later, while the early time-interval presents many incoherent constituents with varying energy. There are quite a few coherent contributions in 1-3 hr, though their energy is relatively weaker than the energy of incoherent ones. Specifically designed short-time windows based on temporal analysis could remove incoherent contributions and their 'contamination' effects. This makes it possible to promote codacorrelation to an early time interval.
Our quantitative separation and identification procedures confirm that coda correlation is dominated by cross-terms between seismic phases that are sensitive to deep-Earth structure and propagate in Figure 6. An 'anatomy' of the cross-term I6-I4, one of the I2 * constituents. (a) Synthetic distribution of sources (stars) and two fixed receivers (triangles) that are used in analysing characteristics of I6-I4. The ray path of I2 (black line) from one receiver to the other. (b) Ray paths of I6 (blue lines) and I4 (orange lines) from the event at the stationary (the same slowness) point to the two fixed receivers. (b) Variations of I6-I4 cross-term as a function of event locations. The correlogram is created using the distribution of events in (a). Red dash line denotes theoretical arrival-time difference between I6 and I4 based on the ak135 Earth model (Kennett et al. 1995). (d) Theoretical slowness curves for I4 (solid line) and I6 (dash line). The black dots denote the event location where I4 and I6 have the same slowness and the I6-I4 differential time curve has the curvature of 0. the great-circle plane. Those seismic phases stem from reverberations between the Earth's free surface and the CMB, whereas the interactions with minor discontinuities and lateral heterogeneity are insignificant. This represents an observational proof to the mathematical conjecture suggested recently Kennett & Pha . m 2018a). Apart from that, it indicates that the late-coda wavefield is not diffuse, a concept equivalent to the unequal energy distribution (Maeda et al. 2006;Sens-Schönfelder et al. 2015;Poli et al. 2017). Instead, the late coda is predominated by deep-Earth phases and hence nearly zero-slowness Figure 8. (a) Differential correlation time curves for different I2 * constituents, with respect to event locations. The colour dots denote the events located in the stationary (the same slowness) points for corresponding cross-terms. (b) Correlation-time variation maxima for different cross-terms due to the selection of events located away from the stationary points. Time variation maxima significantly decrease for cross-terms of higher reverberated body waves.
There are significant differences between the cross-correlation of late earthquake-coda wavefield and the cross-correlation of a diffuse wavefield (e.g. ambient noise wavefield), in which crosscorrelation features are equivalent to reconstructed seismic waves. In a diffuse wavefield, energy is equally distributed in different directions. 'Non-causal' cross-terms in individual cross-correlation functions destructively interfere and cancel each other in stacking, and hence cross-correlation stacks converge to reconstructed seismic waves (Snieder et al. 2008;Snieder & Fleury 2010). However, cross-correlation functions are dominated by cross-terms of deep-Earth phases. Stacking them not only does not cancel them out but results in 'non-causal' features, such as cS-cP.
In the coda-correlation wavefield, seismic-wave cross-terms present different time variation with respect to the selection of events. Cross-terms of seismic waves in late earthquake-coda present nearly invariant correlation time and hence their stacking results in constructive interference to produce stable correlation features. This explains the increasing number of contributions after 3 hr from the event origin time (Figs 4c or 5c). The correlation time of those cross-terms presents less sensitivity to the event locations, and hence they are quite stable when global events are used without selection. In contrast, cross-terms from early time windows present relatively large time variations.
The insensitivity to event locations in the late coda has practical implications. A stable emergence of coda-correlation features does not necessarily mean a stable reconstruction of body waves. The stable emergence results from the nearly invariant timing of crossterms. However, the ray paths vary significantly with event locations (Fig. 9). Therefore, cautions must be exercised in all interpretations or utilizations, and it is a necessary requirement to investigate the ray paths related to event locations.
Coda-correlation tomography is feasible based on our comprehensive analysis. The sensitivity kernel for tomography can be constructed by relating each of the identified constituents to Earth structure based on eq. (3). With informed choice of events, the differential sensitivity between two seismic waves can be obtained. The tomography kernel consists of many differential sensitivities, each of which corresponds to an identified seismic-wave cross-term. This approach is similar to building differential kernels (e.g. Calvet & Chevrot 2005;Ruigrok et al. 2012). The sensitivity and accuracy can be evaluated with either ray-approximation-based routines (e.g. Cervený 2005) or finite-frequency methods (e.g. Marquering et al. 1999;Dahlen et al. 2000a,b). The kernel construction framework is different from the conventional one in which inter-receiver crosscorrelation stacks are treated as Green's function. In the Green's function approach, all of the available events and the full coda time window are used. In contrast, the new approach precisely relates each constituent of correlogram features to Earth structure. Also, our analyses show that there are multiple constituents for each feature in correlograms. Therefore, even a single feature in correlograms may provide a well-posed inference on Earth structure.
Finally, the separation and identification of constituents are labour-intensive steps. Some constituents cannot be separated. For example, PcSPKIKP-PcS and PcSPKIKP-ScP possess nearly identical arrival time. They cannot be separated from each other though Figure 9. Ray paths of different I2 * constituents. (a) Ray paths for cross-term I6-I4 for an event in the stationary (the same slowness) point (left-hand column) and for two events away from the stationary points (middle and right-hand columns). Time difference in seconds for each event is shown above each cross-section. (b)-(d) similar to (a) but for cross-terms of I14-I12, I22-I20 and I28-I26. they have fundamentally different sensitivities to Earth structure. It is necessary to avoid this kind of indistinguishable cross-terms and select other, unambiguous cross-terms for accurate tomographic inference. We can use short-time windows around the arrival times of the phases we select, for instance PKIKPPKIKP and PKIKP, and calculate the cross-term between them (Wang & Tkalčić 2020). Also, the computation and visualization requirements are significant due to the short moving time windows, and modern high-performance computation facilities are necessary. Despite these obstacles, the study presented here is a foundation for a quantitative utilization of coda correlation, and we look forward to further improvements for standard and target-oriented approaches in different seismological problems.

A C K N O W L E D G E M E N T S
The computation and visualization were performed on the ANU Terrawulf cluster, a computation facility developed with support from the AuScope initiative. Waveform data used in this Downloaded from https://academic.oup.com/gji/article-abstract/222/2/1283/5843724 by The Australian National University user on 31 July 2020