Two-station Lg wave attenuation tomography in Eastern Asia

SUMMARY

1 INTRODUCTION

The Lg wave is a high-frequency seismic phase that can propagate efficiently over distant continental paths between ∼2° and ∼30° with a typical group velocity of ∼3.5 km s⁻¹ (Ewing et al. 1957; Nuttli 1973), making it the most prominent seismic phase in regional seismograms. The Lg wave can be modelled as the superposition of multiple supercritical reflected S waves within the crust (Bouchon 1982; Ou & Herrmann 1990) or the sum of higher-mode surface waves (Knopoff et al. 1973; Kennett 2002). The attenuation of Lg wave is quantified by the Lg attenuation coefficient or quality factor (Lg Q), which measures the amplitude decay caused by intrinsic and scattering energy losses. Since Lg waves sample the crustal waveguide they propagate through, Lg Q represents the average shear wave attenuation or Q of the crust.

Values of Lg Q are sensitive to factors such as subsurface temperature, extent of deformation, unconsolidated sediments and the presence of fluids or partial melt of crustal rocks (Mitchell 1995; Mitchell et al. 1997; Xie et al. 2004; Phillips & Stead 2008; Pasyanos et al. 2009; Zhao & Xie 2016). Furthermore, Lg Q values often exhibit a correlation with the length of time elapsed since the latest large-scale tectonic activity (Mitchell 1995; Mitchell et al. 1997). These properties make tomographic mapping of lateral variations in Lg Q a valuable tool for inferring both macroscopic variations of crustal structures and the microscopic rheology of crustal materials. In practice, such a model can also be used to recover source spectra, estimate source size, establish scaling of different source types and predict ground motion to support studies in earthquake and explosion sources, seismic hazard assessment and other relevant fields.

Eastern Asia is a vast area of highly diverse geology, comprised of several major terranes with distinctive crustal structures formed in a geologic time spanning more than 3 billion yr. Among those are the combined Tarim-Sino-Korean craton, the Kazakhstan–Kyrgyzstan block, the Siberian craton, the Yangtze craton, the Cathaysia block and the Tibetan plateau, as well as amalgamated microcontinents and island arcs, surrounded by Phanerozoic orogenic belts (Chen et al. 2007; Zhai et al. 2007; Pirajno 2013, Fig. 1). Volcano-sedimentary sequences and sedimentary basins are often superimposed on the existing tectonic units and occupy large areas across Asia. Three events of particular importance, the Permo-Triassic North China–South China collision, the destruction of the North China Craton (NCC) and the Cenozoic India–Eurasia collision, form the modern tectonic framework of Asia even though many details are still controversial. Therefore, investigating the seismic velocity, attenuation, thermal and rheologic structure in the crust and upper mantle is essential for fully understanding the lithospheric structure and evolution of the Asian continent.

Figure 1.

Map shows the major tectonic blocks and faults in eastern Asia. The blue lines delimit the major tectonic blocks in China. The white lines are the major faults. TNCO, Trans-North China Orogen; S-G, Songpan Ganzi block; SCFB, South China fold belt; WY, Wuyishan Terrane; N-Y, Nanling-Yunkai Terrane; NCP, North China Plain and HZP, Hongzhong Plain.

Open in new tabDownload slide

Numerous studies have explored Lg wave attenuation across Asia, spanning continental to regional scales (e.g. Mitchell et al. 1997; Fan & Lay 2002; Xie 2002; 2003a, b; Xie et al. 2004; Phillips et al. 2005; Pei et al. 2006; Xie et al. 2006; Hearn et al. 2008; 2008; Zhao et al. 2010; Bao et al. 2011; Zhou et al. 2011; 2013a, b; Ranasinghe et al. 2015; Ford et al. 2010; Singh et al. 2015; Thirunavukarasu et al. 2017; He et al. 2021). The Lg Q models from these studies have shown discrepant Q values, lateral variation patterns and resolutions, primarily due to the use of different data sets and Q estimation methods. Lg Q tomography presents unique challenges, making it more complex than, for instance, seismic traveltime tomography. One of the primary difficulties is the modelling errors dominant in the data used for Q inversion. These errors often stem from simplifications in the forward modelling process, such as using a 1-D geometrical spreading correction (Menke et al. 2006; Chen & Xie 2017). Furthermore, direct inversion of Q from Lg amplitude spectra forward function encounters an inherent trade-off between path attenuation and unknown source and site terms. This issue is analogous to the trade-off between velocity and event origin time encountered in traveltime tomography (Xie 1998; Menke et al. 2006). Effectively addressing these challenges is crucial for developing precise and accurate Lg Q tomographic models.

In this study, we use two-station based methods to measure interstation and interevent Lg Q values devoid of source and site effects. These methods include the two-station (TS, Xie & Mitchell 1990), reversed two-station (RTS, Chun et al. 1987) and reversed two-event (RTE) method (Bao et al. 2011). We address the issue of non-unity site response (SR) ratios in TS Q measurements and examine the statistical characteristics of modelling errors in Q measurements. A Singular Value Decomposition (SVD) based technique is adopted to tomographically image the laterally varying Lg Q and quantitively assess the resolution and errors of resulting models. Applying these methods to dense Lg wave data recorded in Eastern Asia, we present new two-station Lg Q tomography models in multiple frequency bands from 0.5 to 4.0 Hz. These models are of significant value for studies in relative geology and seismology and they serve as foundational tools and prior models for future Lg Q tomographic refinements with emerging seismic data.

2 SEISMIC DATA AND PROCESSING

This study utilized a total of 346 permanent broad-band seismic stations in Eastern Asia, including the Chinese National Digital Seismic Network (CNDSN) in China and IRIS stations deployed in Siberia, Mongolia, Kazakhstan, Kyrgyzstan, Tajikistan, Afghanistan, India, Thailand and the Korean peninsula (e.g. II, IU, IM, IN, KN, KR, KZ, RM, TJ and TM, Fig. 2). The CNDSN consists of a backbone network with 144 stations distributed throughout China and 31 provincial networks with more than 1000 stations. For this study, we selected the backbone and the regional stations in the Xinjiang, Qinghai and Tibet provinces, where the backbone stations are relatively sparse. Together, these selected stations provide good coverage of the study area. The seismic instruments installed in CNDSN include the BBVS series (Beijing Gangzhen Instrument & Equipment Co., Ltd), the CTS and JCZ-1 series (Wuhan Institute of Seismic Scientific Instruments Co., Ltd), the KS-2000 series (Geotech Instruments, LLC), the CMG series (Güralp Systems Limited) and the standard GSN STS-1 and STS-2 systems. These broad-band or wideband seismometers have flat responses to ground velocity from at least 0.1–10.0 Hz.

Figure 2.

Distribution of stations and seismic events. The stations are shown as triangles, with various colours representing individual networks. Earthquakes are shown as red dots.

Open in new tabDownload slide

To ensure homogenous distribution of events and stations and maximize interstation paths coverage over the study region, we chose 136 regional earthquakes with focal depths less than 40 km and M ≥ 4.5 that occurred in 2009–2011 and in 2015 (Fig. 2). We made use of events that occurred in Taiwan and were observed on the eastern coast of China. The low-frequency (≤1.0 Hz) Lg waves from these events may be weak but not totally blocked because they propagated across mixed oceanic and continental paths. Lg waves propagating within or across the Tibetan plateau are blocked (Ruzaikin et al. 1977; Xie et al. 2004). Because of high attenuation in this plateau, they can only be observed at 600–700 km epicentral distances (McNamara et al. 1996; Xie 2002). The parameters of the selected events are listed in Table S1.

Over 31 000 vertical-component seismograms were collected and analysed in this study. The Lg phase group velocity window between 3.6 and 3.0 km s⁻¹ was manually adjusted to avoid overlapping with the Sn coda and the beginning of the fundamental mode Rayleigh wave. The first P wave onset time of each seismogram was manually picked to minimize errors in the reported event origin time and hypocentral location. The waveform prior to the first P arrival was considered as the pre-event noise and the Sn coda defined within a group velocity window between 4.2 and 3.7 km s⁻¹ was considered as the pre-phase noise of Lg. Following Xie (1993), preliminary data processing included correcting the waveforms for instrument response, applying a 10 per cent cosine taper to each phase time window, and calculating the amplitude spectral density (square root of power spectral density) for the ground displacements of the Lg phase and noises. The pre-event SNR was estimated by dividing the Lg spectral amplitude by the pre-event spectral amplitude, and the pre-phase SNR was estimated by dividing the Lg spectral amplitude by that of the Sn coda. Only those Lg spectra whose pre-event SNR exceeds 3.0 and pre-phase SNR exceeds 1.3 were selected for further processing. To reduce the effect of ambient noise, we subtracted the pre-event noise power spectra from the Lg power spectra. The pre-phase noise was not subtracted from the Lg power spectra to avoid the underestimation of Lg amplitudes but was used as additional screening to ensure that the measured amplitudes are from the Lg phase instead of the Sn coda and avoid using seismograms in which Lg is absent because of blockage.

Fig. 3 illustrates the process using a featured event occurred in the Sichuan province, China, on 29 June 2009, at18:03:52 and M = 5.1. Fig. 3(a) shows the location of the stations relative to the event epicentre. Evident by the pronounced differences between the vertical component waveforms in Figs 3(b) and (c), Lg waves propagate well eastward but not westward because they are blocked by the Tibetan plateau. Furthermore, Figs 3(d) and (e) display the spectra of Lg, pre-event noise and pre-phase noise of two seismograms recorded at the westward station QML and the northeastward station LIF, respectively. It shows that the Lg spectra of both stations pass the pre-event SNR threshold (3.0), whereas at station QML the waveform is rejected because the Lg spectra relative to the pre-phase SNR threshold (1.3) fails the screening process at f ≥ 1.0 Hz. As Lg waves attenuate more over long distances at higher frequency (f ≥ 2.0 Hz), the number of available Lg spectral amplitude data is expected to decrease.

Figure 3.

(a) Map showing the example earthquake that occurred in Sichuan (yellow pentagon) and selected stations (red and blue triangles) recording the earthquake. The ray paths propagate westward (red) and northeastward (blue) respectively are also shown. (b) Vertical component seismograms of the stations on the westward paths. (c) Vertical component seismograms of the stations on the northeastward paths. The waveforms are sorted by epicentral distance. The onset of the Pn, Sn and Lg phases are marked on the seismograms as vertical bars. (d) Amplitude spectra of Lg wave (blue), pre-phase noise (red) and pre-event noise (green) recorded at station QML, one of the western stations. (e) The same as (d) but for station LIF, one of the northeastern stations.

Open in new tabDownload slide

The Lg amplitude spectra were estimated in nine frequency bands, 0.31–0.81, 0.54–1.04, 0.62–1.62, 1.08–2.08, 1.24–3.24, 1.69–3.69, 2.16–4.16, 2.64–4.64 and 3.12–5.12 Hz with central frequencies of 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5 and 4.0 Hz, respectively. The Lg Q tomographic models were generated at each frequency band individually, so that they can be directly used for relevant seismological studies without assuming a power-law dependence of Lg Q on frequency.

3 MEASURING INTERSTATION AND INTEREVENT PATH Q

3.1 Two-station based methods

The Lg spectral amplitude, |$A( {f,\Delta } )$|⁠, observed at an epicentral distance |${\rm{\Delta }}$| and a frequency f, can be stochastically modelled as (Xie 1993):

where |$S( f )$| represents the Lg source spectrum, |$G( {\rm{\Delta }} )$| the geometrical spreading (GS) and |$R( f )$| the site response (SR) at the receiver. The along-path attenuation is expressed as |$exp( { - \frac{{\pi f{\rm{\Delta }}}}{{VQ( f )}}} )$|⁠, where V is the Lg group velocity, assumed to be 3.5 km s⁻¹ in the Asian continent and |$Q( f )$| is the quality factor at frequency f. To estimate effective attenuation, the GS term must be first modelled and corrected from |$A( {f,\Delta } )$|⁠. Because the true GS is caused by an unknown 3-D crustal structure, we can only use a simplified GS from a 1-D structure with the form |$G( \Delta ) = \frac{1}{{{{{\rm{\Delta }}}_0}}}{{( {\frac{{{{{\rm{\Delta }}}_0}}}{{\rm{\Delta }}}} )}^m}$| for |${\rm{\Delta }} \ge {{{\rm{\Delta }}}_0}$|⁠, where |${{\Delta }_0}$| is a reference distance of 100 km within which Lg propagates as a body wave (⁠|$G( \Delta ) = {{{\rm{\Delta }}}^{ - 1}}$| for |${\rm{\Delta }} < {{{\rm{\Delta }}}_0}$|⁠) and m being a surface wave GS coefficient having a spectral value of 0.5 (Yang 2002).

The two-station (TS) method (Xie & Mitchell 1990) requires two stations, i and j, to record the same event in a geometry where the two stations and the event are aligned along the same great circle arc. The Q value along the interstation path can then be calculated from the ratio of the GS-corrected spectral amplitudes recorded at stations i and j:

where |${\rm{\tilde{\Delta }}} = ( {{{{\rm{\Delta }}}_j} - {{{\rm{\Delta }}}_i}} )$| is the effective distance between stations i and j with a ‘+’ or ‘−’ sign. In practice, we allow a small offset from the collineation (Fig. 4a) since the exact alignment of an event with two stations is extremely rare. We restrict the angular difference (⁠|${{\delta }_{az}}$|⁠) between the azimuths from the source to the two stations and the difference (⁠|${{\delta }_{baz}}$|⁠) between back azimuths within 15° (Bao et al. 2011; Gallegos et al. 2014; 2017). Meanwhile, we impose an additional constraint on the triangle formed by the event epicentre and two receiver sites, such that the relative difference in length between the two shorter sides and the longest side is limited to less than 1 per cent, that is |$( {{{L}_{\textit{shorter}1}} + {{L}_{\textit{shorter}2}} - {{L}_{\textit{longest}}}} )/{{L}_{\textit{longest}}}$| < 1 per cent. The minimum interstation distance is set to 150 km. These criteria have been established to ensure the collection of a greater number of TS paths while minimizing errors that could be introduced by a non-isotropic radiation pattern (Der et al. 1984; Xie et al. 2004), and to avoid improper geometry in which the two stations may be obliquely aligned when the event is located at approximately equal epicentral distances from both stations.

Figure 4.

Examples of actual recording geometry for (a) the two-station (TS) method, (b) the reverse two-station (RTS) method and (c) the reverse two-event (RTE) method. Stations are plotted in blue triangles, and sources are plotted in red stars. A small azimuthal difference |${{\delta }_{az}}$| (≤15°) and small back azimuthal difference |${{\delta }_{baz}}$| (≤15°) are permitted during data selection while preserving geometrical configuration.

Open in new tabDownload slide

The purpose of TS method is to cancel the common source spectrum, |$S( f )$|⁠, but it assumes a unity ratio of the SRs for stations i and j, that is |$\frac{{{{R}_i}( f )}}{{{{R}_j}( f )}} = 1$|⁠, in (2). In reality, the SRs of two distant stations are usually not the same, so that the non-unity SR ratio will raise bias in the estimated Q(f).

The reversed two-station (RTS) method developed by Chun et al. (1987) introduces a reversed two-station path from the opposite direction of a TS path using a second source (Fig. 4b). When eq. (2) is applied to source a and b, respectively, the Q value between stations i and j can be estimated from the cross-spectral ratio as:

where |${\rm{\tilde{\Delta }}} = [ {( {{{{\rm{\Delta }}}_{j,a}} - {{{\rm{\Delta }}}_{i,a}}} ) + ( {{{{\rm{\Delta }}}_{i,b}} - {{{\rm{\Delta }}}_{j,b}}} )} ]/2$| is the effective interstation distance with sign. Note that |$( {{{\Delta }_{j,a}} - {{\Delta }_{i,a}}} )$| and |$( {{{\Delta }_{i,b}} - {{\Delta }_{j,b}}} )$| must have the same sign in the geometrical configuration of RTS method. The major advantage of RTS over TS is that it is unnecessary to assume a unity SR ratio of the two stations in that the SRs have been cancelled in eq. (3) ultimately. The reversed two-event (RTE) method (Bao et al. 2011; Ranasinghe et al. 2015) is modified from the RTS. In the RTE geometry, the two sources that align with the two stations on the same great circle are situated within the station pair (Fig. 4c). The RTE method has the same expression as eq. (3), with |${\rm{\tilde{\Delta }}}$| now being defined as the effective interevent distance.

As the RTS and RTE methods effectively remove the source and site terms from the cross-spectral ratio, they can provide the most reliable path Q measurements. However, they require strict recording geometries, resulting in limited data coverage. In contrast, the TS method offers more flexible geometry and broader coverage, but the Q values are biased due to the presence of non-unity site response ratios in spectral ratios. We will address this issue by estimating SRs of individual stations in the next section.

3.2 Estimation of site responses and application to TS Q

The other application of the RTS method is to determine the SR ratio |$\frac{{{{R}_i}( f )}}{{{{R}_j}( f )}}$| of the two stations (Chun et al. 1987; Gallegos et al. 2017):

where |$\delta \tilde{\Delta } = [ {( {{{{\rm{\Delta }}}_{j,a}} - {{{\rm{\Delta }}}_{i,a}}} ) - ( {{{{\rm{\Delta }}}_{i,b}} - {{{\rm{\Delta }}}_{j,b}}} )} ]/2$| is small enough to be ignored in a proper colinear geometry (Fig. 4b). When we collected all SR ratios that sample different station pairs, we took the logarithm of eq. (4) to set up a linear system:

The matrix in eq. (5) has a dimension of |$( {N + 1} ) \times M$|⁠, where N is the number of measured SR ratios and M is the number of individual log-SRs to be inverted. The top N rows consist of equations with the form of |$ln( {{{R}_i}} ) - ln( {{{R}_j}} ) = {{b}_k}$|⁠, where |${{b}_k}$| is the logarithm of the left-hand side of eq. (4), i, j = 1 … M represent any station pair with SR ratio measurement and k = 1 … N is the index of SR ratios. Because absolute log-SRs are not resolvable from only the log-SR differentials, we assumed that the mean of log-SRs over all M stations is zero and imposed it in the last row in eq. (5) as a constraint for absolute log-SR values. In some past studies (e.g. Zhu & Chen 2012; Zhu 2014; Gallegos et al. 2017) that attempted to estimate SRs using the RTS method, some reference stations were chosen arbitrarily and were assumed to have zero log-SRs (i.e. no amplification or de-amplification). In our study, since the SR ratios are sampled from stations distributed in a large region, it is plausible to assume the mean of log-SR over all stations is zero when the near-station structure is lacked a prior. The only difference between those past and our study is that the resulting solutions in those studies are relative to their chosen reference stations, while our solutions are relative to the zero-mean of all log-SRs.

Eq. (5) can be solved as a damped least-squares problem. The precision of SR solutions is determined by the rank of the matrix in eq. (5), which increases with an increasingly larger number of SR ratio measurements, |${{b}_k}$|⁠. The maximum possible rank P is equal to M, which would allow a unique and precise solution for all M log-SRs. In this study, our dense and ubiquitous distributions for seismic stations and events allow the RTS geometry to connect most stations to form station pairs, yielding much more SR ratios measurements than the unknown SRs. As a result, eq. (5) tends to be very close to being full-ranked (e.g. N > M and P ≈ M). For example, at 1.0 Hz, we collected 5151 |${{b}_k}$| measurements that sample M = 293 stations out of the total 346 stations used in this study. Our calculation shows that the matrix in eq. (5) has a rank of 292, which is very close to 293. The log-SRs can thus be solved with very weak damping. We adopted the ‘L-curve’ criterion to determine the optimal damping parameter from candidate values between 0.1 and 10.0. For all frequencies we are interested in, the optimal damping parameters range from 0.7 to 1.2. The analysis for log-SR solutions of the above example is detailed in the Supporting Information.

The spatial distributions of the inverted SRs at central frequencies of 1.0, 2.0, 3.0 and 4.0 Hz are shown in Fig. 5. The distributions at other central frequencies are provided in Fig. S3. The relative SR values range widely, typically extending from 0.1 to 5.0 at each central frequency. Similar SR ranges have also been observed in different parts of the world (e.g. Drouet et al. 2008; Pasyanos et al. 2009; Bao et al. 2011; Zhu & Chen 2012; Zhu 2014; Gallegos et al. 2017). As the Lg wave is composed of multiple-mode crustal surface waves, local site response for Lg wave reflects a composite effect of velocity anomalies from the free surface to the Moho. This term also encompasses path-independent loss of energy near the station, potentially leading to de-amplification (Boore 2003).

Figure 5.

Geographic distribution of SRs estimated with the RTS method at central frequencies of 1.0, 2.0, 3.0 and 4.0 Hz, respectively. The colours represent SR values. The number of estimated SRs decreases with frequency due to the reduced number of Lg spectral data and RTS paths at higher frequencies.

Open in new tabDownload slide

In this study, we did not attempt to investigate the exact cause of SR variation and its relationship with tectonics and frequency. Instead, our focus was on utilizing individual SRs to rectify the SR ratios retained in the Lg spectral ratios of the standard TS method. To do that, we rewrote eq. (2) as:

where the non-unity SR ratio |$\frac{{{{R}_i}( f )}}{{{{R}_j}( f )}}$| is moved to the left-hand side and treated as a known quantity, leaving a pure Q free of the effects of source and site on the right-hand side. We called eq. (6) the site-corrected TS (SC-TS) method to distinguish it from the standard TS method expressed in eq. (2).

To illustrate the improvement in Q measurements using the new method, we compared the Q values measured by the SC-TS method with those measured by the RTS and standard TS methods for all common paths, respectively. Fig. 6 displays the distribution of the relative difference in Q at 1.0 Hz. On average, the relative difference between SC-TS Q and RTS Q is 10 per cent, while the relative difference between SC-TS Q and TS Q is 18 per cent. More than 75 per cent of the Q values estimated using the SC-TS and RTS methods have a relative difference of less than 10 per cent, whereas only 50 per cent of the Q values estimated using the SC-TS and TS methods have a relative difference of less than 10 per cent. At bins where the relative Q difference exceeds 10 per cent, the number of relative Q differences between SC-TS and TS is consistently greater than that between SC-TS and RTS. This comparison highlights that the difference between SC-TS and RTS Q measurements is much smaller than the difference between SC-TS and TS Q measurements, indicating that the SC-TS Q measurements are more reliable than the standard TS Q measurements.

Figure 6.

Distribution of the relative difference in Q values at 1.0 Hz between the Q estimated using RTS and SC-TS and between the Q estimated using TS and SC-TS for the common interstation paths.

Open in new tabDownload slide

While the RTS and RTE methods provide the most reliable Q estimates, their strictest geometry means that there will be far less ray-path coverage than the TS method with more flexible geometry. As our new SC-TS method provides Q measurements matching the precision of RTS and RTE Q measurements, we integrate SC-TS Q with RTS and RTE Q data to enhance coverage for the regions that cannot be sampled by RTS or RTE paths. The contribution of SC-TS Q is particularly significant for high-frequency bands (≥2.0 Hz). As illustrated in Fig. 7, SC-TS paths substantially expand coverage in areas such as eastern China, Mongolia and the Tibetan plateau, which are sparsely covered by RTS or RTE paths. It is important to note that SC-TS is applicable only to station pairs with SR solutions. Our data set indicates that a majority of stations forming TS paths have SR solutions, and merely about 5 per cent of TS paths are omitted in SC-TS estimations.

Figure 7.

Comparison of ray path coverage obtained from (a) the RTS and RTE methods alone and (b) the RTS, RTE and SC-TS methods. Path colours represent the average Q values at 2.0 Hz estimated for the paths, with high Q values (low attenuation) in blue and low Q values (high attenuation) in red.

Open in new tabDownload slide

4 DATA ERROR ANALYSIS AND SVD-BASED TOMOGRAPHIC METHOD

To conduct tomography using interstation and inter-event Q measurements, we rewrote all types of two-station Q estimates into a unified equation:

where |${{\tilde{A}}_e}( f )$| can be the left-hand side of eqs (2), (3) and (6) for the standard TS, RTS (or RTE) and SC-TS method, respectively, representing the measurements of actual Lg spectral ratios corrected for 1-D GS, |$\tilde{A}( f )$|⁠. |${{Q}_e}( f )$| specifies the estimated |$Q( f )$| from |$ln[ {{{{\tilde{A}}}_e}( f )} ]$| on the interstation (or interevent) path |$\tilde{\Delta }$|⁠. To tomographically map |${{Q}_e}( f )$|⁠, we discretized each path Q measurement in eq. (7) to form a linear equation:

where n = 1, 2, …, N denotes the nth ray path of the total N paths and m = 1, 2, …, M the mth model cell of the total M cells covering the study area. |${{Q}_n}( f )$| is the measured |${{Q}_e}( f )$| for the nth path, which has a length of |${{{\rm{\Delta }}}_n}$| and intersects the mth cell with a distance segment |${{{\rm{\Delta }}}_{nm}}$|⁠. |${{Q}_m}( f )$| is the unknown lateral variation of Q across the mth cell and is assumed to be constant. It should be noted that we invert |$1/{{Q}_m}$| rather than |${{Q}_m}$| values throughout the inversion process, although by convention we use the term ‘Q tomography’ and plot |${{Q}_m}$| values in tomographic maps.

4.1 Error structure in Q measurements

In Q tomographic inversion, the data contain notable modelling errors because the theoretical prediction in eq. (8) is not exact. Primarily, despite our correction for 1-D GS from observed Lg spectra, the effect of unmodelled 3-D GS persists in Q measurements. Additionally, as we assume isotropic Lg radiation and site effects, the effects of non-isotropic Lg radiation and azimuth-dependent SR remain in Q. Finally, the unmodelled amplitude noise and imperfect two-station recording geometry introduce errors into Q measurements. In contrast, the effects of ground motion amplitude measurement error are typically ignorable compared to modelling error. It is crucial to note that because Lg Q tomography operates with the logarithm of Lg amplitudes, termed ‘reduced amplitude data’ (Menke et al. 2006), it would be imprudent to simply presume a Gaussian distribution for the errors in reduced amplitude data.

In the pioneering work of Chen & Xie (2017), they justified that the errors in 1-Hz Lg Q closely follow a Gaussian distribution by rigorously analysing data collected from southeast China. It is interesting to see if the Gaussian distribution is a universal phenomenon for Lg Q data. We therefore extended the analysis of Chen & Xie (2017) to Lg Q data collected in different frequencies across Eastern Asia.

Chen & Xie (2017) described the error of the quantity |$\tilde{\Delta }/{{Q}_e}( f )$| in eq. (7) using a small random variable |$\varepsilon $| under the first-order approximation:

A single |$\varepsilon $| sample can be estimated from repeated Q measurements of the same interstation (or interevent) path and composes a data set of |$\varepsilon $| samples. We denote the |$\varepsilon $| collected using the RTS, RTE, standard TS and SC-TS methods as |${{\varepsilon }_{RTS}}$|⁠, |${{\varepsilon }_{RTE}}$|⁠, |${{\varepsilon }_{TS}}$| and |${{\varepsilon }_{SC - TS}}$|⁠, respectively. The |${{\varepsilon }_{RTS}}$| and |${{\varepsilon }_{RTE}}$| distributions at all nine frequencies from 0.5 to 4.0 Hz are shown in Figs S4 and S5, respectively. The distributions of these |$\varepsilon $| samples can all be characterized by a centre at zero, with long-tail outliers in both the positive and negative directions. The large outliers of the random data |$\varepsilon $| are caused by physically unreasonable |${{Q}_e}( f )$| measurements, that is very large or very small and even negative Q values. Following Chen & Xie (2017), we screened out these large outliers using the two-standard deviation screening criterion (Figs S4 and S5). After screening, the remaining |$\varepsilon $| samples strictly follow a Gaussian distribution. The best-fitting theoretical Gaussian probability density function (PDF) for the screened |$\varepsilon $| samples is plotted in red lines over the observed |$\varepsilon $| distribution in Figs S4 and S5 for each case. The corresponding cumulative distribution functions (CDF) for the best-fitted theoretical Gaussian PDF and the screened |$\varepsilon $| samples are plotted next to each PDF (Figs S4 and S5). For each case, the theoretical CDF and the observed CDF have visually identical shapes, with a cross-correlation coefficient between 0.9900 and 0.9999. As rationalized by Chen & Xie (2017), such a high cross-correlation coefficient ensures that each type of screened |$\varepsilon $| strictly follows a Gaussian distribution.

The same conclusion can be drawn for errors in standard TS and SC-TS Q measurements. But the comparison between |${{\varepsilon }_{TS}}$| and |${{\varepsilon }_{SC - TS}}$| shows that the distribution of |${{\varepsilon }_{SC - TS}}$| is narrower than that of |${{\varepsilon }_{TS}}$|⁠, meaning that the variance of |${{\varepsilon }_{SC - TS}}$| is less than the variance of |${{\varepsilon }_{TS}}$|⁠. The comparisons between the distributions of |${{\varepsilon }_{SC - TS}}$| and |${{\varepsilon }_{TS}}$| at all nine central frequencies can be found in Fig. S6. The variance reduction indicates that the error level in SC-TS Q measurements is lower than that in TS Q measurements. This improvement can be attributed to the SR ratio correction in SC-TS Q, which is the only difference between the two types of Q measurements. Ford et al. (2008) reported that the standard TS Q measurements are often less stable than the Q measurements obtained using other methods due to uncorrected site effects. We also found that more than 10 per cent of standard TS measurements were comprised of unreasonably small and negative Q values due to the uncorrected SR effects. In comparison, the proportion of unreasonable Q values decreased to less than 2 per cent in SC-TS measurements, a level similar to that found in RTS and RTE Q measurements. Once again, we demonstrated that the accuracy of the SC-TS Q measurements was on the same level as the RTS and RTE measurements.

The screening procedure removed less than 15 per cent of the |$\varepsilon $| samples in each case and accordingly removed the same percentage of Q measurements, most of which are repeated measurements. Therefore, the screening procedure can serve as the second round of data quality control. For any path with repeated Q measurements, their average is used as the final Q measurement. As most removed Q values are duplicates, the screening procedure does not significantly reduce the number of interstation (and interevent) paths. Table 1 shows the final number of Q measurements used in the tomographic modelling at 0.5–4.0 Hz.

Table 1.

Open in new tab

The number of two-station Lg Q measurements as function of frequency.

.	0.5 .	0.75 .	1.0 .	1.5 .	2.0 .	2.5 .	3.0 .	3.5 .	4.0 .
RTS	6909	5204	4813	2533	2115	1341	820	516	344
RTE	517	357	322	351	372	299	190	140	105
SC-TS	13 899	10 788	10 748	7 288	6799	4785	3433	2497	2107
Total	21 325	16 349	15 883	10 172	9286	6425	4443	3153	2256

.	0.5 .	0.75 .	1.0 .	1.5 .	2.0 .	2.5 .	3.0 .	3.5 .	4.0 .
RTS	6909	5204	4813	2533	2115	1341	820	516	344
RTE	517	357	322	351	372	299	190	140	105
SC-TS	13 899	10 788	10 748	7 288	6799	4785	3433	2497	2107
Total	21 325	16 349	15 883	10 172	9286	6425	4443	3153	2256

Table 1.

Open in new tab

The number of two-station Lg Q measurements as function of frequency.

.	0.5 .	0.75 .	1.0 .	1.5 .	2.0 .	2.5 .	3.0 .	3.5 .	4.0 .
RTS	6909	5204	4813	2533	2115	1341	820	516	344
RTE	517	357	322	351	372	299	190	140	105
SC-TS	13 899	10 788	10 748	7 288	6799	4785	3433	2497	2107
Total	21 325	16 349	15 883	10 172	9286	6425	4443	3153	2256

.	0.5 .	0.75 .	1.0 .	1.5 .	2.0 .	2.5 .	3.0 .	3.5 .	4.0 .
RTS	6909	5204	4813	2533	2115	1341	820	516	344
RTE	517	357	322	351	372	299	190	140	105
SC-TS	13 899	10 788	10 748	7 288	6799	4785	3433	2497	2107
Total	21 325	16 349	15 883	10 172	9286	6425	4443	3153	2256

The modelling error in the input data for Q tomography is then described by its variance, which can be obtained using the variance of |${{\varepsilon }_{RTS}}$|⁠, |${{\varepsilon }_{RTE}}$| and |${{\varepsilon }_{SC - TS}}$| at individual frequencies, respectively:

4.2 SVD-based tomography

When the Q measurements over N individual paths have been collected, all N equations with formula (8) build the linear system taking the matrix form of |${{\bf Gm}} = {{\bf d}}$|⁠, where the model vector m consists of the |$1/{{Q}_m}$| parameters to be inverted, the data vector d is determined via eq. (7) and G is a kernel matrix sparsely populated with the |${{{\rm{\Delta }}}_{nm}}$|⁠. To stabilize the inversion, we imposed Tikhonov regularization equations, |$\lambda {{\bf Lm}} = 0$|⁠, where |$\lambda $| is the smoothing parameter, and L is the 2-D Laplacian operator, to |${{\bf Gm}} = {{\bf d}}$| (Chen & Xie 2017).

Theoretically, SVD is a powerful tool for analysing the least-squares solution of |${{\bf Gm}} = {{\bf d}}$|⁠. However, its practical application is limited in large-scale data and modelling problems compared to iterative methods such as LSQR (Paige & Saunders 1982) and SIRT (Simultaneous Iterative Reconstruction Technique). In this study, we used the SVD-based 2-D tomographic method developed by Chen & Xie (2017). This method leverages the PROPACK software package (Larsen 1998) for efficient and accurate decomposition of the sparse kernel matrix into singular values and singular vectors with well-maintained orthogonality. As it is not an iterative technique, the inversion does not require any a priori initial model. The main virtue of this SVD-based tomographic technique is that it enables direct computation of model resolution and covariance matrices using inverses of singular values and vectors, allowing quantitative evaluation of the inverted model parameters.

5 RESULTS

5.1 Q_o tomographic model, resolution and error

Lg Q tomographic inversions were performed at frequency bands from 0.5 to 4.0 Hz individually. To ensure appropriate resolution, we divided the study area into 0.8° × 0.8° cells for frequencies up to 2.0 Hz and 1.0° × 1.0° cells for frequencies above 2.0 Hz based on the ray-path density. The inversion results for 1-Hz Lg Q (Q_o) are shown in Figs 8 and 9.

In assessing the spatial resolution of our Q tomographic model, we utilized the resolution matrix |${{{{\bf R}}}_{{\bf m}}}$| calculated in our SVD-based tomographic inversion. The resolution matrix quantifies how the estimated model spatially smears out the true structure. To present the resolution efficiently, we computed the discrete Backus-Gilbert spread function (SF) on each model cell (Menke 2012; Chen & Xie 2017) using |${{{{\bf R}}}_{{\bf m}}}$| as:

where |${{R}_{ij}}$| denotes non-zero elements in |${{{{\bf R}}}_{{\bf m}}}$| and |${{d}_{ij}}$| is the distance between the centres of cells i and j, with the summation extending over all j-th model cells adjacent to cell i. The SF thus measures the distance (km) over which each estimated 1/Q is smeared by the adjacent cells. A smaller SF value indicates better model resolution. For Lg Q_o tomography, the average SF is about 316 km excluding the Siberian craton with SF values larger than 1000 km (Fig. 8a). The average SF across entire China is about 250 km, and in regions with dense ray path coverage, such as the central-eastern parts of China and the Tianshan orogen (Fig. 8c), the SF values drop below 100 km. It is noteworthy that the discrete SF defined in (11) approximates the continuous SF concept introduced by Backus & Gilbert (1970), implying that the actual spreading distance might be smaller than the model cell size. For instance, the minimal SF for Q_o tomography is ∼48 m (Fig. 8a), smaller than the model cell size of 0.8° × 0.8°. This suggests that the model is under-parameterized in these areas given the chosen cell size.

Figure 8.

(a) The resolution of the model represented by Backus-Gilbert SF in kilometres calculated from the model resolution matrix. (b) The relative uncertainty of the model represented by δ(1/Q)/(1/Q) in percentage, where δ(1/Q) is one standard deviation of 1/Q calculated from the diagonal elements of the model covariance matrix. (c) Ray path coverage obtained from the RTS (blue), RTE (red) and SC-TS (green) Q measurements. (d) Checkerboard test with alternative Q values of 100 and 900 at 0.8° × 0.8° grids.

Open in new tabDownload slide

The covariance matrix |${{{{\bf C}}}_{{\bf m}}}$| quantifies how the error in the data propagates to the estimated model through a linear transformation (Aki & Richards 1980; Menke 2012). In this study, we show the first-order relative error of 1/Q, |$\frac{{\delta ( {1/Q} )}}{{1/Q}}$|⁠, which is calculated as the ratio of the square roots of the diagonal elements of |${{{{\bf C}}}_{{\bf m}}}$|⁠, that is the variance of the model error, over 1/Q (Fig. 8b). The off-diagonal elements, which quantify the correlation among errors in neighboring cells, are not shown here because their effects on model errors are secondary. The relative error of our 1/Q_o model remains below 5 per cent in most areas, including those with low resolution, and decreases to about 2 per cent in the well-resolved areas. This overall low error benefits from our rigorous screening for Q measurements.

Additionally, we conducted a checkerboard test using a 0.8° × 0.8° pattern with alternating high Q_o value of 900 and low Q_o of 100. The synthetic data set was generated using the actual ray paths (Fig. 8c) passing through the checkerboard model and was inverted following the same SVD-based tomography procedure. This test successfully reconstructs the pattern across China and Mongolia, but it is less effective at model borders, especially in the Siberian craton (Fig. 8d). The checkerboard test provides a somewhat idealized view of resolution capabilities by visually indicating areas that can reconstruct the pattern; however, it does not mean that these well-reconstructed areas have uniform resolution at 0.8° (∼89 km). The SF, on the other hand, quantifies the extent of smearing for each model cell, showing that it varies widely from ∼50 to ∼400 km, with an average of ∼250 km across different regions in China and Mongolia. The size of the checkerboard grid falls within the range of SF values, but the broader range of SF values reflects the varying density and the geometric distribution of ray paths, providing a more realistic representation of resolution. Although checkerboard tests are commonly utilized to evaluate the resolving power of tomographic methods and data coverage, the SF enables a more detailed assessment of the resolution in different areas, instead of an uniform resolving size provided in the checkerboard test. Consequently, we prefer to utilize the SF over checkerboard tests in this study for a more nuanced understanding of our tomographic results.

The Lg Q_o tomographic map derived using integrated RTS, RTE and SC-TS Q_o measurements is shown in Fig. 9(a). In comparison, Fig. 9(b) shows a model where the SC-TS Q_o data are replaced with standard TS Q_o data, which do not account for non-unity site response ratios. The comparison shows noticeable differences in the lateral variation patterns but highlights the impact of using different data sets on the resulting model parameter values. The model in Fig. 9(b) exhibits a broader range of Q_o values from ∼69 to ∼1000 in comparison with the range from ∼98 to ∼1000 in Fig. 9(a). Because the actual Q_o structure is unknown, and the tomographic inversion itself is influenced by data-driven regularization processes, this comparison is somewhat limited in its conclusiveness. However, the accuracy of Lg Q_o models is a critical factor for their application in seismological and geophysical research. The observed similarities between the Q_o measurements from RTS and SC-TS methods, associated with the reduction in errors in the SC-TS Q_o data, substantiate the enhanced accuracy and reliability of the SC-TS method.

Figure 9.

Lg Q tomographic map at a central frequency of 1.0 Hz using (a) RTS, RTE and SC-TS Q measurements and (b) RTS, RTE and TS Q measurements. Areas of high attenuation (low Q) are shown in red, and areas of low attenuation (high Q) are shown in blue.

Open in new tabDownload slide

In northeast China, low Q_o values are found in the Songliao (∼313), Hailar (∼293) and Erlian basins (∼307), while relatively higher Q_o values are observed in the Greater Khingan (Xing'an) range (∼411) and the Mongolian plateau (∼408). Lg Q_o varies considerably in the NCC, with Q_o as low as 196 in the east and as high as 613 in the west. In the South China Block (SCB), Lg Q_o in the Yangtze craton is lower than that in the Cathaysia block.

Notably low Lg Q_o values are observed in the Tibetan plateau. Its strong crustal attenuation greatly weakens and even blocks the development of Lg waves in the plateau (McNamara et al. 1996; Xie 2002; Xie et al. 2004), limiting the number of seismograms with a fully developed Lg phase. As a result, the average resolution of this plateau is SF = 340 km, lower than the average level of other regions in this study. The Lg Q_o image is poorly resolved in the Himalaya and the western part of Lhasa blocks with SF larger than 600 km.

Western China comprises major basins and active orogenic belts in China. The Junggar basin has a low Q_o of ∼305, while the Turpan basin has a relatively high Q_o of ∼669. The Tarim basin exhibits an alternating Lg Q_o pattern between high Q_o in the northeast and southwest margins and low Q_o in the middle. The Tianshan mountains have an average Q_o of ∼262, lower than the surrounding terranes, such as the Tarim, Turpan, Junggar and Kazakhstan.

5.2 Lg Q models at multiple frequencies

The Lg Q images at lower frequencies (<1.0 Hz) display similar lateral variation patterns to the Q_o image (Figs S7 and S8). At 2.0 Hz (Fig. 10), regions such as Songliao, Hailar, Bohai, Tibet, Qaidam, Tarim, Junggar, Turpan and Tianshan preserve the low-Q or high-Q features shown in Lg Q_o map. In NCC, the area of low Q extends westward, while the Lg Q values in the southeastern Yangtze craton and northwestern Cathaysia are lower than those in the Sichuan basin. In the middle of Tibetan plateau, a relatively high Q belt interludes through the persistent low Q region. As moving to even higher frequencies, the areas with well-resolved Lg Q become more localized (Figs 11 and 12, and Figs S9 to S11). For instance, the well-resolved areas at 4.0 Hz are confined to the Tianshan orogen, the central part of China and eastern China (Fig. 12). The zones with low Q (<750) are confined to Tianshan, Qaidam and central part of China. Conversely, relative high Q zones are observed in the Sichuan, Tarim and Turpan basins (∼1200) and the western NCC region (∼1232).

Figure 10.

(a) Lg Q tomographic map at a central frequency of 2.0 Hz. (b) The resolution of the model represented by Backus-Gilbert SF in kilometres calculated from the model resolution matrix. (c) The relative uncertainty of the model represented by δ(1/Q)/(1/Q) in percentage, where δ(1/Q) is one standard deviation of 1/Q calculated from the diagonal elements of the model covariance matrix.

Open in new tabDownload slide

Figure 11.

The same as Fig. 10 but for a central frequency of 3.0 Hz.

Open in new tabDownload slide

Figure 12.

The same as Fig. 10 but for a central frequency of 4.0 Hz.

Open in new tabDownload slide

To investigate the frequency dependence of Lg Q, we divided our tomographic Q images into distinct regions and analysed the behaviour of Lg Q at different frequencies within each area. We selected six major tectonic provinces within our study area: the Mongolian plateau (including the Greater Khingan range, Hailar basin and Erlian basin), Northeast China (including the Songliao basin and the Changbai mountains), NCC, SCB, Tibetan plateau and Western China (including the Tarim basin, Turpan basin, Junggar basin and the Tianshan mountains). Additionally, we identified specific tectonic blocks or units for analysis (Table 2). At each central frequency from 0.5 to 4.0 Hz, the tomographic Q values within each region or subregion were averaged to represent the Lg Q of that area. Fig. 13 illustrates the frequency dependence of Lg Q for these selected regions and subregions. Notably, as the coverage area of our Q tomographic maps decreases with the increase of frequency, Fig. 13 displays only those tectonic blocks that can still be covered by high-frequency tomographic maps.

Figure 13.

Lg Q frequency dependence for (a) the entire study area and six large regions discussed in the context and (b) selected typical tectonic blocks and units.

Open in new tabDownload slide

The frequency-dependent Q is commonly described in a power law (e.g. Mitchell 1995; Erickson et al. 2004):

where the power η describes how Q varies with frequency and is typically positive. While Fig. 13 indicates an overall increase in Lg Q values with frequency, it does not align well with the expected power-law relationship within the whole frequency range from 0.5 to 4.0 Hz. At 1.0 Hz, the change in Lg Q becomes steeper, dividing the curve into two segments that both exhibit a more linear trend in the double logarithmic plot (Fig. 13a). Consequently, we performed two separate fits of (12) to determine parameters Q_o and η: one using all data from 0.5 to 4.0 Hz and another using data solely from 1.0 to 4.0 Hz. The fitted slope η and the intercept Q_o, along with their associated errors represented by one standard deviation (the numbers in parentheses), are listed in Table 2.

When the parameters Q_o and η are determined using 1.0–4.0 Hz data for each of the six major tectonic provinces, the fitted Q_o values coincide with the averaged 1-Hz tomographic Q and the frequency dependence power η varies between ∼0.45 and ∼1.30 (Table 2). The η values for Mongolia, Northeast China, NCC, SCB and Western China range from 0.693 to 1.091, around the average η value (0.874) for the entire study area. The Western China has the lowest η and the best linear distribution of Lg Q and frequency from 0.5 to 4.0 Hz. While the η values for the major tectonic provinces are roughly consistent with each other, the Tibetan plateau stands out with a higher η compared to other regions (Fig. 13a).

For small tectonic blocks, the relationship between Lg Q and frequency is more complicated, especially at the lower and upper bounds of frequency, where Lg Q may decrease with increasing frequency (such as the Hailar basin, Erlian basin, Sichuan basin, Qilian mountains and Altai orogenic belts; Fig. 13b). It means that η can be negative in a local frequency range. The complexity of Lg Q frequency dependence for individual tectonic blocks is influenced by the specific characteristics of the local tectonic environment.

6 DISCUSSION

6.1 Tectonic features of Lg Q_o tomographic model

The lateral variation in our Lg Q_o map reveals contrasting features between different tectonic blocks in the well-resolved area. In NCC, the thermal destruction (decratonization) has divided this large craton into an eastern block (North China Plain) and a western block (Ordos), separated by the north-trending Trans-North China Orogen (TNCO or the Taihang block, Xu et al. 1994; Zhao et al. 2004; 2005; Kusky et al. 2007; Santosh 2010; Kusky 2011). The eastern block, marked by decratonization characteristics like high seismicity, high heat flow, thin lithosphere and Mesozoic to present-day volcanism, displays low Q_o values (∼196 in the North China plain). The western block exhibits higher Q_o values, with the Ordos thrust belt reaching ∼613, reflecting its cratonic nature and exposure of Precambrian rocks. The TNCO is distinguishable from both blocks with an intermediate value of ∼284, which could be attributed to the heterogeneities introduced by the large-scale thrusts, folds and strike-slip ductile zones in the Late Archaean–Proterozoic crust.

The SCB consists of the Yangtze craton, Cathaysia block and the Phanerozoic South China fold belt, and is separated from the NCC by the Qinling-Dabie orogenic belt. A laterally variable Q_o zone (∼235 to ∼412) coincides with the Qinling-Dabie orogenic belt, which is likely related to the intricate and largely unsolved mechanism of Qinling-Dabie collision and the cessation of the eastward push of the Tibetan plateau by the Yangtze craton. The Cathaysia block consistently shows higher Lg Q_o (∼430) values compared to the Yangtze craton (∼295). The changes in Q_o occur along the tectonic suture zone (the Jiangshan-Shaoxing suture or Jiangnan orogen, Pirajno 2013) that separates the Yangtze and Cathaysia blocks, which indicates different sources and evolutionary histories of these two blocks. The lower Lg Q_o of Yangtze than a typical craton could be attributed to Yanshanian tectono-thermal event (160 Ma), which led to the development of local-scale gentle folds and domes and the emplacement of volcanic rocks in the Yangtze craton (Su et al. 2009).

Within the Cathaysia block, Lg Q_o exhibits strong lateral variations, with high Q_o concentrated in the northeast and low Q_o in the southwest. Previous studies have identified two distinct tectonic terranes in the Cathaysia block separated by a roughly east–west boundary: the Wuyishan terrane in the northeast and the Nanling-Yunkai terrane in the southwest (Xu et al. 2007; Yu et al. 2010). The high Lg Q_o zone (∼591) correlates with the Wuyishan terrane, while the low Q_o zone (∼409) correlates with the Nanling-Yunkai terrane. These contrasting lateral variations in Lg Q_o very likely resulted from the different crustal evolutionary histories of the two terranes. The Nanling-Yunkai Terrane, composed of Meso-Neoproterozoic crust, suffered extensive deformation throughout the Caledonian (450 Ma), Indosinian (240 Ma) and Yanshanian (160 Ma) thermal events, while the Wuyishan Terrane consists of Palaeoproterozoic surface-exposed basement and was only reworked by Yanshanian magmatism.

The Tibetan plateau, formed by the ongoing India–Eurasia collision since ∼55 Ma, manifests the lowest average Lg Q_o value (∼153) out of our study area. The low Q_o zone starts at the Himalaya block, where the India Plate is being underthrusted northeastward beneath Eurasia, stopping at the Yangtze craton and NCC to the east and entering the Qilian fold belt and Qaidam basin to the north and the eastern Kunlun fold belt to the northwest. Tectonically, the Tibetan plateau consists of distinct blocks, such as the Himalaya, Lhasa, Qiangtang and Songpan-Ganzi, in an east–west trend, but the variation pattern of Lg Q_o does not differentiate between these tectonic blocks, probably due to limited resolution.

The abnormally low Q_o values, as observed in this and numerous previous studies (Fan & Lay 2002; Xie 2002; 2003a, b; Xie et al. 2004; Bao et al. 2011; Zhou et al. 2011; Zhao et al. 2013b; Thirunavukarasu et al. 2017; He et al. 2021), indicate the presence of high subsurface temperatures, melt-bearing materials and/or fluid content in the mid- to lower crust of the Tibetan plateau. Other geophysical anomalies, such as high Vp/Vs ratio, low-velocity and high crustal heat flow and electrical conductivities (Nelson et al. 1996; Owens & Zandt 1997; Wei et al. 2001; Klemperer 2006; Bai et al. 2010; Xin et al. 2019; Han et al. 2021; Liu et al. 2021), were also repeatedly found throughout the plateau and can be interpreted as unusually weak and even partial melting layers at different depths. Moreover, the petrological data suggested that the existence of crustal melting has a long history since the mid-Miocene (Hacker et al. 2000; Wang et al. 2012). Partial melting significantly weakens the crust and allows the flow of the mid- to lower crust under the gravity-driven lateral pressure gradients, thus making the crustal channel flow an acceptable mechanism responsible for the uplifting and lateral expansion of the Tibetan plateau (Nelson et al. 1996; Clark & Royden 2000; Clark et al. 2005; Royden et al. 2008).

In Clark & Royden (2000) model, the process of crustal channel flow is dominated by the northward and eastward lower crustal fluxes, which are finally halted by the rigid Tarim, Ordos and Yangtze blocks, respectively. The eastward flux bifurcates into northward and southward directions when resisted by the Sichuan basin, making the Qilian orogen in the north and the Sanjiang region in the south two possible escape paths of lower crustal flux. Low-velocity zones in the mid- to lower crust were found along the bifurcation paths (Zhang et al. 2013; Li et al. 2014; Liu et al. 2021). The shape of our low Q_o image in the eastern plateau appears to coincide with the bifurcation of this crustal flux and implies the flux intrusion into the eastern Kunlun, Qaidam basin and active Qilian orogen.

The India–Eurasia collision rejuvenated the Tianshan mountains, making it active uplift during the Neogene under the north–south compressional stress efficiently transmitted through the rigid Tarim continent. The modern Tianshan stretches about 2500 km long in east–west and 300–500 km wide in north–south, with elevations up to 7400 m. The western section of the Tianshan is largely located in the Kyrgyz Republic, and the eastern section is in China. Lg Q_o values in Tianshan vary from ∼199 to ∼359, generally lower than the Q_o of surrounding terranes, such as the Tarim, Turpan, Junggar and Kazakhstan. The low Lg Q_o beneath the Tianshan mountains is an indicative of partial melting and fluid presence, as well as small-scale convection or mantle upwelling beneath the mountain (Makeyeva et al. 1992; Roecker et al. 1993; Omuralieva et al. 2009). Under the collisional compression field, the molten material was forced upward, leading to the uplift of Tianshan and the development of active thrust faults spreading over the Tianshan range. These active thrust faults contribute to high seismic activity in this region and further reduce Lg Q_o due to scattering energy loss.

6.2 Major basins in China

Basins in China can be classified into two types: rift basins, which subsided due to lithosphere thinning and stretching and flexural basins, which subsided due to lithosphere downwarping under thrust loading in a compression regime (Watson et al. 1987). Rift basins, typically found in eastern China like the Songliao, Hailar, Erlian and Bohai basins are linked to the development of the Pacific oceanic margin. Conversely, flexural basins, including the Tarim, Turpan, Junggar and Qaidam basins in western China, correspond to north–south compression. Generally, the rift basins have lower Lg Q_o, while flexural basins show higher Lg Q_o. However, the Q_o values can be significantly influenced by post-formation mechanisms like sedimentary successions constrained by low velocities and high attenuation at shallow depths (0–10 km, Bao et al. 2015; Chen & Niu 2016; Shen et al. 2016; Xin et al. 2019; Han et al. 2021), mechanism transition and potential crustal flow.

The Songliao basin in Northeast China has an average Lg Q_o value of ∼313, lower than the hard and cold Greater Khingan range (∼411). This low Lg Q_o can be collectively attributed to several tectonic events: (1) The extensional stress-related basin rifting following intense volcanic eruptions in the Late Jurassic and Early Cretaceous: The associated lithospheric thinning and asthenospheric upwelling importantly heated the crust and resulted in the structural depression of the Songliao basin (Ren et al. 2002); (2) The hydrothermal reaction generated during the westward subduction of the Pacific Plate (Mitchell et al. 1997; 2008) and (3) The widespread Cretaceous sedimentary succession and a tectonic reversal, demonstrated as folding and reverse faulting in the Cenozoic or even upper Quaternary (Tang et al. 2009; Wang et al. 2012; Sun et al. 2013; Yu et al. 2015).

The lithospheric thinning and asthenospheric upwelling would be the main reason for the low Lg Q_o of rift basins in the eastern China, supported by high heat flow and geothermal gradients (Yi et al. 1992; Yuan 1996) and low-velocity layers at 20–40 km in the crust and 60–100 km within the mantle (Yuan 1996; Chen & Niu 2016). The presence of widespread sediments would also be responsible for reducing Lg Q_o because deposition took place in most of basins in China including some flexural basins. The Erlian and Hailar basins underwent strong Early Cretaceous rifting and have similar Lg Q_o values to the Songliao basin. The Bohai basin in the eastern block of NCC, formed under extensional stress associated with lithospheric thinning and asthenospheric upwelling and filled with a succession of sediments (Watson et al. 1987; Ren et al. 2002), exhibits an average low Lg Q_o value of ∼278.

The Tarim basin in western China, classified as a flexural basin (Watson et al. 1987), was formed by lithospheric down-warping when the Tarim crust was overthrusted by the Tianshan and western Kunlun orogenic belts. It displays a high Q_o zone (∼853) along the northeast and southwest margins, corresponding to the exposed Archaean basement rocks (3.2 and 2.8–2.6 Ga) of the Tarim craton. Reduced Q_o values (∼222) are observed in the centre due to extensive Triassic and Jurassic fluvial, lacustrine sediments and volcanic successions covering the basin.

The other two large basins in western China, Turpan and Junggar, present divergent Lg Q_o patterns. The Turpan basin, probably formed during the collisional events that created the Tianshan orogen (Windley et al. 1990; Allen et al. 1991), has an average high Lg Q_o (∼669). This overall high Q_o value suggests that the Turpan basin may not be heavily covered by thick sediments. In contrast, the Junggar basin manifests relatively low Q_o values ranging from ∼239 to ∼427. In the Carboniferous, the Junggar Plate collided with the Tarim and Siberian plates by consuming the Palaeo-Asian and Sayan Oceans that respectively separated them. As suggested by the geodynamic model in Zhou et al. (2008), the Junggar Plate transitioned to an extensional regime in the post-collisional period, associated with lithospheric thinning, rifting and crustal melting resulting from the upwelling of the asthenospheric mantle. This transitioned mechanism is likely responsible for reducing the Lg Q_o in this basin.

The Qaidam basin presents intriguing Lg Q_o variations, with high Q_o values up to ∼915 at its margins and low values Q_o down to ∼125 in the middle part. The formation of the modern Qaidam basin began with the structural depression at the onset of the India–Eurasia collision during the Palaeogene. The Precambrian–Silurian metamorphic rocks composing the basin basement are generally exposed along the margins of the basin. The central part of the basin is covered by Eocene to Quaternary sedimentary sequence with a total thickness over 15 km (Yin et al. 2008). The high Q_o corresponds to compressional tectonics and the exposure of ancient rocks at the margins, but the presence of thick sediment alone cannot fully explain the remarkably low Q_o value in the middle of the basin.

This low Q_o in the Qaidam basin is highly likely the result of the injection of a lower crustal flow, which probably originated from the Songpan-Ganzi lower crust, beneath stronger Qaidam crust (Karplus et al. 2011; Wang et al. 2012; Guo et al. 2017). This inference is supported by observations of a 15–20 km Moho offset between the thick northern Tibetan plateau crust (60–70 km) and the relatively thin Qaidam crust (∼45 km, Zhu & Helmberger 1998; Shi et al. 2009; Karplus et al. 2011). The crustal channel flow model proposed by Clark & Royden (2000) assumes that the crustal flow has different viscosities, determining the regional topographic slopes at margins of Tibet. The steep topographic contrast between the eastern Kunlun and Qaidam regions suggests a slow flux with relatively high viscosity. The injection of weaker Tibetan crust beneath the stronger Qaidam crust could serve as a mechanism explaining the active uplifting and northward growth of the eastern Kunlun Mountains into the Qaidam basin. If this model holds true, the pattern of low Lg Q_o could provide insights into the extent of flow expansion beneath the Qaidam crust.

The Sichuan basin in central China, which is classified as a flexural basin by Watson et al. (1987), exhibits a more complex behaviour while resisting the eastward extrusion of the Tibetan plateau. At 2.0 and 3.0 Hz, the Sichuan basin is characterized by distinctively high Q region compared to 1.0 Hz, highlighting its unique tectonic response.

6.3 Comparison with previous crustal attenuation and velocity studies

The lateral variation pattern in our Lg Q_o tomography generally aligns with earlier Q_o models for Eastern Asia (e.g. Phillips et al. 2005; Pei et al. 2006; Xie et al. 2006; Hearn et al. 2008; Ford et al. 2010; Zhao et al. 2010; Bao et al. 2011; 2013a, b; Ranasinghe et al. 2015), demonstrating a broad consensus on the regional Q_o structure. Notably, our model reveals more details associated with minor tectonic features.

Our Lg Q_o map for Northeast China and NCC shows patterns compatible with the models of Zhao et al. (2010), Zhao et al. (2013a) and Pei et al. (2006), respectively, yet the average Lg Q_o values in our model are constantly lower than what reported in those studies. The Lg Q_o values for Songliao, Erlian and Hailar basins (∼300) are consistent with the findings of Ranasinghe et al. (2015), where Lg Q_o values (∼200 to 400) were determined by applying the RTS method to an extensive collection of ray path measurements gathered by a very dense array deployed in Northeast China. Phillips et al. (2005) and Xie et al. (2006) also identified a higher Lg Q_o zone in Cathaysia compared to the Yangtze craton, but their models did not delineate the terrane division within Cathaysia.

There is a consensus among researchers about strong Lg attenuation in the Tibetan plateau (Fan & Lay 2002; Xie 2002; Xie et al. 2004; 2003a, b; Phillips et al. 2005; 2006; Singh et al. 2015; Bao et al. 2011; Zhao et al. 2013b; Thirunavukarasu et al. 2017), supporting the proposal by Fan & Lay (2003a) to lower the baseline for Tibetan Lg Q_o to ∼125. This value falls in the range of our Tibetan Lg Q_o variations (Table 2). Bao et al. (2011) and Zhao et al.(2013b) provided more detailed spatial variations in the Tibetan plateau, which aligns with some features found in our study such as the north–south bifurcation at the eastern plateau boundary. However, their models did not show the low Lg Q_o intrusion into the middle of Qaidam basin that our study reveals. Phillips et al. (2005) and Zhou et al. (2011) also identified Turpan as a high Q_o and Junggar as a low Q_o region. Significant Lg Q_o variations were found in Zhou et al. (2011) for the Tianshan and Tarim basin, but the detailed pattern differs from our results.

Despite the overall agreement with earlier models, we notice discrepancies in absolute Lg Q_o values, suggesting that while there is a shared understanding of the overall Q_o distribution, the specifics of the attenuation properties vary across different models. Such variations could be attributed to differences in methodologies, data collection, or interpretations of the underlying geophysical processes. For instance, joint inversions of Q, source spectra and/or site responses in Phillips et al. (2005), Zhao et al. (2010; 2013a, b) and Zhou et al. (2011) tend to increase the non-uniqueness of the results. Pei et al. (2006) and Hearn et al. (2008) utilized published M_L amplitude data for Q_o tomographic mapping, which also incorporates unknown terms controlling the source magnitude and station corrections. Additionally, how the bulletin data is related to the Lg spectral amplitude might be another factor contributing to the disparity in Q values.

Previous velocity models for Pn and Sn showed high-velocity zones beneath the flexural basins in western China (e.g. Liang et al. 2004; Sun et al. 2004; 2008). However, reading features such as northward crustal flow beneath the Qaidam crust from these velocity maps remains challenging due to resolution limits or the complex nature of the possible ‘double crust’ structure beneath Qaidam. Low upper mantle velocities were found beneath the rift basins with low Lg Q_o, but also observed in the Cathaysia with high Lg Q_o. The correlation between high crustal attenuation and low upper mantle velocities, and vice versa, may not be consistently observed across continents. It indicates the difference in seismic wave propagation characteristics between the mantle and the interior of the crust, as well as in the elastic and anelastic properties of the Earth.

6.4 Relationship of Lg Q with frequency, depth and crustal thickness

The frequency dependence of Lg Q for typical tectonic blocks in Eastern Asia is summarized in Table 2 and Fig. 13. For the entire study area and the six large tectonic provinces, variable rates of increase in Lg Q are observed at frequencies below and above 1.0 Hz (Fig. 13a), which were also noticed in previous studies (Zhao et al. 2010; Zhao et al. 2013b; He et al. 2021). The relationship of Lg Q with frequency is more complex in smaller tectonic blocks. Some tectonic blocks like the Hailar basin, Erlian basin, Sichuan basin, Qilian mountains and Altai orogenic belts (Fig. 13b), even exhibit a decrease in Lg Q as frequency increases, meaning a negative η within a local frequency range. Similar complexity has been observed in other regions worldwide, such as the Northern Arabian platform and Indian shield (Zor et al. 2007; Pasyanos et al. 2009) and the Colorado Plateau and California in Western US (Pasyanos 2013; Gallegos et al. 2017).

Despite the complexity of the frequency dependence of Lg Q in individual regions, we can classify these regions into three groups based on their Q_o and η values (Fig. 14). The first group includes rift basins in eastern China, eastern NCC, TNCO, Yangtze craton and Qinling-Dabie orogenic belts, which have average Q_o values below 300 and η values above 0.9. The Tibetan plateau, although it consists of distinct tectonic blocks, also falls in this group. Additionally, the Erlian basin, Changbai mountains and Sichuan basin are classified into the first group as their η values are above 0.9 although their Q_o values are higher than those of other regions. Conversely, the Wuyishan terrane in Cathaysia, Turpan basin and Altai orogenic belts constitute the second group with Q_o values above 500 and η values below 0.6. We classify the remaining regions into the third group, characterized by medium Q_o and η values. This inverse correlation between η and Q_o has also been observed in the Western US (Erickson et al. 2004; Gallegos et al. 2017), although η values may vary across studies. In summary, strong frequency dependence is generally associated with tectonically active regions with low Q values, while weak frequency dependence is generally associated with stable regions with high Q values (Campillo 1990).

Figure 14.

Relationship between Q_o and η for selected subregions. The Q_o and η values are obtained by fitting eq. (2) using Q samples at central frequencies from 1.0 to 4.0 Hz. The horizontal and vertical bars on each symbol indicate the standard deviation of Q_o and η, respectively.

Open in new tabDownload slide

Understanding the frequency dependence of Lg Q raises a fundamental question regarding its underlying physics. As Lg Q measures the depth-varying attenuation of crustal S wave (Q_s), two end-scenario causes can be considered: (1) the observed frequency dependence of Lg Q is a result of inherent frequency dependence of layered Q_s or (2) it is a manifestation of a depth-varying sampling weight of Lg assuming Q_s is frequency-independent (Mitchell 1991). When Lg waves are modelled as the summation of numerous high-order mode surface waves, hundreds of modes can participate in the contribution to Lg waves, but the number is dependent on frequency. It makes the relationship with depth and frequency much more complicated than the fundamental mode surface waves. As the frequency increases, more higher-order modes contribute incrementally to Lg summation (Kennett 2002). Although the lower-order modes tend to sample the shallow depth heavily with frequency increasing, the behaviour of additional higher-order modes is unclear. How the higher-order modes sample the crust, their relative weight in the summation, and whether high-frequency Lg waves tend to sample the shallow depth more heavily remains an area for future studies.

Zhao et al. (2010; 2013a) proposed a positive correlation between low-frequency Lg Q (0.2–1.0 Hz) and crustal thickness. The crustal structure across China exhibits a dichotomy, characterized by a thin crust in the east and a thick crust in the west, with a North–South transitional zone near 110°E that corresponds to Bouguer gravity anomalies and high seismicity (Liu et al. 2005; Kusky et al. 2007; Chen et al. 2010; Li et al. 2013). This transitional zone parallels to the front of the western Pacific subduction and roughly aligns with the Greater Khingan Range and TNCO. However, this dichotomy does not manifest in our Lg Q maps, as several exceptions disrupt this pattern. For instance, the Wuyishan terrane in the Cathaysia block has high Lg Q_o, but the crust is quite thin. The Tibetan plateau, with the thickest crust, shows the lowest Lg Q_o, which was also noticed in Zhao et al. (2013a, b) and thus this plateau was excluded in their analysis. The Tianshan mountains have a thicker crust compared to the Tarim basin (Chen et al. 2010; Li et al. 2013; Chen & Niu 2016), yet the mountains exhibit lower Lg Q_o than the basin. Our model suggests that there may be no intrinsic relationship between crustal attenuation and thickness. The positive correlation observed by Zhao et al. (2010; 2013a) in Northeast China and NCC should be a consequence of lithospheric thinning and asthenospheric upwelling under the extensional stress field, resulting in a thin crust and low Lg Q.

Finally, it is important to note that the estimated Lg Q encompasses both intrinsic and scattering attenuation effects. While the frequency dependence of intrinsic Q follows a power law, it remains unknown whether scattering Q follows a power law unless the heterogeneity is of fractal nature. Moreover, the contribution of scattering to Lg attenuation has been a subject of long-standing debate although strong small-scale heterogeneities, such as faults and folds, widely exist on the plateaus, orogenic belts and basins in Eastern Asia. Mitchell (1995) argued that scattering plays only a secondary role and may be negligible at low frequencies. Conversely, other studies (Dainty 1981; Campillo 1990) suggested a prominent role of scattering in crustal attenuation. Baumont et al. (1999) attributed the frequency dependence of Lg attenuation in the Bolivian Altiplano to small-scale heterogeneity rather than partial melting. Fan & Lay (2002; 2003a, b) showed that strong crustal heterogeneity scattering within the Tibetan plateau contributes to the strong Lg energy losses in this plateau. Addressing these questions extends beyond the scope of this study but warrants future investigations.

7 CONCLUSIONS

This study presents a comprehensive approach to enhancing the accuracy of Lg Q tomography in Eastern Asia. We introduced an improved method that overcomes the limitations of the standard TS method by correcting non-unity SR ratios in TS Q measurements using inverted SRs of individual stations. This advancement ensures our SC-TS Q measurements are devoid of site effects, yielding precision comparable to RTS and RTE Q measurements. The amalgamation of these refined measurements has facilitated the creation of an extensive and reliable data set for Q tomographic inversions, effectively merging the strengths of two-station methods with broader data applicability.

Our statistical analysis for modelling errors in path Q across a range of frequencies (0.5–4.0 Hz) reveals a general Gaussian error distribution, highlighting the stochastic nature of Lg spectral amplitude modelling. The primary sources of these errors include unmodelled 3-D geometrical spreading, non-isotropic Lg radiation, azimuthal SR variations, ambient noise and data sampling limitations.

Applying this improved methodology, we tomographically map the lateral variations in Lg Q across Eastern Asia using two-station Q data measured from over 31 000 broad-band Lg waveforms. Guided by the quantified model resolution and error, high Q_o (low attenuation) zones are found to correlate with regions of geological stability, ageing and low heat flow, such as the Greater Khingan, Ordos, Alxa and Cathaysia. Conversely, low Q_o (high attenuation) zones correspond to regions experiencing intensive deformation, volcanism and high heat flow, such as the eastern NCC, the Yangtze craton, the Tibetan plateau and the Tianshan orogen. Rift basins in China display low crustal Q_o, while flexural basins usually have high Q_o basements, but subject to reduction in Q_o values duo to post-formation mechanisms. Our analysis also delved into the frequency dependence of Lg Q, revealing complex variations that deviate from a simple power law, especially in smaller regions. An inverse correlation between η and Q_o values emerged, suggesting that regions with low Qo typically have high η. This pattern aligns with the idea that tectonically active areas tend to exhibit strong frequency-dependent behaviour in Lg Q.

In summary, this study advances the understanding of the crustal attenuation properties in Eastern Asia and its relation to tectonic structure and events. By incorporating sophisticated methodologies with an extensive data set, we provide the fundamental information that can be used in relevant seismological and geological studies.

ACKNOWLEDGEMENTS

This research was supported by the Air Force Research Laboratory contract FA9453-16-C-0046 and the United States Department of State contract SAQMMA17M1542 in responding to Broad Agency Announcement. We greatly thank Drs Z. Huang, R. Liu and other staff at CENC for providing data service and technical support. We also thank the handling Editors, Dr Thomas Hearn and the other anonymous reviewers for their insightful and constructive reviews that improved the manuscript significantly.

DATA AVAILABILITY

The waveforms recorded at the GSN stations were collected from the Incorporated Research Institutions for Seismology (IRIS) Data Management Center (DMC, https://ds.iris.edu/ds). The waveforms recorded by CNDSN were provided by the China Earthquake Network Center (CENC, https://news.ceic.ac.cn/). The PROPACK software package (Larsen 1998) is available at http://sun.stanford.edu/~rmunk/PROPACK/. The Lg Q tomographic models with resolution and covariance matrices in digital format and the input data can be downloaded from https://doi.org/10.6084/m9.figshare.25769004.v1.

REFERENCES