Applying the concept of cluster entropy, we demonstrate that the primary P-to-s (Ps) conversions originating from the well-recognized seismological boundaries (Moho, 410 km, and 660 km depth interfaces) in Earth's interior and multiple reflections (reverberations) between the surface and the Moho—the prominent crust-mantle boundary of Earth—can be unambiguously identified. Their attendant information dimension can act as a discriminant. Documentation of shallow mantle stratification at depths ∼150–300 km (L/X layers) using Ps signals is scarce and is sometimes attributed to their weak registration in the seismograms together with interference by reverberations from shallow structure. The slopes of best fits to mean cluster information dimension (as function of epicentral distance) of Ps signals from target boundaries such as the Moho, 410 km, and L/X discontinuities are consistently gentle with a slight negative character. In contrast, those related to multiples show a steep positive behavior. This opposite nature can potentially discriminate between reverberations and direct converted waves, even when they tend to arrive in overlapping time windows. Our approach thus enables unequivocal identification of shallow mantle layering from Ps data recorded at diverse tectonic provenances of wide antiquity.
A stratified Earth model receives support from the presence of prominent seismic velocity discontinuities at average depths of (1) 33 km (termed as the Moho), (2) 410 km (410 km), and (3) 660 km (660 km) detected globally, in addition to the discovery of a shallow mantle boundary at ∼200 km, first recognized and reported by Lehmann (1959), which is restricted to certain regions. Since then, several researchers have documented the presence of shallow mantle layers in diverse tectonic settings using a variety of seismic waves at varied depths (Revenaugh and Jordan, 1991; Li et al., 2002; Gung et al., 2003; Deuss and Woodhouse, 2004). Based on their depth disposition, these boundaries are termed the Lehmann (L, 180–250 km) and X- discontinuities (250–330 km). With the advent of P-to-s (Ps) receiver function technique (Vinnik, 1977), owing to its superior vertical resolution imaging power, unambiguous detection of a stratified mantle by way of well-documented presence of the 410 km and 660 km discontinuities almost everywhere became possible. For the same reason, delineation of shallow mantle layering seemed promising. In spite of their better sensitivity to vertical layering, observations related to the possible presence of the shallow depth boundaries ∼150–300 km in the Ps receiver function data generally remained scarce for several reasons. There is a perception that weaker registrations of Ps conversions from these shallow mantle depth layers, compared to those related to the well-documented 410 km and 660 km discontinuities, often render the former relatively ambiguous (e.g., see Grunewald et al., 2001). On the contrary, the amplitudes of Ps conversions related to the Lehmann discontinuity (∼250 km; PLs) are the highest beneath the northwestern half of the Colorado Plateau–Rio Grande Rift–Great Plains Seismic Transect (LA RISTRA) (Wilson et al., 2005) and eastern segment of the Missouri to Massachusetts Broadband Seismometer Experiment (MOMA) profile (Li et al., 2002). Notwithstanding these data, importantly, primary Ps converted waves from the 150–300 km depth interfaces and reverberations (multiples) due to shallow structure can arrive in overlapping time windows. Several studies can be cited in this context: (1) Li et al. (2002) reported complex scattering originating from shallow structure in the western segment of the MOMA profile that interfered with PLs-like arrivals; (2) the nature of arrivals corresponding to 250–300 km depth beneath the Tanzania craton remains inconclusive (Owens et al., 2000); and (3) the problem of unambiguous detection of PLs arising due to interference by multiples from Moho and sediments was reported by Landes et al. (2006). Traditional methods based on moveout trends (Yuan et al., 1997) remained unhelpful in such situations. Such a dilemma was also expressed recently while dealing with a negative phase in study of the lithosphere-asthenosphere boundary in North America from receiver functions (Abt et al., 2010). Details of the factors that render detection of depth boundaries at ∼150–300 km in the Ps data relatively ambiguous and contentious can be found in Sheehan et al. (2000). Hence, given the proximity of crustal reverberations to possible direct converted phases, such as the PLs in our study, this signal cannot be unambiguously attributed to a direct phase.
The opposing trends of the moveout curves for Ps and the multiples based on slowness criteria (Yuan et al., 1997) together with some innovative data segregation (Ramesh et al., 2005), or with S-to-p (Sp) data as a guide (Farra and Vinnik, 2000; Ramesh et al., 2010), would reduce the ambiguity associated with identification of Ps conversions from shallow mantle boundaries. Though the opposite moveout of primary conversions and multiples can potentially discriminate between them, its success becomes partial when the receiver functions are restricted to a narrow range of slowness values, and also due to possibly small moveout values of PLs-like phase arrival times with ray parameter (Sheehan et al., 2000). Therefore, techniques other than those based on slowness criteria alone deserve attention. One fundamental reason for the limited success of the moveout-based discrimination of seismic phases could be that we are still operating in the measurement space itself. Instead, it may be interesting to view data in the generic space where the inherent differences related to the origin of various seismological phases could be preserved better. One such approach involves application of principles of information theory to our observations.
Use of information theory in general, particularly employing methods based on measures of entropy to solve geophysical problems, is still in its infancy and has gained momentum in recent years. For example, various measures of entropy were successfully applied to investigate dynamical complexity of the Earth's system (Earth's magnetosphere, to be precise) recently by Balasis et al. (2009). Their results provide conclusive evidence for a significant decrease in complexity and accession of persistency in disturbance storm time (Dst) series as a magnetic storm approaches. This finding, according to them, can be used as a diagnostic tool for detection of magnetospheric injury (global instability). To our knowledge, use of measures of entropy in seismology is yet to gain prominence. The present study to discern seismic phases of varied origin by applying a measure of cluster entropy is perhaps a pioneering attempt to understand Earth's stratification.
Documentation of shallow mantle layering assumes significance both in the context of mantle stratification and candidate mineral transformations at these depth regimes (Deuss and Woodhouse, 2004; Ganguly and Frost, 2006). In this work, based on the concept of cluster entropy (Carbone et al., 2004; Carbone, 2007; Carbone and Stanley, 2007; Balasis et al., 2009; see also Shannon, 1948), we determine the cluster information dimension of primary Ps conversions from various depth discontinuities (e.g., Moho:Pms; Lehmann:PLs; 410 km:P410s) and the crustal reverberation phases (e.g., Pps and Pss; see Fig. 1). We further demonstrate that the cluster information dimension can indeed discriminate between direct converted waves and multiply reflected phases registered in Ps receiver functions.
MOTIVATION, METHODS, AND DATA
The reflection and transmission coefficients together with the energy flux associated with converted phases such as Pms, PLs, and crustal multiples (viz. Pps and Pss) as a function of incidence angle were computed and are presented in Figure 1 based on approaches suggested by Aki and Richards (1980) and verified using Young and Braile (1976). The most significant observation is that the conversion coefficients (displacement amplitudes) and the associated energy flux values of Ps conversions (e.g., Pms and PLs phases), and the reverberating phases (Pps and Pss) follow distinctly different distributions (Fig. 1). Also, the varied converted and multiple phases show internal consistency in consonance with their nature. Based on the typical ray geometries of converted and multiple phases and their attendant distinct energy flux distributions shown above (see Fig. 1), it is reasonable to anticipate that their respective average amount of information (entropy) could be different. This possible difference in entropy between the direct conversions and reverberation phases would result in divergent information dimensions of these seismic arrivals. Hence, some measure of attendant information dimension can perhaps potentially discriminate between primary converted waves and multiply reflected phases that might arrive even in overlapping time windows.
Method to Estimate Cluster Information Dimension (DCluster)
In order to establish tangible criteria for the detection of shallow mantle layers and to have better insight into the nature of varied seismic phases registered in Ps receiver function data, we calculated cluster entropy of a time series adapting the approaches mentioned in Carbone et al. (2004), Carbone (2007), Carbone and Stanley (2007), and Balasis et al. (2009).
In information theory, any data string “X” is termed as a source, and the data points constitute the source alphabet. Prior to sampling of this source, there is a certain amount of uncertainty associated with the outcome, which turns into some information on sampling. Thus, concepts of uncertainty and information are related. Cluster entropy is defined as “S” = − Σ Pi × log Pi, where Pi is the probability of occurrence of ith event. The quantity log Pi is the amount of information one gains when an event with probability of occurrence Pi happens. Hence, the average amount of information or the expected value of information for one trial of experiment is the quantity “S.” When all events are equally likely, then Pi = 1/N, where N is the number of options. Hence, the average information will be log Pi = log N. To calculate the cluster entropy “S” of a given data string, the samples are first grouped into histograms of stipulated bin width, and the area under the histogram is normalized to yield the probability distribution. The scaling exponent (factor) that describes the scaling of the computed entropy with bin width is termed the cluster information dimension. Therefore, cluster information dimension (DCluster) is calculated using DCluster = Σ Pi(r) × log Pi(r)/Σ Pi(r) × log ri, where ri denotes the bin width of the histogram.
The important steps used to calculate the cluster information dimension as applicable to our study in a practical recipe form are: (1) imparting stationarity to the time series by applying a moving average to the data (Fig. 2), and (2) obtaining the differential amplitudes between the time series and its moving average (Fig. 2), to result in a stochastic sequence, which is followed by (3) construction of frequency distributions of this amplitude sequence in order to compute its cluster entropy using the corresponding probability distribution functions associated with each class. (4) Later, from the ensemble of entropy of all the classes, the cluster information dimension of the seismic phase under consideration is estimated according to the previous formula for DCluster, and (5) finally, the entire process (1–4) is repeated on each receiver function to enable us to compute DCluster of various seismic phases (e.g., Pms, PLs, P410s, Pps, Pss, etc.) for the total epicentral distance range of the data.
Based on this foundation, we first demonstrate the efficacy of our technique on synthetic receiver function data generated using a model with the “Lehmann” boundary incorporated (see Fig. 3A). The synthetic data were obtained adopting standard methods (e.g., Haskell, 1962; Langston, 1977, 1979). The DCluster values as a function of epicentral distance for various seismic phases (Pms, PLs, Pps, and Pss) along with their best fits are shown in Figure 3B. It is interesting to note that the best-fit slopes of primary converted phases (Pms and PLs) show a negative trend, in contrast to the positive slopes exhibited by the multiples (Pps and Pss phases). This opposing slope-character distilled from the synthetic seismogram data could form the fundamental basis to discriminate between depth converted phases and multiply reflected seismic phases.
Application to the Observed Data
With this premise, a few hundreds of broad band data from six seismic stations distributed mainly in India and North America, and published by us earlier (Ramesh et al., 2002, 2010), were reanalyzed. Two of these stations (HYB—Hyderabad; CUD—Cudappah) are located in southern India, two (FFC—Flin Flon; FRB—Frobisher Bay) are situated in Canada, while the remaining (HRV—Harvard; PAS—Pasadena) are located in the United States. We followed the methods described by Vinnik (1977) to construct and treat the Ps receiver functions. The number of data, their quality, and processing details together with the published main results for these stations from India and North America can be found respectively in Ramesh et al. (2010) and Ramesh et al. (2002). Following Vinnik (1977), to obtain a Ps receiver function, we first rotated the original ZNE recordings to ZRT using back-azimuth information. Subsequently, the Z and R components were further rotated into local ray coordinate system L (P energy) and Q (SV energy) components involving the angle of incidence for better isolation of the P and SV energy from the incident wave field. This was followed by deconvolution of the L component from the Q component to obtain P-to-s (Ps) receiver functions.
The Ps distance moveout-corrected (Yuan et al., 1997) receiver function stack sections for six stations are presented in Figure 4. Apart from the well-recognized Pms, Pps, and Pss arrivals, positive polarity energy corresponding to possible PLs(?) arrivals can be observed in the data at ∼20 s consistently across a large distance range (slowness) at three seismic stations, namely, PAS, CUD, and HYB (Figs. 4B–4D). This additional feature (PLs), though observed at these stations in the aforementioned studies, was not discussed earlier in view of its arrival proximity to the multiple phases (Pps and Pss), and its often inclined nature akin to the multiples (see Figs. 4B–4D) in the Ps moveout-corrected sections. These factors render identification of possible PLs phase as a primary conversion ambiguous. For these reasons, we then preferred to err on the side of caution, especially in the absence of robust primary and multiple phase discriminants, in addition to the lack of additional constraints or information in the form of Sp receiver functions, which are now available at some locations to unambiguously identify these signals (see, for example, Ramesh et al., 2010).
Therefore, our motivation to apply the concept of cluster entropy (Carbone, 2007; Balasis et al., 2009; see also Shannon, 1948) to study the nature of energetic arrivals close to the identified Pps and Pss arrivals, preferably within a 20–30 s delay time window, present in the Ps moveout-corrected stack sections (e.g., Fig. 4) of our previously cited publications seems reasonable. Thus, we computed the cluster entropy and the corresponding cluster information dimension for the well-recognized direct converted phases Pms and P410s, and multiples Pps and Pss, in addition to the ambiguous PLs arrivals.
In computation of the cluster entropy, ∼15% of data did not yield a cluster using the attendant moving averages of the corresponding seismic phases. Temporal variations in the ambient noise beneath a seismic station could perhaps be one contributing factor to this. This minor limitation can however be overcome with availability of moderate amounts of data.
RESULTS AND DISCUSSION
The computed cluster information dimension of the seismic arrivals designated as Pms, PLs, and P410s, in general, share a narrow restricted range, as opposed to the wide scatter exhibited by the multiples Pps and Pss (Fig. 5A). This indicates that the range of uncertainty associated with the ambiguous PLs phase reflects more commonality with the recognized converted waves (Pms and P410s) compared to the well-established reverberated arrivals (Pps and Pss). This could perhaps serve as one diagnostic to discriminate between direct conversions and the multiples.
We constructed the probability distribution functions (pdf) of the cluster information dimension corresponding to the seismic arrivals Pms, PLs, P410s, Pps, and Pss at various seismic stations. Samples of such pdfs at representative stations presented in Figure 5B illustrate the distinct distribution character of the converted waves and the multiples. This observation is valid in most of the cases. Such contrasting behavior of the candidate phases is already documented in their respective energy flux and conversion coefficient distributions (Fig. 1). This observation reinforces the fact that reverberation phases seem to follow a slightly skewed distribution, while primary conversions from target boundaries within Earth are invariably associated with more symmetric distributions (Fig. 5B). This could serve as another diagnostic to discriminate shallow mantle primary conversions (like the PLs phase) from the multiple arrivals.
To further distill our diagnostics and observations in order to arrive at a more robust measure of the discriminant, the cluster information dimension (DCluster) of Pms, PLs, P410s, Pps, and Pss seismic phases, as a function of epicentral distance, was studied in two ways. The first deals with variation of DCluster corresponding to each receiver function with epicentral distance (Fig. 6), which replicates the behavior already documented by the synthetic data presented in Figure 3B.
In the other, the Ps receiver functions obtained over a wide distance range (30°–90°; 1° ≈ 111 km) at each station were grouped into uniform-sized distance bins. The choice of the distance bin width at each station varied in order to maintain a comparable number of receiver functions per distance bin. The mean value of the cluster information dimension associated with the stated direct and multiple phases corresponding to each distance range (bin width) were then computed. A weighted least square fit corresponding to each seismic phase mean information dimension at a few stations sampling diverse tectonic regimes of varied ages is shown in Figure 7. Also, individual mean information dimension values for phases Pms and Pps are shown. Importantly, all the primary converted phases (Pms, P410s) show nearly flat to marginally negative best-fit slopes, while those arising from the multiples (Pps and Pss) are characterized contrastingly by relatively steep positive slopes. Significantly, the weighted best fit for the ambiguous PLs phase shares similar attitude and attributes as the confirmed direct converted waves.
From this, we find that the observed opposing slope-character is indeed an important diagnostic to discriminate between primary converted waves and multiply reflected phases. Our results from individual (Fig. 6) and binned receiver function data (Fig. 7) imply that the nature of observed opposing slopes for direct and reverberated phases would be preserved even when some epicentral distances remain unrepresented or under-represented. We would therefore like to add that a large amount of data is not a stringent requirement in our approach to extract the diagnostics, because even moderate amount of data (Fig. 3B) seems to yield the desired results without influencing our final conclusions significantly.
Another important issue that arises is related to the choice of bin width and its influence on the observed slopes. This is addressed by focusing on the PLs observations recorded at stations PAS, HYB, and CUD, where the distance bin widths at each station were varied in order to maintain a comparable number of receiver functions per distance bin. Figure 8 demonstrates that the nature of the best-fit slope to the PLs phase at these stations is still preserved as in individual station plots (Fig. 7), directly suggesting that the choice of bin width has no significant role in determination of DCluster and on our final results.
Plots similar to these (Fig. 7), combining cluster information dimension results from all the six stations corresponding to each seismic phase, are presented in Figure 9. Considering their location on distant continents of the globe with different age provinces of tectonic diversity, the consistent adherence to gentle-negative slopes by the direct phases (Pms, PLs, and P410s) and to steep positive slopes by the multiples (Pps and Pss) in the mean cluster information dimension versus average epicentral distance bin plots at all locations (Fig. 9) is certainly remarkable. This therefore constitutes a robust measure for discriminating between the competing seismological phases of interest.
CONCLUSIONS AND IMPLICATIONS
This is the first attempt to discriminate between various seismological phases, such as the primary converted waves arising from several depth boundaries in Earth's interior, and the multiply reflected wave types using the concept of cluster entropy and related information dimension. After demonstrating that the distributions of amplitude displacements (conversion coefficients) and corresponding energy fluxes of the direct converted waves and reverberation phases are different, we exploit these dissimilar characteristics to compute their cluster information dimension. Three distinctive diagnostics have emerged to discriminate between primary conversions and multiples upon application of Shannon's theory to the data. These are (1) dichotomy in their DCluster range of uncertainties, (2) distinct DCluster distribution character of corresponding pdfs, which are more symmetric for conversions and skewed for reverberations (Fig. 5B), and, finally, (3) consistently diverse mean DCluster best-fit slopes (Figs. 6, 7, and 9) of the primary converted waves and the multiply reflected wave types registered globally. We therefore successfully address the relatively long-standing issue of existence and unambiguous detection of shallow mantle layering (∼150–300 km depth) from Ps receiver function data by applying concepts based on information theory. The reliability of this new approach is demonstrated by the cluster entropy-delineated PLs phase as confirmed by the Sp receiver function results of Ramesh et al. (2010) at stations HYB and CUD. The confirmed presence of velocity layering beneath the Indian stations, demonstrated by us in this study, at depths akin to the Lehmann boundary, along with the results presented in Ramesh et al. (2010), has at least one important implication. When viewed together with the latest reports on the nature of lithospheric layering in the North American craton, and its attendant implications (Yuan and Romanowicz, 2010), our findings clearly suggest that the depth to the lithosphere-asthenosphere boundary beneath the cratonic regions of SE India definitely exceeds 200 km.
Our new approach may provide a platform to test the validity of the existing phase transition hypotheses on the origin of shallow mantle layers by way of comparison with cluster entropy associated with candidate mineral phase reactions from mineral physics experiments. This has wide implications for geoscientists.
Scope for Future Work
In the near future, as a follow-up to the current study, we wish to attempt the following:
(1) Based on the pierce-point distribution of the receiver functions, the receiver functions would be segregated according to their sampling variable geology in a region. This could enable us test the degree of correlation (e.g., Arianos and Carbone, 2009) of the signals and their origin, leading to better geologic/geophysical explanation of the results. For example, it would be interesting to test our method on data across some of the confirmed suture zones where two distinct geologic provinces are now in welded contact.
(2) Similarly, a correlation analysis of earthquakes recorded by several stations sampling the same geology but originating from different tectonic environments (e.g., convergence, rift-related, transform faulting, and plate interior regions) will be taken up. Such a study of events from diverse tectonic settings will enable us to characterize the similarities and dissimilarities between various source regions on the globe based on cluster entropy.
We are grateful to Pravin K. Gupta (Indian Institute of Technology, Roorkee, India) and Vinod K. Gaur (Indian Institute of Astrophysics, Bangalore, India) for clearing up concepts related to information theory and for their critical review of the manuscript. Rainer Kind (GeoForschungsZentrum, Potsdam, Germany) is thanked for critical reading of the manuscript. Suggestions by all above have improved the quality of this presentation. Thanks are due to the two anonymous reviewers and senior editor of the journal, Raymond M. Russo, who provided clear perception for improvement of our original submission. Help from our colleagues, R.K. Chadha for Indian data and P. Mandal for help with the synthetics, is acknowledged. This work received support from the Council of Scientific & Industrial Research (New Delhi) in the form of the National Geophysical Research Institute Supra Institutional Program (SIP) and in-house main laboratory project MLP-6509-28 (to Das Sharma).