Mantle earthquakes in the Himalayan collision zone

The Himalaya-Tibet continental collision zone has featured prominently in a debate on whether strength in the lithosphere resides in two layers in the brittle crust and uppermost mantle separated by a weak lower crust, or in a single layer largely limited to the crust (Chen et al., 1981, 2012; Chen and Molnar, 1983; Zhu and Helmberger, 1996; Henry et al., 1997; Cattin and Avouac, 2000; Jackson, 2002; Chen and Yang, 2004; Beaumont et al., 2004; Monsalve et al., 2006; de la Torre et al., 2007; Bendick and Flesch, 2007; Liang et al., 2008; Priestley et al., 2008; Craig et al., 2012). The distribution of seismicity with depth in strained regions serves as a proxy for the lithospheric strength profile because the occurrence of earthquakes implies sufficient strength to sustain brittle failure, due to either cold conditions (e.g., Chen and Molnar, 1983; Sloan and Jackson, 2012; Blanchette et al., 2018) or local embrittlement processes (e.g., Incel et al., 2017, 2019; Prieto et al., 2017). The task is therefore to map where exactly local earthquakes occur with respect to layers within the crustal and uppermost mantle column. A difficulty is posed by uncertainty in estimated depths for both local seismicity and the Moho. Crustal earthquake depths are more difficult to constrain than their epicenters because the stations used to locate them are on the surface and there is a strong tradeoff between location depth and P and S velocity structure above the events.

Intralithospheric interfaces such as the Moho can be seen with teleseismic converted waves (receiver functions), but determining their depth suffers from the same tradeoff with velocity structure above the interface as in estimating earthquake depths. Comparing hypocentral depths with structural depths is hampered by the uncertainties in both parameters. The comparison is accurate only if the same velocity models are used for depth determination for both data sets. More typically, results are compared between different studies, and a lack of consistency between the velocity models used exacerbates the uncertainties. We propose a simple workaround to the depth determination uncertainty by comparing seismicity and structure directly as delay times between shear and compressional waves (“S minus P”, or S-P), rather than converting both to depth first.

We used data from the CE 2001–2003 HIMNT (Himalayan Nepal Tibet) seismic broadband experiment (Fig. 1). Of the seismicity located using the network (Monsalve et al., 2006), deep events are seen in Nepal just south of the Lesser Himalaya (Fig. 1, cluster C) and under Tibet (Fig. 1, clusters A and B). Figure 2B shows the same set of events on a previous structural depth profile from receiver functions on a line crossing the collision zone in a N18°E orientation perpendicular to the local range strike (Schulte-Pelkum et al., 2005). The Moho is visible as a contrast in isotropic velocity, and the Main Himalayan thrust as a contrast in anisotropy (Fig. 2B; Schulte-Pelkum et al., 2005). Even in this best-case scenario, with relocations and structural imaging done with the same stations and the same velocity models, the deep seismicity appears diffuse in depth and is difficult to clearly assign to crust or mantle.

Local earthquake depths are usually calculated from the travel-time difference between P and S arrivals from the event at each station. In receiver function studies, delays between teleseismic P and converted S are used to determine depths to interfaces under the station. In both cases, the S-P delay time is mapped to depth by assuming P and S velocity models between the event or converting interface and the station. Teleseismic arrivals and local events close to the station have similar steep incidence angles (Fig. 3C). Direct comparison of such S-P delay times between the two data sets therefore avoids the introduction of bias from using velocity models for depth migration and from unaccounted-for lateral variations, because the ray paths sample similar volumes. The two comparison data sets are constructed as follows.

For the local seismicity used in the S-P delay comparison, we used a previously relocated and published local event set of ∼500 earthquakes (Monsalve et al., 2006; Fig. 1), each located with at least 10 P and S arrivals picked with the HIMNT network; pick statistics are shown in Figure DR1 in the GSA Data Repository¹. For each event in this catalog located within the network footprint, we selected stations with P and S picks at epicentral distances of <35 km to sample a near-station volume. For these event-station pairs, we calculated the time difference between S and P arrival pick times (S-P) and its uncertainty. The S-P pick uncertainty was calculated by propagating the original uncertainty in the P and S arrival picks. Pick uncertainty and location error have long-tailed distributions, and we chose their half-amplitude widths of ∼0.2 s (S-P uncertainty) and 1 km (location error) as cutoff values (Fig. DR2), leaving 136 S-P differential arrival times from 98 events (black outlines in Fig. 1; blue events in Figs. 2B and 2C).

We used a previously published teleseismic Ps receiver function set from the HIMNT network (Schulte-Pelkum et al., 2005) to determine S-P times for the comparison to local seismicity. We additionally processed receiver functions to obtain a new higher-frequency set (Fig. 4) using automated event and receiver function quality-control criteria detailed by Schulte-Pelkum and Mahan (2014). Arrival times for the Moho conversion were picked from the original receiver function set for each station after corrections for slowness, with errors based on the width of the arrival and azimuthal variations seen at the station (Fig. 2C). A comparison with the second higher-frequency receiver function set shows close agreement in picked Moho times (examples in Fig. 4B). A group of stations in southern Tibet shows a positive amplitude arrival preceding the Moho peak (Fig. 4B), which marks the top of a layer with elevated lower-crustal velocities seen previously under this as well as other networks along the Himalaya (Yuan et al., 1997; Kind et al., 2002; Wittlinger et al., 2009; Schulte-Pelkum et al., 2005; Nábělek et al., 2009; Zhang et al., 2014). As in these studies, we picked the later arrival as the Moho at these stations because the resulting Moho is contiguous with that determined at neighboring stations. We also picked Main Himalayan thrust arrival times on the peak amplitude of the first azimuthal harmonic at each station (Schulte-Pelkum and Mahan, 2014).

(where Δt_z is the vertical time difference, Δt is the measured time difference, d is the epicentral distance, and Vp and Vs are the P and S wave speeds, respectively); so our S-P time comparison is not completely velocity free. However, Figure 3 demonstrates that the error introduced by a possible bias in velocity model is of second order because of the small epicentral distances involved, whereas a velocity model bias leads to a first-order error in depth determination.

After correcting all data to vertical incidence, we plotted the vertical S-P delay times from local seismicity and from teleseismic receiver functions on the same range-perpendicular profile (Fig. 2C). The profile is the same as the distance-depth profile in Figure 2B, but the vertical axis is now S-P delay time instead of depth. This representation introduces a distortion of the depth scale but minimizes the effect of the unknown velocity model (Fig. 3) and shows the relative depth relationship between the Moho and earthquake hypocenters. To test for any influence of lateral variations in Moho depth, we compare individual deep events to nearby receiver-function Moho piercing points in Figure 4.

Potential mantle seismicity in the original set of ∼500 events (Fig. 1, dots with black or white outlines) can be roughly grouped into three clusters (Fig. 1). Cluster B, just northwest of Kanchenjunga, falls outside our analysis set because of the lack of nearby stations. Deep events in cluster C near the CE 1988 near-Moho Udayapur earthquake (46–50 km estimated depths; summary in Ghimire and Kasahara, 2007) coincide with downtimes at HIMNT station GAIG or exceed the location or pick error criteria. Five events in cluster A (Fig. 1, dots with black outlines) pass our distance and uncertainty criteria, and all of them plot below the Moho in the S-P profile (Fig. 2C, between 0 and 50 km distance on the profile). These best-constrained events are not the deepest in the depth profile (Fig. 2B). An additional event just south of station RBSH plots within RBSH’s Moho error bar on the S-P profile (Fig. 2C, near –10 km profile distance). The magnitudes of these six deep events range from M_L 2.4 to 2.9 (Table DR1).

The larger error bar on the Moho time picked at station RBSH is due to azimuthal variations of the Moho pick in the receiver functions from that station. To exclude the influence of lateral variations in Moho depths or in velocities along ray paths, we compared deep-event S-P times to those from receiver functions with the closest piercing points. Figure 4A shows the map locations of the five sub-Moho and one near-Moho events in Figure 2C, as well as receiver function piercing points at 75 km depth for nearby stations. Figure 4B compares the S-P delay times of those events to receiver functions from nearby piercing points (larger dots in Fig. 4A). Receiver function stacks for those piercing points are shown for the same frequency band as used in the cross section in Figure 2B as well as for a higher-frequency band to obtain sharper Moho resolution. The shallower event south of RBSH (event 6) is below even shallower Moho picks for nearby piercing points. The considerable difference in Moho pick times between receiver functions from southern compared to northern back-azimuths at RBSH may represent actual Moho topography, but may also be due to lateral variations in velocity structure above the Moho.

The six deep events’ delay times are 0.3–1.1 s larger than the Moho delay times (Fig. 4B). Using uppermost mantle velocities from the standard AK135 model (Kennett et al., 1995) as well as those for Tibet from Monsalve et al. (2006), this delay time range corresponds to event depths of 3–11 km below the Moho. None of the deep events in southern Tibet that pass the error and distance criteria plot above the Moho. In an S-P profile with more events including those with larger location errors and pick uncertainties (Fig. DR3), only one event plots above the error bar of the Moho in southern Tibet despite more scatter overall.

Shallower seismicity is concentrated in the southern Tibetan upper crust, with one event just above the Main Himalayan thrust (Fig. 2C, profile distances 0–150 km). No events are seen in the Indian crust (below the Main Himalayan thrust and above the Moho) under southern Tibet, from ∼ –50 km profile distance to the north. This bimodal depth distribution in seismicity and the concentration of deep seismicity in cluster A persists from the 2001–2003 recording period of the HIMNT network (Monsalve et al., 2006; Huang et al., 2009) through the 2003–2005 period of the HI-CLIMB (Himalayan-Tibetan Continental Lithosphere during Mountain Building) network (Liang et al., 2008; Carpenter, 2010), which covered a footprint including that of HIMNT. Combining this catalog for the 2001–2005 time span with our placement of the deep seismicity below the Moho, we conclude that the seismicity in the Tibetan portion of the profile (north of the High Himalaya) shows a seismogenic upper (Asian) crust, an aseismic middle and lower (Indian) crust, and a seismogenic mantle.

In Nepal, seismicity is focused near the Main Himalayan thrust downdip from the locked portion (profile distances –120 to –20 km), and is distributed through all crustal levels south of the Lesser Himalaya (distances –180 to –120 km). Mantle seismicity in this area is seen in other studies (Monsalve et al., 2006; Huang et al., 2009) and is also suggested by our S-P profile when including events with larger location and pick errors (Fig. DR3). In contrast to the bimodal depth distribution north of the High Himalaya, the Indian plate therefore appears as a single seismogenic layer in this region before it begins its descent under the Himalaya.

Continental lower crust and mantle should normally deform in the ductile regime, and the presence of seismicity at these depths requires either cold material (<∼600 °C; Chen and Molnar, 1983; Sloan and Jackson, 2012) or mechanisms that allow brittle failure under higher-than-normal temperature and pressure conditions (Thielmann et al., 2015; Incel et al., 2017, 2019). Eclogitization of a metastable subducting Indian continental lower crust in the presence of water has been proposed as a source for the deep seismicity under Tibet (Lund et al., 2004; Jackson et al., 2004; Jamtveit et al., 2018). While the presence of higher-than-average velocities in the Indian lower crust under Tibet supports eclogitization (Schulte-Pelkum et al., 2005; Monsalve et al., 2006, 2008; Hetényi et al., 2007; Huang et al., 2009), our observation that the seismicity is in the mantle rather than the lower crust in this location contradicts the interpretation that eclogitization triggers deep crustal seismicity.

While seismicity in the mantle has been interpreted as requiring temperatures <600 °C (e.g., Jackson et al., 2008), there are embrittlement mechanisms in the mantle under ductile conditions that do not require cold temperatures, such as dehydration embrittlement, thermal runaway shear, and grain-size reduction (e.g., Kelemen and Hirth, 2007; Thielmann et al., 2015; Incel et al., 2017). Continental mantle earthquakes have been observed in regions of high temperatures (e.g., in southern California, USA; Inbal et al., 2016) and with rupture characteristics supporting such mechanisms (e.g., for the deep mantle seismicity under the Wyoming province, USA; Prieto et al., 2017). The laterally clustered nature of the deep seismicity in southern Tibet (Fig. 1; Monsalve et al., 2006, 2009; Huang et al., 2009; Carpenter, 2010) may point toward such mechanisms operating locally.

We acknowledge field efforts by R. Bilham and the HIMNT field teams, and discussions with P. Molnar, T. Stern, and K. Mahan. Comments by the editor, two anonymous reviewers, Hélène Lyon-Caen, Alexander Blanchette, and Simon Klemperer helped improve the manuscript significantly. This research was supported by U.S. National Science Foundation grants EAR-0310372, EAR-0538259, EAR-1246287, and EAR-1645009.