Lateral Variation of Crustal Lg Attenuation in Eastern North America

We perform QLg tomography for the northeastern part of North America. Vertical broadband seismograms of 473 crustal earthquakes recorded by 302 stations are processed to extract the Lg amplitude spectra. Tomographic inversions are independently conducted at 58 discrete frequencies distributed evenly in log space between 0.1 and 20.0 Hz. This relatively large dataset with good ray coverage allows us to image lateral variation of the crustal attenuation over the region. Obtained QLg maps at broadband and individual frequencies provide new insights into the crustal attenuation of the region and its relationship to geological structures and past tectonic activity in the area. The QLg shows more uniform values over the older, colder, and drier Canadian Shield, in contrast to higher variations in the younger margins. Results confirm the correlation of large-scale variations with crustal geological features in the area. Existence of low-velocity anomalies, thick sediments, volcanic rocks, and thin oceanic crust are potential sources of observed anomalies. The mean Q values are inversely correlated with average heat flow/generation for main geological provinces.

melting of deep mantle materials 47 . The dense volcanic material is responsible for a major gravity anomaly 48 . The presence of a significant mantle plume contributing to the MCR basalts implies a thinned lithosphere underneath the rift.
The Kapuskasing Uplift is a 500-km long zone of granulite and upper-amphibolite-facies rocks that transects the central Superior province. The sedimentary-volcanic belt of New Quebec (Labrador Trough) orogeny is located at the northeast margin of the Superior province. New Quebec is an 800-km-long thrust-fold belt. This orogeny was formed because of the Early Proterozoic collision between the Archean Superior and Rae provinces, while the Torngat orogeny on the western side is the zone of intense deformation formed as the result of the collision between the Rae and Nain provinces in the early Proterozoic. The Torngat orogeny is possibly younger than 1.81 Ga 38 . The Penokean orogeny at the southern edge of the Superior is a 1.90-1.83 Ga orogeny, crossed by the younger (1.1 Ga) MCR 46,48 .

Data and Pre-Processing
This study initially processed more than 100,000 high-gain vertical seismograms within the region recorded during the time-period January 1990 to January 2017. 1000 sec long vertical seismograms recorded on 642 stations within 1.5° to 30° epicentral distances have been automatically retrieved from IRIS data management center. Based on the catalog of events, 17-minutes long seismograms (~5 minutes before and 12 minutes after onset time) were cut from continuous data recorded by all stations in operation at the time of each event. Since the earthquake signals were not visible at all stations, we first used the STA/LTA technique 49 to distinguish traces containing earthquake signals from those recording just noise. Next, high-quality traces with a signal to noise ration (SNR) of equal or greater than two were used for further processes. All traces have sampling rates of 40 points per seconds and higher and have been checked and corrected for gaps and spikes. After removing the trends and mean, instrument responses were deconvolved using the causal correction method 50 .
Because the seismic signal and background noise are superimposed in the recorded seismograms, studies based on the measurement of observed amplitudes at different distances, such as attenuation estimation, can suffer from varying noise levels at different stations. To address this problem, Zhao et al. 51 proposed an approach to correct the amplitude of Lg-displacement in the frequency domain and reduce uncertainties in the Q estimation. However, in this study, we used time-frequency denoising techniques [52][53][54][55] to decrease the effect of background noise. These methods are data driven and automatically remove the noise from the entire seismogram and frequency bands. These methods assume that the seismic denoising can be treated as a non-linear nonparametric regression problem and use characteristics of sparse time-frequency transform to estimate the underlying regression function (seismic signal) from noisy observations. In these algorithms, noisy data is first transferred into the time-frequency domain. Then time-frequency coefficients in a pre-signal noise section (a few tens of seconds to a few minutes before the first arrival) are used to estimate noise level and find an optimal threshold level. Here, we applied wavelet transform using Morlet mother wavelet and to further constrain the internal arrival time estimation and selection of pre-signal noise we use an algorithm proposed by Kalkan 56 . Next, time-frequency coefficients are thresholded and inverse transformed into the time domain to reconstruct denoised signals.
The final dataset used for the tomographic inversion consists of approximately 20,000 high-quality waveforms recorded from 473 events on 302 broadband stations. All earthquakes are shallower than 36 km and have magnitudes larger than 2.5. Figure 2 presents the distribution of the associated events and stations. Figure 3 shows the plot record section of seismograms from one sample event before and after denoising.
After pre-processing all waveforms, the single-and two-station amplitudes were measured at discrete frequencies between 0.1 and 20.0 Hz. A group-velocity window of 3.6-3.0 km/s is used for the Lg extraction. The joint inversion is performed for individual frequencies.

Lg-wave Q Tomography
The attenuation of Lg-wave displacement is derived from the decay of spectral amplitude with epicentral distance (Δ), as a function of frequency f, based on the following relation 57 is the cumulative effects of other minor factors along the propagation path and computational errors.
is independent of frequency with a reference distance of Δ = km 100 0 58 . The geometrical decay rate γ is known to be 0.5 for the epicentral distance ranges used in this study Δ ≥ .°( 15 ), based on observational studies in eastern Canada 59,60 .
The attenuation factor in Equation (1) is expressed by:   Using this relation, the perturbations of the source S f log ( ) k , quality factor δQ x y f ( , , ), and site term δ P f log ( ) i are inverted in an iterative scheme from observed spectral amplitudes, assumed geometrical spreading, and the source, Q values and site term from initial models or previous iterations (denoted by superscript 0 ). Following 20 we constrain the site terms by assuming ∑ = = P f log ( ) 0 l l 1 .
To deal with the tradeoff between the source and attenuation terms and further constrain the model, the initial Q Lg is obtained using the two-station approach 7,21,22,57 , which improves the reliability of the estimated Q by removing the source term from the inversion model. In this method, the interstation amplitude is calculated from observed amplitudes of one event at two stations that are approximately aligned with each other and the earthquake epicenter 29,57,62 : where A kj and A ki are the spectral amplitudes recorded at stations j and i from event k respectively, and l is a reference point on the ray-path kj (see 63 for more details). Two systems of linear equations can be set up for single-and two-station data based on Equations (3) and (4) respectively: where δQ is a vector composed of the Q-perturbations, δU is a vector composed of the source-term perturbations, Matrix E sets up the relationship between the source-perturbations and the observed Lg-wave amplitudes, δP is a vector for the site-term perturbations, Matrix F is a bridge linking the site-perturbations with the observations, A s and A t are matrices setting up the relations between Q-perturbations and the Lg-wave spectral amplitudes based on single-and two-station relations data, respectively. H s and H t are vectors composed of residuals between the observed and synthetic Lg-spectra.
Zhao et al. 63 proposed a joint tomography approach by combining both dual-and single-station data to further reduce the attenuation-source tradeoffs and improve the resolution at high frequencies as: Based on the given assumption that the sum of all site responses is zero, we add an equation as and simultaneously controlling a relative variation of the site responses, we add: where ε is an empirical value for normalizing the site responses. The initial Q Lg is obtained by linear regression of two-station data; a unit source function is assumed as the initial source model. The LSQR algorithm with regularization, damping, and smoothing is used to solve the linear systems 16 . The damping is designed with a function of iterations as λ λ α = − k k 0 ( 1) , where λ 0 is the initial optimal damping coefficient, α is the attenuation coefficient, usually ranging from 0.5 to 0.9 and k is the iterative number. Here a smoothing with 5 points (a 2 nd order differential function) was selected for the smoothing. Perturbations δQ and δU are simultaneously inverted at individual frequencies by minimizing This hybrid approach has been successfully applied to map spatial variation of Q Lg in North China 63 , Tibetan Plateau 64 , and the Middle East 65 .

Results
Using the aforementioned inversion procedure, we obtained the crustal attenuation model of the region, consisting of Q Lg variations at 58 individual frequencies between 0.1 to 20 Hz. Q Lg maps at three selected frequencies (0.5 Hz, 1.0 Hz, and 3.0 Hz) and associated checkerboard resolution analysis 63 and ray coverage are presented in Fig. 4. Ray density is very high over most of the region. As can be seen from Fig. 4c,f,i, most of the area is covered by crossing ray paths providing relatively good azimuthal coverage. Among these frequencies, ray path coverage at frequencies lower than 1 Hz, is poor in western Superior Province. For the areas with better coverage of the two-station paths, the inverted Q would be independent of the input model and have certain values; whereas, for the regions covered with fewer two-station paths and more single-station paths, the uncertainty would be larger due to the possible source-Q and site-Q tradeoffs. The checkerboard resolution analysis shows that the maximum spatial resolution is higher than 2° in well-covered areas and for frequencies between 1.0 and 10.0 Hz.
From these maps, we see that Q values increase with frequency. The Canadian Shield (Superior province plus northern Grenville) is characterized by broad, uniformed low attenuation, while younger marginal areas that have experienced more intense deformation during past tectonic processes exhibit greater lateral variation and higher attenuation. Low-Q-anomalies (high attenuation) are observed at all frequencies in the New Quebec orogeny (anomaly #2 on the map with Q 0 ~= 600), coastal area of Labrador Sea (anomaly #3 on the map with Q 0 ~= 550), northwest of Hudson Strait (anomaly #1 on the map with Q 0 ~= 300), Belcher belt southeast of Hudson Bay (anomaly #11 on the map with Q 0 ~= 100-200), Wisconsin area south of Lake Superior (anomaly #10 on the map with Q 0 ~= 100-400), southwest of CHV (anomaly #7 on the map with Q 0 ~= 100), southeast of KP (anomaly #9 on the map with Q 0 ~= 100), southeast of SGL including the New York metropolitan area and Delaware Valley (anomaly #8 on the map with Q 0 ~= 100), southeast and northeast of Nova Scotia (anomalies #6 & #5 on the map with Q 0 ~= 300), and north of St. Lawrence Gulf (anomaly #4 on the map with Q 0 ~= 400). In contrast, the Great Lake Basin (anomaly #13 on the map with Q 0 ~= 2000) and NAP (including the state of Maine) (anomaly #14 on the map with Q 0 ~= 1500), exhibit high-Q-values (low attenuation).
The uncertainty in the tomography can be due to two sources. First, the final residuals, which could be used for measuring how many observed Lg amplitudes were not interpreted. The second issue is whether the observations were properly spread out between the source, Q and site terms, which usually depends on the survey system with the given ray coverage, including single-station and two-station paths, as well as some constraints for source 63 , and site 20 terms. In Fig. 5, the residuals and interpretation scales at 58 individual frequencies are provided. This plot shows that ~70% amplitudes (red dots) are interpreted by joint inversion including the source and Q, and ~80% (blue dots) can be explained if additionally considering site response. The residuals are relatively lower between 1 Hz to 10 Hz and reach ~20% which could be the random effects.
We also examined the mean and standard deviation of Q for all frequencies by resampling the original dataset using the bootstrapping technique. To obtain an error estimate, we randomly selected 85% of raypaths from the total number of observations to create 100 new bootstrap datasets and then inverted each bootstrap dataset to determine 100 Q values at each point. From these 100 Q values, we calculated a standard deviation and mean value for the region at each frequency. The directly inverted Q (Fig. 4) and Bootstrap Mean Q (Fig. 6) agree very well which indicate the stability of the results. The standard deviation is larger at very low (<1.0 Hz) and high (>13.0 Hz) frequencies. These agree with previous observations in the region 61 . High uncertainties at low frequencies can be due to the presence of microseismic noise and poor signals of the small events. While high uncertainties at high frequencies can be due to the stronger attenuation of Lg at higher frequencies or/and contamination of Lg window by the Sn coda.
Singh and Herrmann 66 conducted one of the earliest works addressing the crustal attenuation in southeast Canada, studying the coda attenuation of 250 local and near regional earthquakes with magnitudes between 3.0 and 5.0. They determined a crustal Q 0 (Q at 1.0 Hz) value of 700 to 900. Atkinson 71   Pasyanos 25 provided the attenuation models of crustal and upper mantle P and S waves in the North America for a 2-4 Hz passband. The high-Q anomalies in the eastern Great Lakes area and the low-Q anomaly south of Lake Superior that we observed in our study have been mapped 25 . He also estimated lower attenuation values for the cratonic Canadian Shield, compared to the younger surrounding areas. Gallegos et al. 34 performed a Q-tomography study using USArray data from 39 events (occurring between 2010 to 2012) using the two-station and the reverse two-station method, and provided 1.0-Hz-Q models for the central and eastern United States. However, since their study did not cover our study area, it is difficult to provide a comparison. Mitchell et al. 14 , who performed the latest attenuation study over the entire continent of North America, mapped the Lg coda Q 0 at 1 Hz. Mitchell et al. 14 obtained higher Q Lg over the Canadian Shield, in agreement with our results. Their model also identified relatively high Q in the eastern Great Lakes area and lower Q south of Lake Superior, consistent with our estimates.
In summary, estimates of crustal Q from previous studies in this region have a frequency dependent value of 0.19 to 0.65 and Q 0 varying from 525 to 1100 (Table 1), which are within the range of our results (Fig. 4). This variation in the previous results presumably comes from differences in the specific regions being analyzed, earthquake data with different magnitude and epicentral distance ranges, different frequency ranges, and/or the different methods used for Q estimation. The large dataset used in this study provides much higher ray path coverage, resulting in higher resolution of the model. We confirm the reliability of our results by comparing the range of our obtained Q values and verifying dominant anomalies with values from previous studies. In addition, finer attenuation structure obtained in this study can help improve our understanding of crustal structure of this region.
Correlations with geological structures. The age, physical properties, thermal status, degree of heterogeneity, and/or crustal thickness all affect Q Lg . Hence, crustal attenuation can reflect large-scale crustal features and the intensity of recent tectonic activity in a region 13,72 . Q Lg exhibits higher values for stable regions, in contrast to lower values in tectonically active regions. Q Lg values have a direct relationship with crustal thickness 62,63 and an inverse relationship with either strong scattering by small-scale heterogeneity or partial melting in the crust 22 . Thick unconsolidated sediments 12 volcanic and geothermal regions 6,27 have been associated with low-Q-values. We present our broadband Q Lg model (average Q values of individual frequencies between 0.5 to 5.0 Hz) in Fig. 7a. This average crustal model, which has similar features as the Q 0 model (Q Lg at 1 Hz), is used to further explore the correlation of the attenuation model with geological structures in this region. Generally, Q Lg increases with the age of a region since its last major tectonic event 65 , because crustal rocks cool over time. The Canadian Shield is the exposed part of the craton that is the oldest and most stable part of North America, which explains the relatively high and uniform Q values observed over the Canadian Shield. Moreover, several studies showed that seismic velocities are significantly higher than average under the Canadian Shield 44 , while the surrounding Paleozoic areas have slower S wave velocity 44 (Fig. 6d).
The Canadian Shield also has low average heat flow (Fig. 7e), indicating that the cratonic lithosphere must be thick and cold 73 . The intrinsic attenuation of shear waves is strongly correlated with the temperature of the Earth. Average crustal heat generation is estimated to be 0.6 μW/m 2 in the Canadian Shield 74 . The New Quebec orogeny, northwest of Hudson Strait, and the Belcher belt are parts of a passive margin, composed of Proterozoic volcanic materials. Hence, relatively lower Q values observed at these regions correlate with a slightly younger geological setting and major thrust-faulting in these areas.
Relatively low Q-values observed in coastal areas of Labrador Sea and Nova Scotia can be attributed to their younger age relative to the neighboring areas, and/or Lg-blockage and strong attenuation of Lg energy in oceanic regions 27 . Mousavi et al. 61 have reported Lg-blockage for southeastern Nova Scotia. However, these areas have relatively low shear wave velocities (Fig. 7d) and thick sedimentary layers (Fig. 7f). Mitchell and Hwang 12 mentioned the thick unconsolidated sediments as a source of high Lg-attenuation.
A marked increase in heat flow has been observed in the young Appalachian Province. Shapiro et al. 75 mapped high mantle heat-flow and thin lithosphere to the southwest of the CHV. Moreover, low shear wave velocities have been reported in this zone (Fig. 7d) by 76 and 44. These results can explain the observed low-Q anomaly to the southwest of the CHV. Volcanic rocks of Quaternary age are located north of St. Lawrence Gulf 76 . Xie 6 and Zor et al. 24 have reported low Q Lg values for volcanic and geothermal regions.
The low-Q anomaly south of Lake Superior, reported by both 25 and 34, is confined inside the 1.1-billion-year-old MCR, a lithospheric-scale feature. The lithosphere beneath the rift is different from the surrounding continental material because of thinning, modification, or removal of the continental lithosphere and its replacement by less-depleted mantle material 77 . The upper lithospheric layer is thinnest, with a thickness of about 50 km, at the MCR 45 The crust in the west arm of the MCR, thicker than its surroundings, is composed of sedimentary rocks underlain by layered volcanic rocks 47 . The middle crust in this zone, associated with low shear wave velocities 44,45,77,78 and higher temperatures 2002 , is interpreted to be influenced by the Great Meteor hotspot track 77 .
The low-Q anomaly at southeast of SGL (#8), including the New York metropolitan area and Delaware Valley, is close to the area of the intense low-velocity anomaly known as Northern Appalachian Anomaly (NAA) 44 . This low-velocity anomaly has been associated with the Great Meteor hotspot 79 , the combined effects of repeated rifting processes and northward extension of the hotspot-related Bermuda low-velocity channel 80 , and a small-scale asthenosphere upwelling 81 .
Anomalously high shear-wave velocity has been reported 44,78,80 for the crustal structure of the Great Lake basin in the northwestern SGL, where dominant high-Q anomalies are observed (#13). The thickness of the continental lithosphere in this zone is about 115 km beneath the Michigan basin 80 . Yuan et al. 80 proposed the existence of deep-rooted high-velocity blocks in this area representing the Proterozoic Gondwanian terranes of pan-African affinity, which were captured during the Rodinia formation but left behind after the opening of the Atlantic Ocean. A big portion of this high-Q anomaly is within the Granite-Rhyolite province (1.55-1.35 Ga) 43 that was formed through intracratonic magmatism. 25 and 34 obtained similar high-Q values in this region.
The frequency dependence of crustal attenuation. The Q Lg variations in different frequency bands sheds light on the features of crustal structure in response to varying depth ranges and wavelengths. Crustal attenuation as a function of frequency has been investigated along three lines, two almost parallel in NW-SE direction and one orthogonal to these two (Fig. 8). Two apparent low-Q Lg anomalies along profile-I correspond to the low-Q Lg under the Belcher belt (#11) and north of the Gulf of St. Lawrence (#4), respectively. The absorbing frequency band of the anomaly under the Belcher belt is relatively wider: 0.6-1.6 Hz, while it is narrower at 0.6-0.8 Hz under the NGSL region.
Three low-Q Lg anomalies along profile-II correspond to high attenuation under the east Kapuskasing Uplift (EKU), southwest of the CHV, and southern Nova Scotia (SNS), respectively. Dominant absorbing frequencies for all anomalies are less than 1 Hz with a decreasing trend from northwest to southeast toward the Atlantic coast. Anomalies along profile-III are associated with NGSL and southwest of the CHV. The high-Q anomaly along this profile has a relatively wide frequency band (0.3-2.2 Hz). We also observe relatively lower Q Lg within higher frequencies at the end of each profile where it crosses off-shore regions, which can be due to the high attenuation of Lg within oceanic crusts 61 .
Average Q. To investigate any systematic differences between attenuation in different geological provinces, we calculated the average Q Lg for major geological blocks, both at 1.0 Hz and at broadband (0.5-5.0 Hz). No clear relationship between the age and average Q value of geological provinces can be found ( Table 2), because of significant variation in Q values within each province. For instance, southern Grenville has a totally different range of Q values compared to the northeastern Grenville. However, an inverse relationship can be seen between average Q and heat flow/production (Fig. 9). Mitchell et al. 14 also showed that an inverse relation exists between crustal Q 0 and upper mantle temperature for North America. Curves in Fig. 9a-e exhibit a kink at 0.5 Hz that may result from a change between crustal trapped Lg and other surface waves (fundamental and/or higher modes together). Interestingly, the lowest frequency Q does not show much regional variation. In addition, the high frequency curve upward, that might be an effect of increasing Sn coda in higher bands.

Conclusion
High-resolution Q-models have been developed for eastern North America in the frequency range of 0.1 to 20.0 Hz and have been explained in the context of the tectonics and geology of the region. Our results not only agree with those from previous studies of this region, but also provide higher resolution and more reliable results due to higher azimuthal coverage of ray paths. Our Q estimates confirm the existence of heterogeneities in the crust related to the large-scale geological features in the area. Low-Q anomalies (high attenuation) are observed in New Quebec orogeny, coastal area of Labrador Sea, northwest of Hudson Strait, southeast of the Hudson Bay, south of Lake Superior, southwest of the CHV, southeast of Kapuskasing Uplift, New York metropolitan area and Delaware Valley, southeast and northeast of Nova Scotia, and north of St. Lawrence Gulf. In contrast, the Great Lake Basin and north Appalachian AP regions exhibit high-Q values (low attenuation). Thick sediments, volcanic rocks, and thin oceanic crust are potential sources of observed low-Q and low-velocity anomalies. The average Q of the main geological provinces is inversely correlated with average heat flow/generation.