This electronic supplement contains a table of teleseismic events and figures of receiver functions (RFs), H-κ stacks, common conversion point (CCP) stacks, shear-wave splitting measurements, robustness estimates, and available focal mechanism solutions.
This electronic supplement contains selected results of the modeled Moho depths beneath the Kashmir basin using joint inversion of RFs and surface-wave dispersion, as well as Moho depths and Poisson’s ratio estimates for selected stations that were independently obtained using the H-κ stacking algorithm (Figs. S1–S12). The CCP RF image along a northeast–southwest profile (L2) across the valley (Fig. S13) is also included. A shear-wave splitting measurement at station BAR (Fig. S14) and tests regarding sensitivity of inverted models to different initial velocity models (Fig. S15) are also presented.
RF generation workflow, available focal mechanism solutions in Kashmir Himalaya, and adjacent areas followed by tests to constrain uncertainties in the estimates of Moho depths and tests for robustness of the low-velocity zone (LVZ) observed at a depth of 12–16 km are briefly discussed below.
RFs were generated using the time-domain deconvolution method of Ligorria and Ammon (1999). Deconvolution was continued for 200 iterations or terminated earlier if the misfit between the estimated and observed RFs was found to be <0.001. H-κ stacks were computed using the Eagar and Fouch (2012) package; the error estimation procedure used is discussed in Eaton et al. (2006).
The presence of an LVZ is qualitatively indicated by a prominent negative phase at ~2 s in most RFs (Fig. 3 of the main article; Figs. S11 and S12), and quantitatively better highlighted in their inverted models (see Fig. 7 of the main article). We tested the robustness of our inference concerning the existence of this layer by calculating synthetic RFs of a suite of rationalized inverted models for AHR, with and without the LVZ, and by comparing these to the corresponding ones calculated from data. We performed a similar test for GND and found the same effect which leads to the inference that the existence of this layer, perhaps at varying depths beneath the valley, is real (see Figs. S16–S18).
We examined 47 moment tensor solutions available for the study area, in the Global Centroid Moment Tensor (CMT) catalog, for the period 1976–2015 and 73 well-defined solutions for M >4.6 events during the same period in the International Seismological Centre (ISC, 2013) catalog. Both sets of solutions are almost identical for similar events (see Fig. S19).
Error estimates in a majority of the Moho depth determinations are ±2 km (say group 1, G1), with very few up to ±4 km (G2), as shown in Table 2 of the main article. Here, we present examples from both these groups, first from G2 (Fig. S20), for which inverse and forward modeling successfully highlighted the Moho, even though the Ps phase was not quite clear. In these cases, a perturbation of ±4 km around a median value in Moho depth was found to satisfactorily fit the RFs calculated from the data. For the G1 group of RFs, perturbations beyond ±2 km in Moho depths caused easily recognizable misfits (Fig. S21).
Triangulated maps of the Moho depths were obtained from forward modeling of inverted results and H-κ stacking (Fig. S22).
Table S1. List of the selected teleseismic events used in this study (blue dots in Fig. 2 of main article, headers are self-explanatory), with the back azimuth and great circle distances calculated at the station AHR.
Figure S1. (Left) Radial RFs (RRFs) and (right) tangential RFs (TRFs) for the station AHR. Note the gap in back azimuth (BAZ) from 160° to 280°. The polarity reversal of the Ps phase on the TRFs and its smaller amplitudes near the 100° BAZ, suggesting a possible role of anisotropy (Fontaine et al., 2013), are also notable.
Figure S2. H-κ stack for station AHR for events from BAZ 285° and 294°. These events pierce the Moho beneath the Pir Panjal range (see Fig. 1 of the main article); reported Moho depth is 57.5 ± 2.5 km.
Figure S3. Forward-modeling result for station AHR showing (a) stack of four RFs from BAZ ~289° and delta 39°–43°. The ±1σ error bounds are shown by black dashed lines. The RF corresponding to the jointly inverted model with p = 0.1 is shown by the blue line, and the synthetic RF corresponding to a rationalized version of the inverted model is shown by the red line. RFs corresponding to the latter with ±2 km perturbation in Moho depth are shown by red dashed lines. (b) Fundamental-mode Rayleigh group wave velocity from Mitra et al. (2006) is shown in a black line with ±1σ error bounds. The synthetic dispersion curve corresponding to the jointly inverted model is shown by the blue line and that of the rationalized model with perturbed Moho by the red line. (c) Shear-wave velocity profile beneath AHR up to a depth of 100 km. The blue line corresponds to the inverted model and the red line to that of the rationalized model. The initial model with a constant upper-mantle velocity of 4.48 km/s for initiating the inversion is shown by a black dashed line (for details on the initial model see the main article). The black arrow denotes the Moho depth.
Figure S4. H-κ stacks for the station BTP for events having BAZ between 99° and 101°. The Moho depth, VP/VS, and Poisson’s ratio are shown on top of the figure.
Figure S5. Forward-modeling results for a set of six RFs calculated for the station BTP, for BAZ of ~100°, and for delta between 82° and 85°. A sharp LVZ can be seen at a depth of ~12–16 km, for which thickness is larger than the resolving limit of ~2.6 km; this limit of resolution also places an uncertainty of about ±2.6 km on the inferred depth of the LVZ. Red and black arrows denote the depth of LVZ and Moho, respectively.
Figure S6. H-κ stacks for station GND for events having BAZ between 62° and 64°. The Moho depth, VP/VS, and Poisson’s ratio are shown on top of the figure.
Figure S7. Forward-modeling results for a set of RFs calculated for the station GND, details being the same as in Figures S3 and S5. A LVZ can be seen at a depth of ~10–13 km.
Figure S8. H-κ stacks for the station GUL for events having BAZ between 138° and 141°. The Moho depth, VP/VS, and Poisson’s ratio are shown on top of the figure.
Figure S9. H-κ stacks for the station NIL (an Incorporated Research Institutions for Seismology [IRIS] station in Pakistan). The Moho depth, VP/VS, and Poisson’s ratio are shown on top of the figure.
Figure S10. Forward-modeling results for a set of six RFs calculated for the station NIL, details being the same as in Figures S3.
Figure S11. The RRFs at station GUL in increasing order of back azimuth, which shows the distortion of the first P phase caused by Ps arrivals from a very shallow interface, most likely the base of the sedimentary layer in the valley, in the back azimuth range of 22°–90°. This is the only instance of the RFs used for which the signature of valley sediments is clearly identifiable. Careful location of our stations on extensive rock exposures most likely ensured this advantage.
Figure S12. The 0–5 s RFs for the stations AHR and GND; note the variation of arrival times and change in polarity of the negative phase at ~1–3 s. For details, see the main article.
Figure S13. The CCP RF image (Dueker and Sheehan, 1997; Zhu, 2000) along profile L2 (Fig. 1 in the main article). The distance along the x axis starts from 33.4°/74.4° and ends at 34.5°/75.6°. The unusual breadth of the positive amplitude at 50–60 km may be due to the LVZ multiples. A negative phase around 10–20 km represents the LVZ that can be traced only up to HAR/BTP starting from AHR. The data used are from June 2013 to April 2016. Interpretation of this image is tentative because LVZs beneath sedimentary basins have been reported to create errors in CCP stacks (Dueker and Sheehan, 1997).
Figure S14. SKS splitting measurement for an event recorded at station BAR using the SplitLab package (Wustefeld et al., 2008). SplitLab employs three methods for SKS measurements: the rotation correlation (RC) method, the minimum energy method, and the eigenvalue method. Here, we present the result obtained using the RC method that has tighter error bounds. A time delay of 0.5 s between fast- and slow-split components along 55° of azimuth is calculated using this method. Event description appears on top of the figure. (Top left) Transverse (blue line) and radial (red) components before the RC correction is applied. (a) Seismogram components in fast (red) and slow directions (blue line) after the RC normalization (i.e., correction of the delay between the fast and slow arrivals of the SKS phase). (b) Transverse (blue line) and radial (red) components after correction. (c) Particle motion before (blue line) and after (red) the correction is applied. (d) Map of the azimuth of fast axis versus delay time; gray region indicates the error bound on both values.
Figure S15. Tests of sensitivity of inverted models to the choice of the input-velocity model represented by a 100-km-thick homogeneous layer incorporated in the ak135 model of the continental lithosphere. The blue dashed line marks the value of the upper-mantle shear-wave velocity equal to 4.48 km/s up to a depth of 100 km. The red, blue, and green lines, respectively, show the inverted models in which the above value is perturbed by ±0.2 km/s. The closeness of the Moho depths and the velocity structures inverted with VS = 4.68, 4.48, and 4.28 km/s up to the first 40 km of the crust indicates that the inversion is not significantly affected by the choice of the initial model. We, therefore, chose the median value of 4.48 km/s for inverting all RFs, because this also happens to be the upper-mantle velocity in the seismic ak135 model of the continental lithosphere and has been used for inversion of crustal structures in other segments of the Himalaya (Rai et al., 2006; Acton et al., 2011).
Figure S16. Summary of three tests to investigate the existence of the LVZ. The bottom RF (blue line) is a forward-modeled RF from the inversion of a stack of RFs at BAZ 36° for station AHR; the corresponding blue line on the left shows the velocity model from which it is generated using the Herrmann (2013) package. The velocity model shown by the red line contains only the LVZ. The corresponding RF shows the Ps phase from this layer at ~2 s and the PpPs and PsPs multiples at ~5 and ~7 s. The velocity model in green contains only the Moho, and the corresponding RF shows the Ps phase followed by the Moho multiples PpPs and PsPs. RF of the model in black, which contains both the LVZ and the Moho, shows an enhanced Ps phase compared to the green one and broadened, obviously because of superposition of the PsPs multiple generated by the LVZ, on the Moho Ps phase. Red and black arrows denote LVZ and the Moho, respectively.
Figure S17. Summary of tests similar to those in Figure S16 for the station GND. Here, the PsPs multiple of the LVZ RF (red) arrives at ~5.5 s, much earlier than the Moho Ps phase at ~7.5 s, which is, therefore, not reinforced by the positive PsPs multiple of the LVZ, as in the case of AHR. It is therefore unbroadened and smaller in amplitude, as brought out by the Moho along Ps phase for AHR. Forward modeling with selected features of the inverted model thus enables one to clearly identify the Moho Ps phase that would ordinarily be indistinguishable from the LVZ-generated PsPs multiple. With the ambiguity in the identification of the Moho Ps phase at GND thus resolved, we tested the resulting inverted solution by forward modeling to constrain the Moho depth and its uncertainty bounds at GND to be ~60 km. This is shown in Figure S18.
Figure S18. Forward modeling for a set of RFs calculated for the station GND for BAZ of ~101° and delta between 81° and 85°, details being the same as in Figure S3 and S5.
Figure S19. The available focal mechanism solutions in the Kashmir basin and adjoining areas (ISC, 2013); the figure beneath each focal mechanism plot is in the format YYYYMMDD-MAG (year, month, day, and magnitude in Mw). T and P axes are shown by white and black circles, respectively. The scale (shown outside of main figure) shows the depth of earthquakes in kilometers. Most of these events have a thrust mechanism, as observed elsewhere in the Himalaya and that occurred at depths of ~10–25 km. Recent thrust events, such as the 2005 Mw 7.6 Muzaffarabad event, occurred at a depth of ~20 km and with a dip angle of 29° (U.S. Geological Survey [USGS]) on a shallower fragment of the Main Himalayan thrust (MHT; Avouac et al., 2006). The 2013 mb 5.7 Kishtwar event occurred at a depth of 16 ± 3 km with dip of 26° and might have occurred on Main Boundary thrust (MBT; Mitra et al., 2014).
Figure S20. A typical forward-modeling example in respect to a four-RF stack at the station BTP; details are the same as in Figure S3. (a) [Error bound ± 2 km (EB2)]: (top) Note that the RF Ps phase is not well resolved by inversion (blue line), because it incorporates a shift of Ps phase with respect to the observed RF, because, although the forward modeling (bottom right) positions the Moho at ~64 ± 2 km, a reasonable option to resolve such cases is to perturb the Moho depth by ±4 km (or more), so as to ensure that the Moho can be more definitively inferred to lie within the corresponding depth range (64 ± 4 km), albeit with a larger uncertainty. (b, EB4) The same example as in (a) with the Moho depth perturbed by ±4 km. Note that this depth range encompasses both the earlier estimates (inversion and forward model) given in EB2.
Figure S21. A set of forward-modeling results for the station AHR, details being the same as in Figure S3. (a) [Error bound ± 4 km (EB4)]: here, forward modeling is first carried out by assuming an uncertainty of ±4 km on the Moho. Note that in (bottom left) and (bottom right) the estimates look reasonable, but (top) suggests a possible tighter constraint to ensure a better fit of the Ps phase. (b, EB2) The same as in (a) with Moho perturbation, shown by a red line, reduced to ±2 km.
Figure S22. The triangulation of Moho depths beneath the Kashmir basin obtained independently by (1) using joint inversion and (2) H-κ stacking. (Left) Joint inversion: At each station, RFs from close-enough sources, 10° in BAZ and 5° in delta, are stacked and an average position of a representative Moho piercing point calculated using the TauP package. Stacking smoothens random noise and increases the signal-to-noise ratio (SNR) by , in which N is the number of RFs used for stack, and joint inversion of RFs and surface-wave dispersion data yields an optimized velocity model using the constraints of average crustal velocity imposed by the latter in defining the seismic discontinuities highlighted by the former. The uncertainty of an inverted solution for the velocity structure is then constrained by forward modeling. Thus, Moho depths were determined for most stations with an uncertainty of ±2 km (see Table 2 of the main article). Here, Delaunay triangulation is applied such that only points which are within the convex hull are interpolated. Note the crustal thinning along the basin and further deepening along the northeast direction. (Right) Moho depths are determined by forming H-κ stacks, in respect to selected RFs with clearly identifiable multiple phases. For error estimates on these data, see Table 2 of the main article. Here also, crustal thinning beneath the valley is observed, despite a greatly reduced data set (restricted to RFs with clear multiples). It should be noted that uncertainties in the results yielded by H-κ stacking over sedimentary basins and tectonically complex areas are subject to high errors (Dueker and Sheehan, 1997).
Acton, C. E., K. Priestley, S. Mitra, and V. K. Gaur (2011). Crustal structure of the Darjeeling–Sikkim Himalaya and southern Tibet, Geophys. J. Int. 184, 829–852.
Avouac, J.-P., F. Ayoub, S. Leprince, O. Konca, and Don V. Helmberger (2006). The 2005, Mw 7.6 Kashmir earthquake: Sub-pixel correlation of ASTER images and seismic waveforms analysis, Earth Planet. Sci. Lett. 249, 514–528.
Dueker, K. G., and A. F. Sheehan (1997). Mantle discontinuity structure from midpoint stacks of converted P to S waves across the Yellowstone hotspot track, Geophys. J. Int. 102, 8313–8327.
Eagar, K. C., and Fouch M. J. (2012). FuncLab: A MATLAB interactive toolbox for handling receiver function datasets, Seismol. Res. Lett. 83, 596–603, doi: 10.1785/gssrl.83.3.596.
Eaton, D. W., S. Dineva, and R. Mereu (2006). Crustal thickness and VP/VS variations in the Grenville orogen (Ontario, Canada) from analysis of teleseismic receiver functions, Tectonophysics 420, 223–238, doi: 10.1016/j.tecto.2006.01.023.
Fontaine, F. R., H. K. Tkalcic, and B. L. N. Kennett (2013). Crustal complexity in the Lachlan orogen revealed from teleseismic receiver functions, Aust. J. Earth Sci. 60, 413–430, doi: 10.1080/08120099.2013.787646.
Herrmann, R. B. (2013). Computer programs in seismology: An evolving tool for instruction and research, Seismol. Res. Lett. 84, 1081–1088, doi: 10.1785/0220110096.
International Seismological Centre (ISC) (2013). On-line Bulletin, International Seismological Centre, Thatcham, United Kingdom, available at http://www.isc.ac.uk.
Ligorria, J. P., and C. J. Ammon (1999). Iterative deconvolution and receiver-function estimation, Bull. Seismol. Soc. Am. 89, 1395–1400.
Mitra, S., K. Priestley, V. K. Gaur, S. S. Rai, and J. Haines (2006). Variation of Rayleigh wave group velocity dispersion and seismic heterogeneity of the Indian crust and uppermost mantle, Geophys. J. Int. 164, 88–98, doi: 10.1111/j.1365-246X.2005.02837.x.
Mitra, S., S. K. Wanchoo, and K. Priestly (2014). Source parameters of the 1 May 2013, mb 5.7 Kishtwar earthquake: Implications for seismic hazards, Bull. Seismol. Soc. Am. 104, 1013–1019, doi: 10.1785/0120130216.
Rai, S. S., K. Priestley, V. K. Gaur, S. Mitra, M. P. Singh, and M. Searle (2006). Configuration of the Indian Moho beneath the NW Himalaya and Ladakh, Geophys. Res. Lett. 33, L15308, doi: 10.1029/2006GL026076.
Wustefeld, A., G. Bokelmann, C. Zarolib, and G. Barruol (2008). SplitLab: A shear-wave splitting environment in Matlab, Comput. Geosci. 34, 515–528.
Zhu, L. (2000). Crustal structure across the San Andreas fault, southern California from teleseismic converted waves, Earth Planet. Sci. Lett. 179, 183–190.
[ Back ]