TY - JOUR AU1 - Sanborn, Christopher J AU2 - Cormier, Vernon F AB - SUMMARY Comprised of S waves trapped in the Earth’s crust, the high-frequency (2–10 Hz) Lg wave is important to discriminating earthquakes from explosions by comparing its amplitude and waveform to those of Pg and Pn waves. Lateral variations in crustal structure, including variations in crustal thickness, intrinsic attenuation and scattering, affect the efficiency of Lg propagation and its consistency as a source discriminant at regional (200–1500 km) distances. To investigate the effects of laterally varying Earth structure on the efficiency of propagation of Lg and Pg, we apply a radiative transport algorithm to model complete, high-frequency (2–4 Hz), regional coda envelopes. The algorithm propagates packets of energy with ray theory through large-scale 3-D structure and includes stochastic effects of multiple scattering by small-scale heterogeneities within the large-scale structure. Source-radiation patterns are described by moment tensors. Seismograms of explosion and earthquake sources are synthesized in canonical models to predict effects on waveforms of paths crossing regions of crustal thinning (pull-apart basins and ocean/continent transitions) and thickening (collisional mountain belts). For paths crossing crustal thinning regions, Lg is amplified at receivers within the thinned region but strongly disrupted and attenuated at receivers beyond the thinned region. For paths crossing regions of crustal thickening, Lg amplitude is attenuated at receivers within the thickened region, but experiences little or no reduction in amplitude at receivers beyond the thickened region. The length of the Lg propagation within a thickened region and the complexity of over- and underthrust crustal layers can produce localized zones of Lg amplification or attenuation. Regions of intense scattering within laterally homogeneous models of the crust increase Lg attenuation but do not disrupt its coda shape. Coda waves, Earthquake monitoring and test-ban treaty verification, Seismic attenuation, Wave propagation, Wave scattering and diffraction 1 INTRODUCTION The Lg wave is often the largest feature in high-frequency (>2 Hz) seismograms observed at regional (200–1500 km) distance. Lg follows the Sn wave, initially arriving with a typical group velocity of 3.5 km s−1 and extends in a complex coda up to 100 s or more. Lg propagation can be represented either by sums of higher mode Love and Rayleigh waves confined to the Earth’s crust or by sums of multiple S waves critically reflected from the Moho and scattered by small-scale heterogeneities in the crust (Kennett 1986; Keers et al.1996). Since S waves are not radiated by idealized models of underground explosions, observations of Lg and the amplitude ratio of Lg to regional P waves (Pn and Pg) are important to discriminating the seismic signatures of underground nuclear tests from those of natural earthquakes. The efficiency of Lg propagation, measured by the peak amplitude of its coda, can be affected by lateral variations in crustal thickness, scattering and intrinsic attenuation. These propagation effects diminish the consistency of Lg/P amplitude ratios to separate earthquake from explosion populations. Because of its importance to source discrimination, the efficiency of Lg propagation has been extensively investigated, primarily by empirical observation by mapping and contouring path efficiency (e.g. Campillo 1987; Rodgers et al.1997; Furumura & Kennett 2001; McNamara & Walter 2001), but also by synthesizing waveforms using analytic, numerical and hybrid numerical/approximate methods applied to well-known and canonical crustal models (e.g. Campillo et al.1993; Wu & Wu 2001; Furumura et al.2014). General agreement among these studies is that zones of crustal thinning, for example, transitions from continental to oceanic crust, are associated with Lg blockage (Gregersen 1984; Mendi et al.1997). The efficiency of Lg propagation across zones of crustal thickening, for example, collisional mountain belts, however,is mixed, often as efficient as propagation across laterally homogeneous regions of continental crust, with blockage controlled by either intrinsic attenuation, scattering or by other large-scale complex structure (McNamara & Walter 2001; Sens-Schönfelder et al.2009). Fig. 1 summarizes these Lg behaviours from a few well-studied paths in Europe and Japan. Figure 1. Open in new tabDownload slide Example efficient (solid lines) and inefficient propagation paths of Lg in Europe and Japan assembled from Gregersen (1984), Campillo et al. (1993), Chazalon et al. (1993), Shapiro et al. (1996), Mendi et al. (1997), analysis by the authors of FDSN records of an earthquake in France (49.34°N, 6.62°E on 2003 February 22 at 20:41:04) and an analysis of a North Korean earthquake from G. Wagner 2017 (personal communication). Crustal thicknesses are from Grad et al. (2009) and the Cornell Moho (http://atlas.geo.cornell.edu/ctbt/geoid.html). Prior computational studies of Lg blockage in crustal thinning models include Maupin (1989), Cao & Muirhead (1993), Vaccari & Gregersen (1998), Zhang & Lay (1995), Mendi et al. (1997) and Wu & Wu (2001). The majority of these studies are limited to 2-D propagation and S waves polarized as SH. A radiative transport study that examined localized heterogeneity as an explanation for Lg blockage for paths crossing the western Pyrenees can be found in Sens-Schönfelder et al. (2009). Although 3-D, that study considered a relatively simple crustal model consisting of a layer over a half-space, locally perturbed by a velocity anomaly specified by a rectangular solid. Applying a similar radiative transport algorithm, this paper examines the effects of crustal thickness variations, scattering and intrinsic attenuation on the 3-D propagation of high-frequency Lg and Pg. The modelling is implemented by our Radiative3D code (Sanborn et al.2017), which includes the effects of multiple scattering in 3-D, mode conversion (P to S and S to P) and polarization changes induced by statistically described small-scale (≤ wavelength) heterogeneity. Our code also includes the effects of source radiation pattern, intrinsic attenuation, reflection/transmission by layered structures, generalized variations in crustal thickness, and the focusing and defocusing by large-scale deterministic structure. Since the North Sea Central Graben has been a well-studied path for Lg blockage investigations, the canonical models for crustal thinning effects are approximately scaled to the cross-sectional dimensions of the North Sea graben shown in Mendi et al. (1997) for a path from Norwich, UK, to station BLS near the Blåsjøreservoir in southern Norway. The dimensions of this structure are also close to those of the oceanic basin of the Ligurian Sea crossed by Lg paths from France to Corsica in the studies by Shapiro et al. (1996) and Contrucci et al. (2001). 2 TECHNICAL APPROACH 2.1 Modelling strategy The majority of previous modelling studies utilized 2-D (range and depth) numerical simulations. In order to compare with these prior studies with a 3-D method, we have chosen to make our models azimuthally symmetric, with the crust modification region (pinch or bulge) occupying a constant range from the source. This allows us to look at 2-D cross-sectional views of our model and wave-front snapshots by plotting range and depth only (collapsing out azimuth). The modelling thus simulates the propagation of seismic phases normal to the strike of a modified region, similar to previous 2-D studies, but includes the effects of out-of-plane scattered arrivals. Another advantage of the azimuthally symmetric model design is that, for a hypothesized strike-slip earthquake source, we can capture and isolate particular propagation modes based on the azimuth selected for a linear array of virtual seismometers. For a strike-slip mechanism the apportionment of radiated energy by polarization mode is a function of azimuth. For example, a seismometer array aligned with either the compressional or tensional principal axes of a strike-slip source will isolate a P–SV system of polarized waves, whereas an azimuth intermediate between the compressional and tensional axes will isolate SH polarized waves, as these azimuths coincide with the null-planes for P-wave radiation. Note, however, that for a strike-slip source, there are no null planes for S-wave radiation, and thus it is impossible to completely isolate P radiation. For this purpose, however, we can simulate a pure isotropic explosion source. Thus it is possible to produce simulations that illustrate the model effects on Pg in isolation from horizontally polarized Lg, and vice-versa. The gross structure of our models are fan-shaped to catch outward-radiating wave fronts within an azimuthal span of 90°. This is wide enough to enclose both an expanding SH wave front and an expanding P–SV system of wave fronts. The radiance maxima of these two systems have a 45° angle between them, so the 90° angular spread of the model means both azimuths of interest are well within the model bounds. We place three linear arrays of virtual seismometers. One array catches SH, one catches P–SV, and a third is intermediate between them and catches a mixture of all modes. The model extends to a range of 1020 km from the source and to a depth of 360 km below the Earth’s surface. A range of 1020 km is large enough that the Earth’s curvature can affect timing and amplitudes. As such, we include Earth-like curvature in the design of our model. The general modelling strategy is illustrated in Fig. 2. In all models except the baseline model, the wave fronts must cross a region of crustal variation (RCV), indicated by the brown band in the plan view (Fig. 2: top), and shown as a pinch region in the side profile view (Fig. 2: bottom). We consider RCVs in the form of crustal pinches, bulges or regions of high scattering due to heterogeneity. The layers of the model follow Earth-like curvature, as can be seen on the side profile. Figure 2. Open in new tabDownload slide Models are fan shaped (plan view, top) to accommodate outward radiating wave fronts within a 90° azimuthal arc from the source, which is indicated by a red dot in the handle of the fan. Choice of source parameters and seismometer arrays (illustrated by green discs representing gather area of every 10th seismometer) allows selective isolation or mixing of Lg and Pg phases. A crust variation region (pinch, bulge or high scatter region) begins at a range of 310 km from the source and concludes at 530 km. The model is azimuthally symmetric so that expanding wave fronts are perpendicularly incident on the structural modification. Model layers follow Earth-like curvature, as shown in the profile view (bottom). The event we simulate is alternately an earthquake or an explosion at a depth of 10 km. For the earthquake, we use strike, dip and rake parameters of 22.5°, 90° and 0°, respectively. This is a strike-slip source oriented so that the compressional and tensional principal axes are parallel to the ground and the null axis is oriented vertically. With these parameters, the compressional principal axis extends outwards at an azimuth of 67.5° east of north, which becomes a favoured direction for sampling the Pg phase. Offset 45° from that, at an azimuth of 112.5° east of north, is a maximum for SH. This becomes a favoured direction for sampling Lg. The radiation pattern is illustrated in Fig. 3. The strike orientation of the event is chosen to align with the virtual seismometer arrays such that one array catches SH polarized energy, which predominantly manifests itself as the Lg phase, and another array catches P polarized energy along with a small contribution of SV, which predominantly manifests itself as the Pg phase. Figure 3. Open in new tabDownload slide Radiation pattern of simulated earthquake source, showing directional dependence of polarization mode (P, SV or SH). The chosen strike-slip focal mechanism isolates SH radiation on azimuths of 22.5°, 112.5°, etc., which we use to situate seismometer arrays to observe Lg effects separate from Pg effects. Similar to classic beach ball diagrams, we plot polarization on the lower half focal sphere, using colour density to indicate composition of polarization mode. Solid colouration indicates outward first motion, whereas checkerboard colouration indicates inward first motion. 2.2 Models For the purposes of investigating different crust and bulge structures, we enumerate a set of model identifiers based on the labelling pattern NSCP nn, where nn is a numeric code, and the identifier can be read as North Sea Crust Pinch model, variant nn. (The model dimensions and positioning of the structural variations are loosely based on a path crossing the North Sea Central Graben.) We use NSCP 00 to refer to an unpinched, flat crust model, which can serve as a baseline to observe effects of pinches and bulges, which we will label with numeric codes 01 and above. The models are illustrated in Fig. 4 and layer thicknesses and attributes are tabulated in Tables 1–3. Figure 4. Open in new tabDownload slide Profiles of crustal thickness variations for variants of the North Sea Crust Pinch (NSCP) model series. Four models perturb the layer structure within a localized region (pinches and bulges), and two add a localized high-scattering region in the crust layer. Each model is divided into four major structural layers: a low-velocity sediment layer; a crust layer; a thin (5 km) high-gradient Moho transition layer; and the mantle, extending to a bottom of 360 km. The layers support velocity gradients, and reflection and transmission are simulated (including mode conversion) at layer interfaces due to velocity discontinuities. Table 1. Layer thicknesses in pinched/bulged regions of the NSCP models. NSCP 00 is the baseline condition in which the variation region is unchanged from the remainder of the model. For mantle upwell, negative values imply a model in which the Moho descends into the mantle, as in a mountain root or crust ‘bulge’. . NSCP 00 . NSCP 01 . NSCP 02 . NSCP 03 . Sediments 2 km 11 km 2 km 2 km Crust 30 km 11 km 11 km 50 km Moho transition 5 km 5 km 5 km 5 km Mantle upwell 0 km 10 km 19 km −20 km . NSCP 00 . NSCP 01 . NSCP 02 . NSCP 03 . Sediments 2 km 11 km 2 km 2 km Crust 30 km 11 km 11 km 50 km Moho transition 5 km 5 km 5 km 5 km Mantle upwell 0 km 10 km 19 km −20 km Open in new tab Table 1. Layer thicknesses in pinched/bulged regions of the NSCP models. NSCP 00 is the baseline condition in which the variation region is unchanged from the remainder of the model. For mantle upwell, negative values imply a model in which the Moho descends into the mantle, as in a mountain root or crust ‘bulge’. . NSCP 00 . NSCP 01 . NSCP 02 . NSCP 03 . Sediments 2 km 11 km 2 km 2 km Crust 30 km 11 km 11 km 50 km Moho transition 5 km 5 km 5 km 5 km Mantle upwell 0 km 10 km 19 km −20 km . NSCP 00 . NSCP 01 . NSCP 02 . NSCP 03 . Sediments 2 km 11 km 2 km 2 km Crust 30 km 11 km 11 km 50 km Moho transition 5 km 5 km 5 km 5 km Mantle upwell 0 km 10 km 19 km −20 km Open in new tab Table 2. Velocities (km s−1) and densities (gm cm−3) for each layer in NSCP models. Ranges imply a gradient from top to bottom of layer. . Density . P wave speed . S wave speed . Sediments 2.20 4.50 to 4.52 2.60 to 2.61 Crust 2.80 6.20 to 6.24 3.58 to 3.60 Moho transition 3.39 7.70 to 8.00 4.44 to 4.46 Mantle 3.40 8.0 and up 4.46 and up . Density . P wave speed . S wave speed . Sediments 2.20 4.50 to 4.52 2.60 to 2.61 Crust 2.80 6.20 to 6.24 3.58 to 3.60 Moho transition 3.39 7.70 to 8.00 4.44 to 4.46 Mantle 3.40 8.0 and up 4.46 and up Open in new tab Table 2. Velocities (km s−1) and densities (gm cm−3) for each layer in NSCP models. Ranges imply a gradient from top to bottom of layer. . Density . P wave speed . S wave speed . Sediments 2.20 4.50 to 4.52 2.60 to 2.61 Crust 2.80 6.20 to 6.24 3.58 to 3.60 Moho transition 3.39 7.70 to 8.00 4.44 to 4.46 Mantle 3.40 8.0 and up 4.46 and up . Density . P wave speed . S wave speed . Sediments 2.20 4.50 to 4.52 2.60 to 2.61 Crust 2.80 6.20 to 6.24 3.58 to 3.60 Moho transition 3.39 7.70 to 8.00 4.44 to 4.46 Mantle 3.40 8.0 and up 4.46 and up Open in new tab Table 3. Heterogeneity parameters and scattering statistics in model layers. Weak scattering is simulated throughout the model, but a few test conditions call for enhanced scattering within the zone of crust variation (the pinch/bulge/scatter region). Parameters ν, ε, a and κ (defined in the text) describe the von-Kármán heterogeneity spectrum of the medium. Mean free paths for P and S scattering at 2.0 Hz are shown and compared with the absorption lengths of intrinsic attenuation. An effective ‘scattering Q’ is computed from the mean free paths and compared with intrinsic Q’s. Scattering directionality is summarized by mean cosine of scattering angle, which we refer to as ‘dipole projection’ or DP. Absorption lengths and scattering Q’s are computed as follows: , and (and likewise for the P mode), where VP,S are the P or S elastic velocities and ω is the angular frequency of the simulated propagating wave. . ν . ε . a . κ . MFPP . MFPS . |$\Lambda _{\mathrm{P}}^{\mathrm{(absorb)}}$| . |$\Lambda _{\mathrm{S}}^{\mathrm{(absorb)}}$| . Sediments 0.8 1  per cent 0.2 km 0.2 8932 km 4632 km 161 km 41 km Crust (baseline) 0.8 1  per cent 0.2 km 0.3 13669 km 6802 km 1665 km 427 km Crust (enhanced scattering) 0.8 10  per cent 0.2 km 0.3 137 km 68 km 1665 km 427 km Moho transition 0.8 1  per cent 0.2 km 0.4 18774 km 8979 km 2073 km 530 km Upper mantle 0.8 1  per cent 0.2 km 0.5 16128 km 7364 km 1383 km 319 km (continued) |$Q_{\mathrm{P}}^{\mathrm{(scat.)}}$| |$Q_{\mathrm{S}}^{\mathrm{(scat.)}}$| |$Q_{\mathrm{P}}^{\mathrm{(int.)}}$| |$Q_{\mathrm{S}}^{\mathrm{(int.)}}$| DPP DPS Sediments 24943 22390 449 200 −0.087 0.180 Crust (baseline) 27704 23877 3374 1500 −0.186 0.060 Crust (enhanced scattering) 277 239 3374 1500 −0.186 0.060 Moho transition 30640 25412 3384 1500 −0.236 −0.007 Upper mantle 25334 20748 2172 900 −0.238 0.007 . ν . ε . a . κ . MFPP . MFPS . |$\Lambda _{\mathrm{P}}^{\mathrm{(absorb)}}$| . |$\Lambda _{\mathrm{S}}^{\mathrm{(absorb)}}$| . Sediments 0.8 1  per cent 0.2 km 0.2 8932 km 4632 km 161 km 41 km Crust (baseline) 0.8 1  per cent 0.2 km 0.3 13669 km 6802 km 1665 km 427 km Crust (enhanced scattering) 0.8 10  per cent 0.2 km 0.3 137 km 68 km 1665 km 427 km Moho transition 0.8 1  per cent 0.2 km 0.4 18774 km 8979 km 2073 km 530 km Upper mantle 0.8 1  per cent 0.2 km 0.5 16128 km 7364 km 1383 km 319 km (continued) |$Q_{\mathrm{P}}^{\mathrm{(scat.)}}$| |$Q_{\mathrm{S}}^{\mathrm{(scat.)}}$| |$Q_{\mathrm{P}}^{\mathrm{(int.)}}$| |$Q_{\mathrm{S}}^{\mathrm{(int.)}}$| DPP DPS Sediments 24943 22390 449 200 −0.087 0.180 Crust (baseline) 27704 23877 3374 1500 −0.186 0.060 Crust (enhanced scattering) 277 239 3374 1500 −0.186 0.060 Moho transition 30640 25412 3384 1500 −0.236 −0.007 Upper mantle 25334 20748 2172 900 −0.238 0.007 Open in new tab Table 3. Heterogeneity parameters and scattering statistics in model layers. Weak scattering is simulated throughout the model, but a few test conditions call for enhanced scattering within the zone of crust variation (the pinch/bulge/scatter region). Parameters ν, ε, a and κ (defined in the text) describe the von-Kármán heterogeneity spectrum of the medium. Mean free paths for P and S scattering at 2.0 Hz are shown and compared with the absorption lengths of intrinsic attenuation. An effective ‘scattering Q’ is computed from the mean free paths and compared with intrinsic Q’s. Scattering directionality is summarized by mean cosine of scattering angle, which we refer to as ‘dipole projection’ or DP. Absorption lengths and scattering Q’s are computed as follows: , and (and likewise for the P mode), where VP,S are the P or S elastic velocities and ω is the angular frequency of the simulated propagating wave. . ν . ε . a . κ . MFPP . MFPS . |$\Lambda _{\mathrm{P}}^{\mathrm{(absorb)}}$| . |$\Lambda _{\mathrm{S}}^{\mathrm{(absorb)}}$| . Sediments 0.8 1  per cent 0.2 km 0.2 8932 km 4632 km 161 km 41 km Crust (baseline) 0.8 1  per cent 0.2 km 0.3 13669 km 6802 km 1665 km 427 km Crust (enhanced scattering) 0.8 10  per cent 0.2 km 0.3 137 km 68 km 1665 km 427 km Moho transition 0.8 1  per cent 0.2 km 0.4 18774 km 8979 km 2073 km 530 km Upper mantle 0.8 1  per cent 0.2 km 0.5 16128 km 7364 km 1383 km 319 km (continued) |$Q_{\mathrm{P}}^{\mathrm{(scat.)}}$| |$Q_{\mathrm{S}}^{\mathrm{(scat.)}}$| |$Q_{\mathrm{P}}^{\mathrm{(int.)}}$| |$Q_{\mathrm{S}}^{\mathrm{(int.)}}$| DPP DPS Sediments 24943 22390 449 200 −0.087 0.180 Crust (baseline) 27704 23877 3374 1500 −0.186 0.060 Crust (enhanced scattering) 277 239 3374 1500 −0.186 0.060 Moho transition 30640 25412 3384 1500 −0.236 −0.007 Upper mantle 25334 20748 2172 900 −0.238 0.007 . ν . ε . a . κ . MFPP . MFPS . |$\Lambda _{\mathrm{P}}^{\mathrm{(absorb)}}$| . |$\Lambda _{\mathrm{S}}^{\mathrm{(absorb)}}$| . Sediments 0.8 1  per cent 0.2 km 0.2 8932 km 4632 km 161 km 41 km Crust (baseline) 0.8 1  per cent 0.2 km 0.3 13669 km 6802 km 1665 km 427 km Crust (enhanced scattering) 0.8 10  per cent 0.2 km 0.3 137 km 68 km 1665 km 427 km Moho transition 0.8 1  per cent 0.2 km 0.4 18774 km 8979 km 2073 km 530 km Upper mantle 0.8 1  per cent 0.2 km 0.5 16128 km 7364 km 1383 km 319 km (continued) |$Q_{\mathrm{P}}^{\mathrm{(scat.)}}$| |$Q_{\mathrm{S}}^{\mathrm{(scat.)}}$| |$Q_{\mathrm{P}}^{\mathrm{(int.)}}$| |$Q_{\mathrm{S}}^{\mathrm{(int.)}}$| DPP DPS Sediments 24943 22390 449 200 −0.087 0.180 Crust (baseline) 27704 23877 3374 1500 −0.186 0.060 Crust (enhanced scattering) 277 239 3374 1500 −0.186 0.060 Moho transition 30640 25412 3384 1500 −0.236 −0.007 Upper mantle 25334 20748 2172 900 −0.238 0.007 Open in new tab All model variants consist of seven structural layers. These are a sediment layer, a layer representing the solid bulk of the crust, a thin Moho transition layer, and four layers spanning the upper mantle down to a depth of 360 km. The four mantle layers follow AK-135 (Kennett et al.1995), and the three crust layers are loosely patterned after the crustal structure in the North Sea region. A velocity gradient in each layer is prescribed, such that velocity increases with depth, starting from some initial value, and ending at a final value, which may or may not be continuous with the next layer. Interlayer velocity discontinuities are utilized between the sediments and the crustal layer and between the lower crustal layer and the Moho transition layer. As in AK-135, a discontinuity is included at 80 km depth. All layers simulate scattering effects appropriate for a von Kármán heterogeneity spectrum. We utilize a fractional fluctuation ε of 1 per cent in local P and S velocities and a spectrum scale length corner a of 0.2 km in all layers to provide a weak background level of scattering throughout the model. The scattering formulation is that of Sato et al. (2012) and we specify also the density fluctuation Δρ/ρ = ν ΔVP/VP as a fraction ν of the P velocity fluctuation, as well as the Hurst parameter κ which controls rate of decay below the corner scale length. Scattering is simulated as a stochastic process controlled by angular probability distributions precomputed for each model cell. These distributions are dependent on the local average background velocities, and in layers with velocity gradients, the scattering behaviour is uniform throughout the layer and assumes the background velocity as specified at the top of the layer. Some models additionally feature a localized increase in the heterogeneity to 10 per cent fractional velocity fluctuation in the crust elements located within the structural variability region. Heterogeneity parameters and scattering statistics are tabulated in Table 3. 2.2.1 NSCP 01: Graben zone with sedimentary overlay Model variant NSCP 01 includes a pinched region with a sedimentary overlay and a slight mantle upwelling. In the pinch region, the lower crust layer tapers from an initial thickness of 30 km down to a final thickness of 11 km. The tapering is coupled with a broadening of the sediment layer from its initial thickness of 2 km to a final thickness of 11 km, and an upwelling of the mantle layer of 10 km. (The Moho transition layer thickness is unchanged; however the layer migrates upwards to remain intermediate between the ocean basalt and mantle layers.) In this configuration, the surface elevation remains fixed throughout the pinch region. The width of the fully pinched region is 100 km, and is bracketed on either side by taper regions of 60 km each, such that the complete width of the crust variation region is 220 km. The variation region begins at a range of 310 km from the source and concludes at 530 km. No changes are made to the velocities in variation region, only to the layer thicknesses. (The velocity gradients, however, are enhanced or diminished with the narrowing or broadening of the layers, so as to preserve the velocity values at the layerboundaries.) 2.2.2 NSCP 02: Thin crust due to mantle upwelling Model variant NSCP 02 includes a similar pinch structure to NSCP 01 (lower crust tapering from 30 km down to 11 km). This model, however, does not include a sedimentary overlay. (The sediment layer remains at a constant thickness of 2 km throughout the model.) Instead, the mantle upwelling is increased to 19 km. As in NSCP 01, the surface elevation remains unchanged throughout the variation region. 2.2.3 NSCP 03: Crust bulge protruding into mantle NSCP 03 is a crust bulge model. As in NSCP 01 and NSCP 02, we do not modify the elevation of the model surface, but instead allow the crust to bulge into the mantle. In the variation region, we thicken the crust from its initial thickness of 30 km to a final thickness of 50 km. Keeping the sediments thickness, Moho transition thickness, and surface elevation unchanged, this translates to a 20 km intrusion into the mantle. As in models NSCP 01 and NSCP 02, the variation region begins at 310 km and spans 220 km, which includes 100 km fully bulged span bracketed on either side by 60 km taperregions. 2.2.4 Models featuring scattering zones For each model structure, we additionally simulate a case where the heterogeneity in the region of crust variation (from 370  to 470 km, corresponding to the fully pinched or fully bulged regions) is increased from 1 per cent to 10 per cent fractional velocity fluctuation. In the case of structural baseline NSCP 00, in which there are no changes in layer thickness in the crust variation region, the enhanced scattering variation tests the effects of a high-scatter zone separate from the effects of structural variation. There is no additional heterogeneity (above the 1 per cent background levels present throughout the model) in the regions that correspond to the taper zones of the pinch/bulge models. The scattering zones are indicated by hatch markings in Fig. 4, and the heterogeneity parameters are listed in Table 3. 3 RESULTS AND DISCUSSION The simulation outputs primarily take the form of traveltime curves and energy curves. Both are representations of the signals collected by single-azimuth linear arrays of virtual seismometers, thus giving information about the wave fronts propagating along those azimuths. The traveltime curves are plots of energy amplitude observed at the surface as a function of range and time. Visual inspection of these plots allows easy identification of individual seismic phases and their relative strengths and temporal profiles. The energy curves illustrate time-integrated total energy recorded by the receivers as a function of range from the source. In this sense, they correspond to integrating the time axis on the traveltime plots, while retaining the range axis. As with the traveltime plots, these plots represent signal along a single azimuth as recorded by a single virtual seismometer array. The energy curves are one way to illustrate and quantify the effects of model structure on the outward propagating wave fronts. Towards this end, they are most useful when used with arrays that isolate a single seismic phase, such as Pg or Lg, as it enables us to observe the effects on those phases individually. 3.1 Wave-front time-series The time-series plot in Fig. 5 shows profile views (range and depth) of energy packets (hereafter referred to as phonons) propagating through model NSCP 01 as a function of time. This gives an indication of the bulk propagation characteristics of the wave-front propagation, and provides some context for interpreting the more quantitative traveltime and energy curves. In the figure, blue dots represent S-polarized phonons and red dots represent P-polarized phonons. The temporal snapshots are selected still frames from an animation that can be viewed by downloading files supplied in a supplement. Figure 5. Open in new tabDownload slide Earthquake time-series in crustal pinch model (NSCP 01). The time-series shows phonon propagation through a crust pinch earth model and illustrates how wave fronts evolve with time. Red markers represent P-phonons and blue markers represent S-phonons. 3.2 Traveltime curves Traveltime curves are presented in Figs 6–8 and indicate received energy flux as a function of range and time. Plotted image density is proportional to the square-root of the energy flux signal so that the visual presentation is amplitude-like, which better reveals coda detail to the eye. The energies are also scaled by a distance-dependent reference curve determined by a power-law fit to a reference model. The reference curve serves to normalize the image intensity so that phase structure remains easily visible despite range based decay, and in essence corresponds to the removal of the attenuating effects of geometric spreading and intrinsic attenuation. In fact, the decay constant of the reference curve can serve as an empirical measure of the combined effects of these attenuation mechanisms, and in this capacity is most useful in the baseline case where these mechanisms are laterally invariant. The reference curve also helps to highlight how one test case deviates from another. For example, in Fig. 6, the baseline (non-pinched, non-scattering) test case serves as the reference model both for itself and for the pinch model, so that the attenuation posterior to the pinch zone is made plainly visible in the pinch model by the clear reduction in image intensity. By contrast, in Figs 7 and 8, each model serves as its own reference, which enhances visibility of the detailed phase structure at long range, even though it makes the attenuation with respect to the baseline model less obvious. Figure 6. Open in new tabDownload slide Traveltime curves for the baseline model (left) and a crustal pinch model (right) showing disruption and attenuation of major phases resulting from crustal structure. Colour density indicates all-component energy amplitude (square root of energy flux across all channels) relative to a distance-dependent reference curve computed by power-law best fit to the time-integrated energy of the baseline model. The common reference curve (dashed line) and time-integrated energy signal (solid line) of each condition are inset at the top of the traveltime curves on a 4th-root scale to accommodate compressed vertical space. The region of crustal thickness variation is outlined by dashed vertical lines, and major regional phase velocities are indicated via velocity slope lines. Source event is an ideal strike-slip earthquake and the array azimuth is intermediate between the preferred azimuths for SH and P–SV radiation, such that Pg and Lg phases are both well represented in the images. Figure 7. Open in new tabDownload slide Crustal thickness and heterogeneity effects on all-component energy originating from Lg phase in a variety of Earth structures (cf. Fig. 4). Source is a hypothetical strike-slip earthquake, and energy signal is measured along an azimuth aligned with pure SH radiance. Image density is proportional energy amplitude (square root of energy signal) relative to a reference curve specific to each plot established by a power-law best-fit to the time integrated energy as a function of distance. Figure 8. Open in new tabDownload slide Crustal thickness and heterogeneity effects on all-component energy originating from Pg phase in a variety of Earth structures (cf. Fig. 4). Source is an isotropic explosion. Image density is a proportional energy amplitude (square root of energy signal) relative to a reference curve specific to each plot established by a power-law best fit to the time integrated energy as a function of distance. The energies we display in these plots are summations of all three components of motion, rather than restricting to a single component. This prevents the masking of converted phases which may otherwise fail to appear if they excite a different component of motion. Additionally, since our simulation reports energy flux, which is not identical with squared particle velocity without application of a free-surface correction (a planned but not-yet-implemented software feature), comparison of single-component signal to real data can be non-intuitive. In our experience the channel-summed energy provides more interpretive insight. Within each traveltime plot, the reference curve used for normalization and the time-integrated energy of the current plot are shown as dashed lines and solid lines, respectively, as an overlay plot. A fourth-root scale is used for this to fit narrow vertical space. Vertical lines outline the crustal variation region. Major regional phase velocities are indicated via velocity slope lines. In the example shown in Fig. 6, array azimuth is intermediate between the preferred azimuths for SH and P–SV radiation, such that Pg and Lg phases are both well represented in the images. Here we show only the baseline case and the pinch with sedimentary overlay case. Subsequent figures show the full set of test cases along azimuths selected to isolate Pg and Lg phases to facilitate separate analysis of Pg and Lg effects. The time-integrated energy curves for each model condition are shown on a logarithmic scale in Fig. 9, with separate plots for Lg and Pg effects, and facilitate a comparative analysis of the range-dependent decay of each model structure. These time-integrated energy curves are unwindowed about their respective phases, and thereby include energy from phases that convert from the incident Pg or Lg phase. Note however that the curves for one phase do not include cross-effects from the curve for the other phase, as the azimuthal isolation is sufficient to preclude this. Thus, the energy curves remain a good characterization of energy that began as either Pg or Lg. Figure 9. Open in new tabDownload slide Energy decay curves of Lg phase (left) and Pg phase (right) by crust structure variation. Shown here are whole-recording (unwindowed) time-integrated energies as a function of range from source. Lg phase is selected via choice of azimuth from an ideal strike-slip source. Pg phase is selected by use of a pure isotropic source. Unwindowed energies retain the contribution from signal diverted from originating phase into converted phases (e.g. Pn), which can be understood by reference to Figs 7 and 8. Crust structures investigated include: a flat-crust baseline; two structures that include crust pinch regions; and one that includes a crust bulge protruding into the mantle. A further variation on two crust profiles in which the variable region includes enhanced scattering via an increase in fractional velocity fluctuation from |$\epsilon =1\hbox{ per cent}$| to |$\epsilon =10\hbox{ per cent}$| are also presented. The red dashed lines show power-law best fits to the baseline conditions and afford an empirical measure of the effect of geometric spreading (albeit not isolated from effects of intrinsic attenuation) on energy. The curves fit energy to E(r) ∝ r−q with q = 2.12 in the Lg case and q = 2.51 in the Pg case. 3.3 Pinch versus Scattering Effects on Lg Here we examine the results from models NSCP 00 with scattering and NSCP 01 without scattering and compare and contrast the ways in which they differ from the baseline model. The two types of crust variation being considered here are a zone of high scattering and a zone of pinched crust with an overlaid sedimentary basin. In the next subsection, we compare the two pinch models (with and without sedimentary overlays) and the bulge model against each other. First, consider Figs 7(a)–(c). These show the traveltime plots along the azimuth favouring SH-polarized Lg for the baseline model (a), the pinch with sedimentary overlay model (b), and scattering-region model (c). This enables us to investigate the differing ways in which the crust pinch and scattering region disrupt the travel of the Lg phase through the variation region. Both test cases (b) and (c) are associated with a reduction of energy at long range compared with the baseline case (a). This reduction can be made quantitative by inspecting the energy curve comparison in the left panel of Fig. 9, which shows the time-integrated energy curves (surface energy flux as a function of distance from source event) for a selection of model conditions recorded along the Lg azimuth. In that figure, along with the quantitative summary in Table 4, we can see that case (b) (pinch with basin) resulted in an energy reduction of approximately −8 dB at the terminal distance of 950 km. Similarly case (c) (baseline plus scattering) was associated with a reduction of approximately −4.5 dB at the same distance. The greatest attenuation was achieved by the mixed pinch and scatter condition, a reduction of approximately −10.4 dB, showing the more-or-less additive quality of the two attenuation mechanisms. Terminal range attenuations are tabulated in Table 4. Table 4. Energy attenuation with respect to baseline of Lg and Pg phases at 950 km from source (420 km beyond conclusion of crust variability region), as summarized from Fig. 9. Energies are whole-recording (un-windowed) time-integrated three-component energies, and thus include contribution from converted phases and scattered energy resulting from interaction of the originating phase (Pg or Lg) with the model structure. . Lg . Pg . Baseline (NSCP 00) 0 dB 0 dB + Scatter (NSCP 00) −4.53 dB −1.49 dB Pinch with Basin (NSCP 01) −8.11 dB −4.73 dB + Scatter (NSCP 01) −10.38 dB −5.44 dB Pinch without Basin (NSCP 02) −2.99 dB −1.94 dB Bulge (NSCP 03) −1.54 dB +0.13 dB . Lg . Pg . Baseline (NSCP 00) 0 dB 0 dB + Scatter (NSCP 00) −4.53 dB −1.49 dB Pinch with Basin (NSCP 01) −8.11 dB −4.73 dB + Scatter (NSCP 01) −10.38 dB −5.44 dB Pinch without Basin (NSCP 02) −2.99 dB −1.94 dB Bulge (NSCP 03) −1.54 dB +0.13 dB Open in new tab Table 4. Energy attenuation with respect to baseline of Lg and Pg phases at 950 km from source (420 km beyond conclusion of crust variability region), as summarized from Fig. 9. Energies are whole-recording (un-windowed) time-integrated three-component energies, and thus include contribution from converted phases and scattered energy resulting from interaction of the originating phase (Pg or Lg) with the model structure. . Lg . Pg . Baseline (NSCP 00) 0 dB 0 dB + Scatter (NSCP 00) −4.53 dB −1.49 dB Pinch with Basin (NSCP 01) −8.11 dB −4.73 dB + Scatter (NSCP 01) −10.38 dB −5.44 dB Pinch without Basin (NSCP 02) −2.99 dB −1.94 dB Bulge (NSCP 03) −1.54 dB +0.13 dB . Lg . Pg . Baseline (NSCP 00) 0 dB 0 dB + Scatter (NSCP 00) −4.53 dB −1.49 dB Pinch with Basin (NSCP 01) −8.11 dB −4.73 dB + Scatter (NSCP 01) −10.38 dB −5.44 dB Pinch without Basin (NSCP 02) −2.99 dB −1.94 dB Bulge (NSCP 03) −1.54 dB +0.13 dB Open in new tab Although both the scattering region and the pinch-with-overlay test conditions resulted in significant attenuation on the far side of the variation region, there is a difference in the qualitative character of the attenuation, as can be seen again in the traveltime plots of Fig. 7. First, observe that in the baseline condition, the Lg phase is actually a superposition of many contributing multiple SmS paths as a result of the crust layer acting as an effective waveguide. This is seen in the plot as a striated grouping of many parallel phase curves propagating as a wave train within a broad velocity window surrounding the nominal Lg phase velocity. The banding structure is very cleanly visible in our simulated baseline traveltime plot due to several factors: (1) the relatively weak background heterogeneity spectrum of the bulk model, which makes the layers of our model a very clean propagation medium, in contrast to real Earth structure characterized by a rich spectrum of small-scale heterogeneity, (2) the idealizations of the ray-theory method of propagation at the heart of the radiative transport method, which does not simulate diffraction effects that might serve to wash out or blur the banding structure in real Earth propagation and (3) the precise regularity of the idealized, simplified and uniform layer structure of our model, which preserves the coherence of many multipaths due to the lack of fine lateral variations in the thicknesses of the layers or the bulk and shear moduli of the media. The result is an artificially clean-looking Lg wave train in the baseline traveltime plot. These plots, however, become useful in differentiating the disruption mechanisms of pinch versus scattering variation regions. In particular, we note in the scattering case that the energy amplitude is significantly attenuated, but that the multipath structure of the wave train and the width of the velocity window in which it arrives are largely unaffected. The scattering seems to result merely in a reduction of energy reaching the far side of the variation region, with little change in the character of the arriving signal. The mechanism of action here is almost certainly simple deflection of energy. The traveltime plot would seem to indicate a measurable degree of lateral scattering and back-scattering, as indicated by two distinct graphical features, to be described in what follows. Lateral scattering would include a deflection of the propagating phonons either left or right (w.r.t their propagation direction) or up-or-down. As we are probing an SH-polarized Lg phase, we would expect the left-right scattering to dominate, although an up-down fraction would exist to a lesser extent. In either case, both are loss mechanisms. The left-right scattering will divert energy away from the probed azimuth, which represented one of the four amplitude magnitude directional maxima in the SH radiation pattern, and into the Pg azimuths where SH radiation is at a minimum. Thus, more energy should be scattered away from this azimuth than into this azimuth, and the result will be an overall attenuation. The fraction of phonons that get scattered into up or down directions, on the other hand, will be lost due to other mechanisms. Those trending downwards may make it into the mantle, and thus be lost. Those trending upwards will reflect off the surface, and will then trend downwards. In the latter case, though, they will have interacted with the surface, and this should be visible. Indeed, in the distance range of about 340 km to about 460 km of Fig. 7(c), the Lg phase manifests a visible halo, appearing as orange in the colour scale used in the plot. We believe this may be the result of phonons being scattered into the surface before being reflected downwards and being lost. The second visible feature in the traveltime plot is a reflected phase originating at the onset of the scattering region at 370 km, and extending backwards with an Lg-like retrograde velocity in the time window between 120 and 200 s. The feature appears yellow on the plot and is faint, but visible. We believe that this represents back scattering. The mechanisms are different for condition (b) (pinch-with-overlay), however. In this condition, no abnormally high scattering is occurring in the crust variation region. In contrast to condition (a), wherein the phase structure is attenuated but not significantly disrupted, in condition (b) the banded multipath structure of the signal is severely disrupted in the variation region, and does not survive as a clear multipath signal on the far side of the region. There is still a clear Lg phase, and the timing of the phase is similar (the energy still falls within a very similar velocity window), but the signal has become amorphous in character. This must be a result of a perturbation of the waveguide character of the crust layer in the crust variation region. We can get some clue as to what is happening by again looking at the energy curves in the left panel of Fig. 9. The pinch-with-overlay condition is represented by the solid orange line, and the baseline condition by the solid blue line. The blue line follows a more-or-less uniform trend of decaying energy with increasing distance, as would be expected. The pinch line, however, experiences a major downward shift within the variation region and resumes a more gradual decaying trend after the region. An interesting aspect of the trend inside the variation region is that it begins with an initial increase in the energy interacting with the surface that starts slightly anterior to the onset of tapering leading to the pinch region, reaches a maximum as the tapering approaches the fully pinched region, then trends sharply downwards throughout the pinched region, and finally recovers some energy in a narrow region near the conclusion of crustal re-broadening and the region immediately posterior to the variation region. These three phenomena need explaining: the sharp decline of energy within the variation region, the initial amplification of energy, and the final amplification (recovery) of energy. Some possible explanations are suggested here, and in the next subsection we will attempt to corroborate them by investigating different pinch and bulge structures. The sharp decline of energy in the pinch region is most likely the result of phonon ray paths that transmit into the 11 km thick sedimentary basin overlaying the pinch region, wherein they become trapped until they exit back into the crust on the posterior side of the variation region. The sediments layer is a low-velocity, high intrinsic attenuation (high Q−1) zone, with also a shortened scattering mean free path (compared to the crust layer in the baseline model). As such, the energy content of each phonon will decay much more rapidly in the sediments layer than in the rest of the model. Because of the sharp velocity discontinuity between the sediments layer and the crust layer, down-trending phonons are very likely to reflect upwards, and multiple reflections between the surface and sediment-crust boundary mean that the sedimentary basin will be a waveguide for a slowly moving, multiply reflected, highly attenuating, Lg-like phase that will separate from the faster moving Lg phase still propagating through the narrowed crust. Another loss mechanism would be phonons that transmit from the crust layer into the mantle as a result of the reduced angle-of-incidence (which would tend to increase the transmission coefficient) that would occur where the crust floor angles upwards to produce the tapering region. The initial upswing in energy in the anterior tapering region is harder to explain, though two possible mechanisms may play a role. One is that the upward grading of the crust floor in the anterior tapering zone would bend rays reflecting off of it from a more glancing to a more directly upward trajectory towards the surface. Furthermore, as the ray transmits from the crust to the sediments layer, there is a decrease in velocity that would bend the ray towards the interface normal, which would be an additional adjustment towards vertically upwards propagation. These two bending events would have the effect of focusing the energy towards the surface, and might be enough to substantially increase the energy flux rate at the surface in that narrow region, before the attenuating effects of the sedimentary layer begin to dominate. The fact that the increase begins about 10 km prior to the onset of the tapering zone is puzzling, but may actually be a measurement artefact rather than a physical effect, as the virtual seismometers at 300 km range have a gather radius of 14 km, and thus are, in fact, sensitive to the leading edge of the crust variation region. Another mechanism that may explain the initial energy upswing is the possibility of rapid multiple reflections between the surface and sharp velocity discontinuity between the sediments layer and the crust layer, in the tapering zone where the sediments layer is still thin enough that the attenuation affect has not yet taken over. In essence, the energy might be rattling around in the sedimentary basin layer, serving to amplify the surface signal in the anterior regions, but ultimately succumbing to intrinsic attenuation as the sedimentary Lg fraction propagates towards the posterior region of the basin. The energy recovery that happens at the very end of the crust variation region (upswing in the pinch-with-basin series centred at about 540 km in Fig. 9, left panel) is perhaps the hardest to explain, and what is offered here is very conjectural. Again, two specific mechanisms are proposed. First, the rate of multiple reflections, or rattling, that occurs in the basin-shaped sedimentary layer, would increase as the sedimentary Lg fraction reached the tail end of the basin, due to the thinning of the layer back to its nominal thickness. This would increase the interaction of the phase with the surface, and is essentially the same mechanism as is proposed as a partial explanation of the upswing at the beginning of the basin. This effect would be observed in the region slightly before the conclusion of the crust variation region. It should not, however, contribute significantly beyond the conclusion of the region (at 530 km). The second mechanism may explain the remainder of the upswing, which occurs mostly beyond the conclusion of the variation region. Within the pinched-crust region, the Lg phase has separated into a slower-moving sedimentary Lg phase, trapped in the basin and the faster moving crustal Lg, where the thinned crust is still acting like a waveguide at depth. The energy trapped in this waveguide is not interacting with the surface. When the crust re-broadens, however, the energy expands back into the full-width crust and regains the ability to interact with the surface in the same way as it did prior to the crust variation region. So in this sense, it may be that the at-depth crustal Lg, which is shielded from the surface while trapped in the pinch, is returned to the surface after the pinch, and this could be recorded as an increase in surface energy flux in the region posterior to the crust variation region. In fact, a possible visible indicator of this can be seen in the traveltime Fig. 7(b) as an energy-dense temporally sharp pulse almost exactly aligned with the Lg velocity line that is present between 530 km and 680 km. 3.4 Pinch models compared To further explain the phenomenology of the pinch-with-overlay model, it is helpful to compare against other models of crustal thickness variation which differ structurally from the pinch-with-overlay model. To that end, we here compare models NSCP 01, NSCP 02 and NSCP 03, which are, respectively, (1) the pinch-with-overlay model, (2) a pinch model with the same degree of crustal narrowing, but without the sedimentary overlay (the mantle upwelling is exaggerated to compensate), and (3) a crust bulge model, in which crustal thickness increases and protrudes into the mantle, in a manner similar to a mountain root. Still keeping our focus on Lg phenomenology, the results we will be looking at here appear in Figs 7(b), (d) and (f), and in the left panel of Fig. 9. First, we will compare pinch results with and without sedimentary overlay. The qualitative differences in phase structure observed between Figs 7(b) and (d) are difficult to characterize, though they are visually distinct. It is easier to draw conclusions from the energy curves in Fig. 9, left panel. In this figure, the pinch-with-overlay results are represented by the orange line, and the pinch-without-overlay results are represented by the green line. The first thing to note is that the overall attenuation at the distal end of the model (950 km range) is less severe for the pinch-without-overlay model than for the pinch-with-overlay. In particular, we see that the pinch-without-overlay model experiences an energy reduction of approximately −3 dB with respect to the baseline model, compared with approximately −8 dB for the overlay model. This supports the notion that the sedimentary basin in NSCP 01 is a major component of the loss mechanism in that model, accounting, perhaps, for as much as 5 dB of attenuation. The remainder of the attenuation must be explainable by other mechanisms. The second thing to note in comparing the pinch models is that both pinch models show an initial upswing in surface energy flux in the anterior tapering wing of the crust variation region, and that both pinch models show a trend of decreasing energy throughout the remainder of the variation region, with the overlay model exhibiting the sharper downtrend, and the non-overlay model exhibiting a more gradual decrease in this region. Most interesting, though, is that whereas the overlay model exhibits a final upswing or recovery of energy at the completion of the crust variation region, the non-overlay model exhibits no such recovery, manifesting a smooth continuation from the posterior tapering wing of the variation region into the region beyond the crustal variation. In the preceding subsection, we proposed that rattling (rapid reverberation) in the very beginning and the tail end zone of the sedimentary basin played a role in the initial and final energy upswings. The absence of this reverberation in the model without the sedimentary overlay, and corresponding lack of final upswing in the energy curve, could be seen as support for the hypothesis that this plays a role in producing the upswing in the overlay model. A second mechanism, however, to explain the final upswing in the overlay model was also proposed. This second mechanism was the division of the Lg phase into two subphases each following a different propagation modality: a sedimentary Lg propagating slowly through the sediment basin (while rapidly attenuating), and the deeper crust Lg propagating with normal velocity, but isolated from interaction with the surface, until reaching the end of the basin where it could again interact with the surface, explaining the final upswing. This mechanism would also be absent from the non-overlay model, since the thinned crust in that model is positioned at the top of the model, rather than being pushed down to deeper depths by the presence of the sedimentary basin. Thus there would be no sudden energy return event at the conclusion of the crust variation region to produce an upswing of energy. Thus between both proposed mechanisms, we can say that both are likely contributing, though it is not clear without further investigation whether one mechanism is dominant, or whether both contribute significantly to the observed outcomes. The upswing is present in both models at the beginning of the variation region, however. Here again, two mechanisms were proposed. The first was basin rattling, and the second was a focusing effect from the upward inclination of the crust floor, which would bend rays reflecting off of the floor from a more horizontal orientation to a more vertical orientation, intensifying their interaction with the surface in that region. Basin-rattling is not an option in the non-overlay model. (Although a similar increase in multiple reflections is likely in play.) This leaves the change in incidence angle as the most viable explanation. Note that because the incline of the crust floor in the taper zone is more extreme in the non-overlay model, the intensifying effect could be more extreme. Indeed, the peak (at approximately 350 km) of the non-overlay energy curve is slightly higher than that of the overlay energy curve. It also occurs at a slightly more distant range than in the overlay curve. This could be because both mechanisms (rattling and focusing) are occurring in the overlay model, but only focusing is occurring in the non-overlay model, and the basin-rattling mechanism might be more efficacious in the nearer range window. The next thing to examine is what happens to the non-overlay energy curve within the fully pinched central region of the crust variation region. As is the case with the overlay model, here we see a generally downwards trend. Unlike the overlay model, however, the downward trend is both (1) not as steep as the overlay model and (2) the energy remains above the baseline curve until the very end of the fully pinched zone (the beginning of the posterior tapering wing). In fact, except for edge effects, the downward trend seems to parallel the baseline within the pinched zone. This may indicate that no additional attenuation mechanisms are in play beyond those also present in the baseline: namely geometric spreading, intrinsic attenuation and a certain amount of continual leakage into the mantle. The fact that the energy level rises above the baseline in this zone may be due to an increase in the rate of surface reflections on account of the shallower crust floor. This is similar to the rattling mechanism of the basin, only it is in basaltic rock media, not lossy sedimentary media. At the conclusion of the pinched region, as the crust begins to broaden back to baseline thickness, the energy level drops again. Posterior to the crust variation region, the non-overlay energy curve remains below the baseline curve, but above the pinch-with-overlay curve. This indicates that while the pinch-without-overlay had an overall attenuating effect on the Lg phase, it was not nearly as efficient at blocking the Lg as the model with the sedimentary overlay was. This makes sense, since a significant amount of the Lg-blockage mechanism of the overlay model was attributed to the intrinsic attenuation acting on the fraction of the Lg wave that propagated through the sedimentary basin. The remainder of the attenuation was attributed to mantle leakage as the phase was incident on the anterior tapering wing of the crust modification region, and that mechanism is probably the dominant loss mechanism in the non-overlay model. 3.5 The bulge model Last, we turn our attention to the bulge model (NSCP 03). The most striking thing to note is that this is the least disruptive of all the crust variation models. The bulge model traveltime results for Lg propagation appear in Fig. 7(f). At long range, the phase timing, structure and intensity are nearly identical to the baseline, except for a few small pockets of anomalous quiescence. The most striking difference occurs in the middle zone of the crust modification region, (the zone where the crust has expanded to its maximal thickness between the two tapering wings). Here, there is a very sizable quiescence with almost a complete vanishing of energy interacting with the surface. In the energy curve in Fig. 9, the plum-coloured line represents the bulge result. This figure shows that in the regions beyond the variation region, the bulge result parallels the baseline very closely and manifests only a very minimal attenuation of approximately −1 dB at the terminus of the plot. This tells us that the bulge structure was a very inefficient mechanism for Lg blockage. The majority of the energy propagates through the bulge and retains the majority of its structural and timing characteristics as well. Within the variation region, a very sharp dip, or notch, is seen in the energy curve. This corresponds with the quiescence zone in the traveltime plot. The mechanism for this is not hard to understand. Since the Lg phase is composed of a wave train of several multiply reflected wave fronts channelled between the surface and the crust floor, there is a substantial amount of interaction with the surface in the region anterior to the variation region. When the Lg phase enters the variation region and the crust thickness expands, the individual wave fronts, after reflecting off the surface, must bottom at a deeper depth before they can return to the surface. The extra time needed to bottom at the deeper depth creates a time window in which the energy is unavailable to interact with the surface. The length of the bulge region is such that when the energy returns to the surface, it has traversed the variation region and continues its original propagation modality. This would seem to implicate the length of the bulge region as an important parameter in controlling the signal in the bulge region, and perhaps beyond. For example, if the bulge region were sufficiently lengthened, it would be possible for the Lg phase to bottom twice before traversing the region. This should then produce two quiescent zones in the bulge region on the traveltime plot, with an active zone in between, and correspondingly two notches in the energy curve. It is interesting that the bulge model produced very little terminal-range attenuation, in contrast to the pinch models, especially since mantle leakage as the Lg phase was incident on the inclined crust floor was proposed as a significant attenuation mechanism. An inclined crust floor is not absent from the bulge model, but rather occurs at the end of the variation region rather than the beginning. The results in our model would seem to indicate that the Lg phase never substantially interacts with this incline. But if the bulge region were lengthened, but not so much as to allow for a double-bottoming, then a substantial interaction with this inclined floor might be possible, and corresponding attenuation of the post-variation region might be observed. Thus we hypothesize that the attenuation measured at the terminus of the model might be controllable and be a cyclical function of the length of the bulge region. This of course could be tested quite easily (but has not as yet beentested). 3.6 Effects on Pg compared to Lg The effects on Pg (Figs 8 and 9 right) are for the most part pretty similar to the Lg results (Figs 7 and 9 left). In the baseline traveltime curves, three general differences are notable, and are obviously unrelated to the pinch, bulge or scattering-region structures. One is that the Pn phase is notably excited in the Pg case, whereas the Sn is all but absent in the Lg case. This could indicate that the Moho transition in our model is more leaky to Pg at near-critical incidence angles, allowing a greater fraction of energy into the transition region and upper mantle (Shaw & Orcutt 1984), where the energy can propagate at the faster Pn velocity. The second difference to note is that the series of path multiples making up the Pg wave train seems truncated at a small number beyond approximately 400 km range, in contrast to the much longer wave train that develops in the Lg phase. This may be related to the aforementioned increased transmission into the Moho transition region, which may prevent the Pg phase from supporting the higher-order path multiples. Thirdly is that the Pg phase develops a low energy, long-tailed coda region that expands into a much wider velocity window than what is seen in the Lg case. This is a result of scattering combined with the higher intrinsic Q of the P propagation mode. Every model we simulated, even if absent any regions of high scattering, does have a very small amount of background heterogeneity (1 per cent) specified in each layer, and although the Pg will experience proportionally less scattering events than the Lg per distance travelled (MFPP = 13669 km, compared with MFPS = 6802 km in the baseline crust), those phonons that do scatter will sustain much longer owing to improved intrinsic Q of the P mode (⁠|$Q_{\mathrm{P}}^{\mathrm{(int.)}}=3374$|⁠, |$Q_{\mathrm{S}}^{\mathrm{(int.)}}=1500$| in the crust). The effect of that is seen in the mild, but long, coda in the Pg results. Additionally, the P to S scattering conversion channel is generally much stronger than the S to P conversion channel (Sato et al.2012), which explains the energy seen in the Lg velocity window, even absent a distinct Lg phase. Despite these overall differences in phase presentation, very little differs qualitatively or quantitatively between the Pg and Lg pinch, bulge and scatter results. There is, for the most part, a one-to-one correspondence of graphical features between the traveltime plot and energy curve results from both sets, except for two notable distinctions, to be discussed below. The similarity between Pg and Lg results is perhaps seen most clearly in the energy curve results (Fig. 9). Here, the energy trends for each model case can be seen to be quite similar between Pg and Lg in most respects. Perhaps the most distinct difference between Pg and Lg visible in the energy curve results is in the overall attenuation at the terminal range of the model (950 km), which can be seen in Fig. 9 and is also tabulated in Table 4. The attenuation of the Pg phase was less than the attenuation of the Lg phase in each model case. For example, in the pinch-with-overlay model, the Lg phase was attenuated by −8.11 dB at 950 km, whereas the Pg phase in the same model experienced an energy attenuation of only −4.73 dB. Similar reductions in the end range attenuation are seen in all other model cases as well. (Indeed, in the bulge model, the Pg phase actually experiences a very slight amplification at 950 km, compared to −1.54 dB attenuation for Lg.) A simple conclusion to be drawn from this is that the model structures studied here are significantly more efficient Lg blockers than they are Pg blockers. The second distinct difference to note is visible in the traveltime plots. In the Pg results we see a phenomenon that is absent in the Lg results: mode conversions. Particularly, we see very clear Pg to Lg conversions in the pinch-without-overlay, bulge and scattering-region results, originating from within the variable region. Again, a preferential conversion channel from P to S for scattering and a significant P to S reflection coefficient seem to be in play. Even the retrograde backscatter phase in Fig. 8(c), (which is also seen in the Lg case in Fig. 7(c), appearing as a weak reflection off the high-scattering crust region), matches an Lg-like phase velocity, despite being excited by an incident Pg wave. Inclusion of these converted phases in the non-temporally windowed Pg energy curves in Fig. 9 somewhat veils the Pg-blocking activity of the crustal structure. Nevertheless, since the energy is not lost but rather shifted in arrival time, we believe the whole-recording energy curves still provide useful insight, especially when coupled with the traveltimeplots. 4 COLLISIONAL MOUNTAIN BELT WITH UNDERTHRUSTING Surficial geologic interpretation and reflection profiling of collisional mountain belts often reveal a significantly greater degree of structural complexity associated with collisional mountain belts than that represented by our simple canonical models. Campillo et al. (1993), for example, demonstrate that complex over- and underthrusting of crustal and mantle wedges may explain zones of narrow, overlapping, corridors in efficient and inefficient Lg for paths crossing the Alps. To understand at least some of these effects, we designed the models shown in Fig. 10, consisting of an underthrust crustal sheet sampled by wavefields incident from either direction of the dipping underthrust crustal layer. Lg efficiency, inferred from the traveltime plots of Fig. 10 and energy–distance plots in Fig. 11, is similar to the simpler canonical model of crustal thickening (NSCP 03), with end range (950 km) attenuation relative to baseline of −1.94 dB for the case of Lg crossing the collision from the underthrust side and −1.71 dB for Lg crossing the collision from the overthrust side, as compared to −1.54 dB for the crust bulge. Directly above the anomalous zone, however, sharp differences in efficiency are seen between the wavefield sampling the structure from the under- or overthrust side, as well as a fundamental difference in both cases from that of the simple canonical model. For receivers directly above the anomaly with a wavefield arriving from the underthrust side, Lg is predicted to be less efficient than one arriving from the overthrust side. It should be noted in passing that this experiment with a strike-slip source is not an appropriate test of the elastic reciprocity theorem, in which the input source and receiver must instead be strictly given in terms of vector components of forces and observed motions (Gangi 1970). Together with the previously cited effects of the length of a zone of crustal thickening, the complexity of Lg amplification and attenuation for sources and receivers within a realistic mountain and the lack of reciprocity of source and receiver positions outside the mountain belt will degrade the reliability of Lg–Q tomographic inversion as a tool for predicting Lg propagation efficiency. Figure 10. Open in new tabDownload slide Traveltime curves and model profiles for a collisional mountain belt, for conditions where upthrusting slab drift is towards the source (left) and away from the source (right). Source event is an ideal strike-slip earthquake and azimuth is chosen to detect both Pg and Lg phases. Gross model structure is the same as the NSCP models with the exception of additional mesh points through the crust to enable modelling of distinct over- and underthrusting slabs. Figure 11. Open in new tabDownload slide Crustal thickness effects of collisional mountain belt on Lg energy by upthrust direction. ‘Upthrust left’ refers to case where upthrusting slab is migrating leftwards (towards source) and corresponds the left panel in Fig. 10, and ‘Upthrust right’ to the right panel. The baseline condition for the upthrust models is schematically identical to that of the NSCP models, although with a more intricate mesh structure. Curves represent unwindowed time-integrated all-component energies as a function of distance from source for an azimuth that selects the Lg phase as the initial wave train. End range attenuation relative to baseline for the Upthrust-left case is −1.94 dB and for the Upthrust-right case is −1.71 dB. 5 CONCLUSIONS Previous numerical and observational studies of high-frequency regional seismic propagation have found that phase attenuation and blockage can be due to a combination of factors, for example, Kennett (1986), Maupin (1989) and Sens-Schönfelder et al. (2009). These factors include intrinsic attenuation, scattering attenuation, crustal thickness variations, source depth and source mechanism. In this paper we have described the application of a fast, efficient, radiative transport algorithm to the modelling of Lg propagation that allows all of these factors to be included and tested against observations. Some general effects of crustal thickness variations have been demonstrated for several source depths and mechanisms. For earthquakes and explosions at shallow depths in the crustal waveguide having propagation paths crossing crustal thinning regions, Lg is amplified at receivers within the thinned region but strongly disrupted and attenuated at receivers beyond the thinned region. For the same source depths and mechanisms having propagation paths crossing regions of crustal thickening, Lg amplitude is attenuated at receivers within the thickened region, but experiences little or no reduction in amplitude at receivers beyond the thickened region. Localized regions of intense scattering within laterally homogeneous models of the crust increase Lg attenuation but do not disrupt its coda shape. The effects of crustal thickness variations on Pg correlate with those on Lg, but are much stronger on Lg. In future work, we plan on systematically testing the effects of a broad spectrum of the earth model and source parameters for regional phase paths crossing at least one specific region each of crustal thinning and thickening. These may include regions that have been studied by other investigators, including the North Sea, Barents Sea, Sea of Japan, Tibet Plateau, and the Pyrenees. SUPPORTING INFORMATION Supplementary data are available at GJI online. Earthquake_NSCP01.mp4.gz Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. ACKNOWLEDGEMENTS This research was supported by the U.S. Air Force Research Laboratory under contract FA9453-15-C-0069. We thank two anonymous reviewers for helpful comments and corrections. REFERENCES Campillo M. , 1987 . Lg wave propagation in a laterally varying crust and the distribution of the apparent quality factor in central France , J. geophys. Res. , 92 ( B12 ), 12 604 – 12 614 .. 10.1029/JB092iB12p12604 Google Scholar Crossref Search ADS WorldCat Campillo M. , Feignier B., Bouchon M., Béthoux N., 1993 . Attenuation of crustal waves across the Alpine Range , J. geophys. Res. , 98 ( B2 ), 1987 – 1996 .. 10.1029/92JB02357 Google Scholar Crossref Search ADS WorldCat Cao S. , Muirhead K.J., 1993 . Finite difference modelling of Lg blockage , Geophys. J. Int. , 115 ( 1 ), 85 – 96 .. 10.1111/j.1365-246X.1993.tb05590.x Google Scholar Crossref Search ADS WorldCat Chazalon A. , Campillo M., Gibson R., Carreno E., 1993 . Crustal wave propagation anomaly across the Pyrenean Range. Comparison between observations and numerical simulations , Geophys. J. Int. , 115 ( 3 ), 829 – 838 .. 10.1111/j.1365-246X.1993.tb01495.x Google Scholar Crossref Search ADS WorldCat Contrucci I. , Nercessian A., Béthoux N., Mauffret A., Pascal G., 2001 . A Ligurian (Western Mediterranean Sea) geophysical transect revisited , Geophys. J. Int. , 146 ( 1 ), 74 – 97 .. 10.1046/j.0956-540x.2001.01418.x Google Scholar Crossref Search ADS WorldCat Furumura T. , Kennett B.L.N., 2001 . Variations in regional phase propagation in the area around Japan , Bull. seism. Soc. Am. , 91 ( 4 ), 667 – 682 .. 10.1785/0120000270 Google Scholar Crossref Search ADS WorldCat Furumura T. , Hong T.-K., Kennett B.L.N., 2014 . Lg wave propagation in the area around Japan: observations and simulations , Prog. Earth Planet. Sci. , 1 ( 10 ), 1 – 20 .. 10.1186/2197-4284-1-10 Google Scholar OpenURL Placeholder Text WorldCat Crossref Gangi A.F. , 1970 . A derivation of the seismic representation theorem using seismic reciprocity , J. geophys. Res. , 75 ( 11 ), 2088 – 2095 .. 10.1029/JB075i011p02088 Google Scholar Crossref Search ADS WorldCat Grad M. , Tiira T., ESC Working Group, 2009 . The Moho depth map of the European Plate , Geophys. J. Int. , 176 ( 1 ), 279 – 292 .. 10.1111/j.1365-246X.2008.03919.x Google Scholar Crossref Search ADS WorldCat Gregersen S. , 1984 . Lg-wave propagation and crustal structure differences near Denmark and the North Sea , Geophys. J. Int. , 79 ( 1 ), 217 – 234 .. 10.1111/j.1365-246X.1984.tb02852.x Google Scholar Crossref Search ADS WorldCat Keers H. , Nolet G., Dahlen F.A., 1996 . Ray theoretical analysis of Lg , Bull. seism. Soc. Am. , 86 ( 3 ), 726 – 736 . Google Scholar OpenURL Placeholder Text WorldCat Kennett B. , Engdahl E., Buland R., 1995 . Constraints on seismic velocities in the Earth from travel times , Geophys. J. Int. , 122 ( 1 ), 108 – 124 .. 10.1111/j.1365-246X.1995.tb03540.x Google Scholar Crossref Search ADS WorldCat Kennett B.L.N. , 1986 . Lg waves and structural boundaries , Bull. seism. Soc. Am. , 76 ( 4 ), 1133 – 1141 . Google Scholar OpenURL Placeholder Text WorldCat Maupin V. , 1989 . Numerical modelling of Lg wave propagation across the North Sea Central Graben , Geophys. J. Int. , 99 ( 2 ), 273 – 283 .. 10.1111/j.1365-246X.1989.tb01687.x Google Scholar Crossref Search ADS WorldCat McNamara D.E. , Walter W.R., 2001 . Mapping crustal heterogeneity using Lg propagation efficiency throughout the Middle East, Mediterranean, Southern Europe and Northern Africa , Pure appl. Geophys. , 158 ( 7 ), 1165 – 1188 .. 10.1007/PL00001219 Google Scholar Crossref Search ADS WorldCat Mendi C.D. , Ruud B.O., Husebye E.S., 1997 . The North Sea Lg-blockage puzzle , Geophys. J. Int. , 130 ( 3 ), 669 – 680 .. 10.1111/j.1365-246X.1997.tb01861.x Google Scholar Crossref Search ADS WorldCat Rodgers A.J. , Ni J.F., Hearn T.M., 1997 . Propagation characteristics of short-period Sn and Lg in the Middle East , Bull. seism. Soc. Am. , 87 ( 2 ), 396 – 413 . Google Scholar OpenURL Placeholder Text WorldCat Sanborn C.J. , Cormier V.F., Fitzpatrick M., 2017 . Combined effects of deterministic and statistical structure on high-frequency regional seismograms , Geophys. J. Int. , 210 ( 2 ), 1143 – 1159 .. 10.1093/gji/ggx219 Google Scholar Crossref Search ADS WorldCat Sato H. , Fehler M.C., Maeda T., 2012 . Seismic Wave Propagation and Scattering in the Heterogeneous Earth , 2nd edn, Springer-Verlag . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Sens-Schönfelder C. , Margerin L., Campillo M., 2009 . Laterally heterogeneous scattering explains Lg blockage in the Pyrenees , J. geophys. Res. , 114 ( B7 ). Google Scholar OpenURL Placeholder Text WorldCat Shapiro N. , Béthoux N., Campillo M., Paul A., 1996 . Regional seismic phases across the Ligurian Sea: Lg blockage and oceanic propagation , Phys. Earth planet. Inter. , 93 ( 3–4 ), 257 – 268 .. 10.1016/0031-9201(95)03069-7 Google Scholar Crossref Search ADS WorldCat Shaw P. , Orcutt J., 1984 . Propagation of PL and implications for the structure of Tibet , J. geophys. Res. , 89 ( B5 ), 3135 – 3152 .. 10.1029/JB089iB05p03135 Google Scholar Crossref Search ADS WorldCat Vaccari F. , Gregersen S., 1998 . Physical description of Lg waves in inhomogeneous continental crust , Geophys. J. Int. , 135 ( 2 ), 711 – 720 .. 10.1046/j.1365-246X.1998.00681.x Google Scholar Crossref Search ADS WorldCat Wu X.-Y. , Wu R.-S., 2001 . Lg-wave simulation in heterogeneous crusts with surface topography using screen propagators , Geophys. J. Int. , 146 ( 3 ), 670 – 678 .. 10.1046/j.1365-246X.2001.00489.x Google Scholar Crossref Search ADS WorldCat Zhang T.-R. , Lay T., 1995 . Why the Lg phase does not traverse oceanic crust , Bull. seism. Soc. Am. , 85 ( 6 ), 1665 – 1678 . Google Scholar OpenURL Placeholder Text WorldCat © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society. TI - Modelling the blockage of Lg waves from three-dimensional variations in crustal structure JF - Geophysical Journal International DO - 10.1093/gji/ggy206 DA - 2018-08-01 UR - https://www.deepdyve.com/lp/oxford-university-press/modelling-the-blockage-of-lg-waves-from-three-dimensional-variations-0TQDc34e89 SP - 1426 EP - 1440 VL - 214 IS - 2 DP - DeepDyve ER -