Supporting information for: Imaging the Indian Subcontinentbeneath the HimalayaVera Schulte-Pelkum, Gaspar Monsalve, Anne Sheehan, Som Sapkota, M. R. Pandey, RogerBilham, & Francis WuReceiver function calculation and data selectionAll recorded teleseismic events with mb≥5.0 within P , P P , and P KP distance ranges basedon the NEIC catalogue were checked automatically for their suitability for receiver functionanalysis. For all records where the direct P arrival showed a minimum signal-to-noise ratio of3 on the vertical component, receiver functions were calculated using a time-domain iterativemethod [1] and those with a variance reduction of at least 70% were retained for analysis. Theautomated selection eliminated the majority of the initial events (including most of the events inthe P P , and all but a handful in the P KP distance ranges) and left only high-quality receiverfunctions and a few outliers per station, which were removed at visual inspection. Typically,transverse component receiver functions are noisier, and fewer are retained based on the vari-ance reduction criterion during the automatic processing than the number of radial receiverfunctions at the same station. Three southern Nepal stations situated on thick sediments (BIRA,JANA, GAIG) that displayed receiver functions strongly dominated by sediment multiples wereexcluded from the common point conversion stack. Most stations in Tibet yielded sparser datathan those in Nepal due to limited accessibility (less frequent servicing).Anisotropic layerWe discuss the azimuthal sections in Figure 3 in the main paper in further detail. The timingof the arrivals is slowness corrected to represent vertical incidence. In Figure 3a, note thephase with polarity reversal at 1.6 s on the radial component. The transverse component showsarrivals with polarity reversals at 0.6 and 1.6 s. Azimuthal variation of receiver functions can becaused by several different subsurface features, such as dipping isotropic interfaces, scattering,anisotropy, and other lateral variations. A polarity reversal on the radial component requiresstrong anisotropy; the periodicity of 360 degrees indicates a steeply plunging symmetry axis.While an arrival at 0.6 s on the radial component is mostly obscured by direct P , paired arrivalswith reversed polarity are clearly visible on the transverse component and originate from the topand bottom of a layer of only a few km thickness. Figure 3b shows synthetics [2] that match theobserved pattern at stations SIND in Figure 3a. The synthetic receiver functions are calculatedfor the average observed slowness and a moveout correction is applied in the same way as tothe data. Parameters for the synthetic model are as follows (for the symmetry axis, in this casea slow axis as defined by the negative P and S anisotropy, strike is clockwise from North, andplunge is positive down from strike azimuth):1Supplementary Table 1: model with anisotropic and sediment layers used to produce syntheticsin Figure 3bthickness ρ vpvs% vp% vsstrike plunge ηkm g/cm3km/s km/s deg deg2 2.5 4.0 1.6 0 06 2.9 5.4 2.8 -20 -20 18 -50 0.4∞ 3.2 6.5 3.7 0 0The forward modelling is nonunique, as tradeoffs exist between layer thicknesses and veloc-ities, and between anisotropy magnitude and velocity contrasts. Required features to model thearrivals between 0 and 2 s are strong anisotropy (≥ 15%), slow symmetry (i.e., planes of highvelocity perpendicular to a slow symmetry axis), alignment of the azimuth of the symmetryaxis with the NNE direction normal to the Himalayan arc, a steep dip of the fast planes downtowards NNE, and an anisotropic shape parameter η (which describes the distortion of the phasevelocity surface from an ellipsoid) at the lower end of the range typical for anisotropic crustalmaterials (0.4-0.9) [3]. While the polarity pattern in the radial and tangential components isgenerated by the thin anisotropic layer in an otherwise homogeneous halfspace, we have addeda sedimentary layer above it to improve the match in the waveforms. The sedimentary layer isalso required by geology (station SIND is located on the Siwaliks sediments) and for matchingSIND arrivals to nearby stations. Arrivals after 2 s display additional azimuthal complexity notexplained by the influence of the shallow anisotropic layer alone, and they are not modeled inthe synthetics.Figure 3c clearly shows the conversion from the bottom of the anisotropic layer at stationBUNG, situated 70 km north of station SIND. The topside conversion at BUNG is less clearthan at SIND, where there is a large contrast to overlying sediments (BUNG has fast shallowvelocities rather than a surface sedimentary layer, see section below). Similar azimuthal patternsare seen at all stations in Nepal. The areal extent of station locations with similar arrivalsfurther supports the presence of a continuous anisotropic layer as the cause, rather than lateralheterogeneities. Our data are insufficient to confirm or exclude continuation of the layer intoTibet, mostly due to reduced data coverage, but also because the arrivals are later and moreprone to modification by shallower structure. However, Figure 2c suggests that the decollementshear zone may persist under Tibet, and the hint of its presence we see corresponds well to thedecollement imaged by INDEPTH reflection data (Figure 2d).High-velocity anomaly at BUNGThe azimuthal variations are observed at all stations and render a standard analysis of moveoutwith azimuthal averaging nearly meaningless. Only a few stations have a sufficient number ofreceiver functions with a range of slownesses in a narrow azimuthal range to allow a moveoutanalysis. We show an example for station BUNG in Figure S1. Moho and midcrustal arrivalsappear early at BUNG compared to surrounding stations. Matching the midcrustal multiple2requires a high-velocity anomaly above 20 km depth, which brings the single midcrustal andMoho conversion in line with that observed at surrounding stations.Common conversion point stack and velocity model for depth migrationSIND, located on the Siwaliks sediments, and BUNG, on top of the high-velocity anomaly,were the only stations we applied shallow velocity corrections to in the common conversionpoint stack in Figure 2. Depth migration for all other stations were performed with separate1-D models for Nepal and Tibet, obtained from tomographic inversion of local and regionalevents (Monsalve et al., manuscript in preparation).