Little Geodetic Evidence for Localized Indian Subduction in the Pamir‐Hindu Kush of Central Asia

Geodetically derived velocities from Central Asia show that Northern Afghanistan, the Tajik Pamir, and northwestern Pakistan all move northward with comparable large velocities toward Eurasia. Steep velocity gradients, hence high strain rates, occur only across the Main Pamir Fault zone and with lesser magnitude between the northernmost Hindu Kush and the south and southeast margins of the Tajik Depression. Localized shortening is not apparent on any active India‐Hindu Kush crustal boundary; hence, crustal convergence between India and Eurasia in Central Asia is absorbed primarily on the northern and western margins of the Pamir. This concentrated strain on the Pamir margins is consistent with one, geometrically complex, interface between subducting Asian lithosphere and the Pamir. That interface might curve westward such that the Hindu Kush seismic zone is a continuation of the Pamir seismic zone, or alternatively, Hindu Kush earthquakes might occur in convectively unstable mantle lithosphere mechanically detached from surface faults.


Introduction
The Pamir-Hindu Kush region of Central Asia (Figure 1) serves as the best present-day example of ongoing subduction of continental lithosphere. It is interpreted as a case of initiation of subduction in continental materials (Burtman & Molnar, 1993;Hamburger et al., 1992;Jay et al., 2017;Kufner et al., 2016;Negredo et al., 2007;Pegler & Das, 1998;Schneider et al., 2013;, in contrast to the more common scenario where continental lithosphere follows oceanic lithosphere into a subduction zone. Because subduction is generally attributed to gravitational foundering of negatively buoyant intact lithosphere into the asthenosphere (Isacks et al., 1968), while continental lithosphere is usually thought to be less dense than the mantle regardless of age, such subduction is unusual (McKenzie, 1969(McKenzie, , 1977Turcotte et al., 1977).
Beneath the Pamir and the Hindu Kush, a seismogenic zone reaching to approximately 350-km depth  has several subduction features, including crustal thrust faults (Burtman & Molnar, 1993), localized seismicity (e.g., Pavlis & Hamburger, 1991;Prieto et al., 2012;, and a zone of high seismic wave speeds and low seismic attenuation extending to mantle depths (e.g., Khalturin et al., 1977;Mellors et al., 1995;Schneider et al., 2013;Sippl, Schurr, Tympel, et al., 2013). The distribution of earthquake hypocenters below the crust is conventionally divided into two parts, a northeastern part beneath the Pamir and a southwestern part beneath the Hindu Kush (Figure 1), distinguished by the Pamir zone dipping shallowly to the south and the Hindu Kush zone dipping nearly vertically to the north (Billington et al., 1977;Chatelain et al., 1980;Roecker, 1982;Roecker et al., 1980;, with the two separated by a gap in hypocenters at depth. Fan et al. (1994) (2016), and Liao et al. (2017), among others, attribute both the reversal in dip direction and the gap in seismicity at depth to different sources of lithosphere for each part of the seismic zone, with the Pamir zone of Eurasian origin and the Hindu Kush zone of Indian origin. This "two-sided" scenario corresponds to the two-boundary models discussed below (model family 2). Because of the absence of a gap in the distribution of seismic hypocenters shallower than 80-100 km , as well as the narrowness of the deeper gap below~100-km depth, others interpret Pamir-Hindu Kush   (Pegler & Das, 1998), or Asian origin . The latter "one-sided" subduction corresponds to the one-boundary models discussed below (model family 1). An alternative model for the one-boundary models allows a downgoing, intact lithospheric slab hosting Pamir seismicity while Hindu Kush seismicity occurs in a negatively buoyant blob of mantle lithosphere, analogous to the interpretation that Houseman and Gemmer (2007) and Lorinczi and Houseman (2009) gave for the Vrancea zone beneath the Carpathians.
In contrast to the Pamir, although the high seismic velocity anomaly beneath the Hindu Kush ( Figure 1) is spatially coincident with hypocenters (Koulakov & Sobolev, 2006;Kufner et al., 2016;Mohan & Rai, 1995;Negredo et al., 2007), it does not project to a unique thrust fault at the surface (Sippl, Schurr, Tympel, et al., 2013). Moreover, strain rates inferred from seismic moments at intermediate depth are much higher than average horizontal convergence rates at the surface Zhan & Kanamori, 2016).

Data and Methods
We compile data collected from 2008 to 2016 from 66 Geodetic Positioning System (GPS) stations throughout Tajikistan, Afghanistan, Kyrgyzstan, and Pakistan, including 32 publicly available and 9 restricted campaign sites and 25 regional continuous sites including International GNSS Service (IGS) reference stations (Table S1 in the supporting information). We process these data using the GAMIT/GLOBK software package (Herring et al., 2015) following the procedure described in Reilinger et al. (2006). GAMIT is used to calculate initial daily position estimates of each station. These are edited, averaged, and weighted over approximately 2-week long intervals. GLOBK's Kalman filter is used to estimate linear horizontal velocities from the position averages, incorporating a random walk noise model to account for systematic errors. The velocity solution is tied to the International Terrestrial Reference Frame 2008 (ITRF08) and then transformed into a stable Eurasian reference frame using the ITRF08-Eurasia angular velocity calculated by Altamimi et al. (2012 ; Table S1).
Following the determination of the regional velocities, we define five families of models (0, 1, 2, 3, and 4) based on the number of tectonic boundaries in the study area ( Figures S1-S4). We use TDEFNODE (McCaffrey, 2009) to calculate angular velocities for rigid crustal domains that minimize misfit to the observed velocity field. We use the misfit calculated in TDEFNODE as a measure of the likelihood function of possible tectonic boundaries on the north margin of the Pamir; the north, east, and south margins of the Tajik Depression; within the central Pamir; and south of the Hindu Kush ( Figure 1). Because we only use the TDEFNODE results for misfit related to the presence or absence of particular tectonic boundaries, the resulting block models do not necessarily describe the full regional kinematics.
We set the specific location of boundaries using major faults in the Central Asia Fault Database (Mohadjer et al., 2016). Lacking detailed 3-D structural information about these faults, we approximate them as 20°dipping fault planes for thrust faults, 45°dipping planes for normal faults, and near-vertical fault planes for strike-slip faults. The TDEFNODE inversions are not sensitive to the structural geometry (such as 3-D fault shape, dip angle, or locking depth) because very few of the observed geodetic velocities are located within the elastic length scale of the block bounding faults. The zero-boundary case (model 0) corresponds to a rigid rotation of the whole study area, the one-boundary case (model 1a or 1b, Figure S1) to convergence between Eurasia and India on the Main Pamir Fault and its westward extension along the Darvaz-Karakul' Fault, the two-boundary case (model 2b, Figure S2) to convergence between Eurasia and the Pamir at the Main Pamir Fault and between the Pamir-Hindu Kush and India on a structure south of the Hindu Kush, and the three-and four-boundary approximations to additional active boundaries on the north side of the Tajik Depression and within the central Pamir, respectively (Figure 1). We present alternative configurations for each number of boundaries in the supporting information (Table S2 and Figures S1-S4). Boundaries defining the far-field edges of the model domain are specified for geometric simplicity rather than geologic or kinematic accuracy (Figure 1 and Figures S1-S4).
We then compare model favorability using the Akaike Information Criterion (AIC) and a modified version for limited observational data (AICc; Akaike, 1974). These tools are designed for model selection in nonunique problems by penalizing both large misfits to data and large numbers of free parameters. We specifically use the least squares case of the AIC and the AICc, , respectively (Burnham & Anderson, 2002), where n is the number of GPS velocity observations, b σ 2 is the rootmean-square misfit determined from the TDEFNODE inversion, and K is the number of free parameters. K increases by 4 with each additional boundary in the kinematic model, because the existence of the boundary and three parameters defining the angular velocity of the resulting additional block are added.

Kinematic Results
The GPS velocities (Figures 2 and 3 and Table S1) show ongoing convergence between India and Eurasia, with the steepest velocity gradients localized near the northern margin of the Pamir, consistent with other GPS results (Ischuk et al., 2013;Reigber et al., 2001;Zhou et al., 2016;Zubovich et al., 2010). The shortening rate on and adjacent to the Main Pamir Fault system is 18-22 mm/year (65-80% of the total India-Eurasia relative rate of 28 ± 4 mm/year, measured from Karachi; Mohadjer et al., 2010), increasing along strike eastward (Figures 2 and 3). Convergence across the Alai Valley, the topographic margin of the Pamir, accommodates~15 mm/year of this total N-S convergence across less than 50 km centered on or near the Main Pamir Fault (e.g., Burtman & Molnar, 1993;Ischuk et al., 2013;Liao et al., 2017;Zubovich et al., 2016) and Vakhsh Thrust. Between the Western Pamir and the Tajik Depression, there is approximately 15 mm/year of NW-SE convergence (Figure 3). Sites in the western Tajik Depression move with negligible velocities relative to Eurasia (Figures 2 and 3).
The southern margin of the Pamir-Hindu Kush-Karakorum from Peshawar to Chitral accommodates not more than 10 mm/year of India-Asia convergence over >250 km (Figures 2 and 3). Little convergence is allowed south of Peshawar, such as on the Salt Range Front Fault, since Karachi (KCHI in Table S1) moves north toward stable Eurasia at about 28 mm/year, while Peshawar (NCEG in Table S1) moves north at about 26 mm/year (Ischuk et al., 2013;Mohadjer et al., 2010). The Central and Northern Pamir also move only slightly less rapidly than Peshawar, at 18-24 mm/year relative to Eurasia  (Ischuk et al., 2013), precluding rapid north-south shortening in the Wakhan Corridor or southernmost Tajikistan (Figures 2 and 3).

Model Selection Results
The ΔAICc (Tables 1 and S2) compares the favorability of plausible tested models. The smallest AICc score corresponds to the model with the most empirical support; small ΔAICc values indicate models of similar favorability. Larger ΔAICc values indicate either that total misfit is large, that there are many underconstrained free parameters, or both.
Unsurprisingly, based on the AICc scores, a rigid rotation model with no boundaries is the least favorable, with the maximum ΔAICc of the configurations tested. Adding a curving boundary at the location of the steepest velocity gradients (green line in Figure 1) minimizes the AIC and AICc (Tables 1 and S2), indicating the most favored (in the AIC sense) model among those that describe the deformation in terms of boundaries. An alternative boundary on the north side of the Tajik Depression gives a similar AICc (ΔAICc = 1.7; Table S2), indicating that the observations support either variant, and we cannot distinguish between a boundary on the north or south side of the Tajik Depression with the available observations. However, a model with boundaries on both sides of the Tajik Depression is less favorable (Figures S1-S4 and Table  S2). Adding an additional boundary on the south side of the Hindu Kush increases the AIC and AICc by 5.3 and 8.9, respectively, for a model that is considerably less favorable but still with some empirical support (Tables 1 and S2). Adding additional boundaries through the Central Pamir produces unfavorable models, with large ΔAICc. Therefore, the presence of Eurasian lithosphere underthrusting the Main Pamir Fault footwall is well constrained by the velocity observations, but a slab of Indian lithosphere attached to the surface and underthrusting the Hindu Kush either south of or within the zone of anomalous seismicity is less so.

Discussion
Three explanations have been proposed for the Pamir-Hindu Kush intermediate depth earthquake zone: a contorted subducting slab of Eurasian material (one-sided subduction; Figure 4a), distinct slabs of Eurasian and Indian origin subducting beneath the Pamir and Hindu Kush, respectively (two-sided subduction; Figure 4b; e.g., Burtman & Molnar, 1993;Kufner et al., 2017;Liao et al., 2017), and foundering lithosphere (one-sided subduction plus convective overturn; Fillerup et al., 2010;Houseman & Gemmer, 2007;Lorinczi & Houseman, 2009). The first of these options is consistent with available constraints from prior geodesy (Ischuk et al., 2013;Mohadjer et al., 2010;Zhou et al., 2016;Zubovich et al., 2010) and the additional results presented here, plus slip on known faults (Arrowsmith & Strecker, 1999;Bernard et al., 2000;Coutand et al., 2002;Cowgill, 2010;Kuchai & Trifonov, 1977;Nikonov et al., 1983;Sobel et al., 2011;Strecker et al., 1995;Trifonov, 1978Trifonov, , 1983) and seismic imaging (e.g., Schneider et al., 2013; Tympel, et al., 2013), except for the inferred rapid sinking of the deepest part of the Hindu Kush anomaly from moment summation Zhan & Kanamori, 2016). Such a single slab model implies relatively simple dynamics: subduction needs to initiate and persist only once in continental lithosphere, presumably in an area of previous crustal thinning (Leith, 1982   Indian material must be mostly or entirely detached from the surface, such that there is no longer an active, high strain rate thrust system separating downgoing Indian lithosphere from overriding Hindu Kush crust. A third option, convective instability arising from lithospheric thickening, has been invoked in the Carpathians (Houseman & Gemmer, 2007;Lorinczi & Houseman, 2009;Fillerup et al., 2010) and also matches the observed lack of a surface boundary as well as the spatial distribution of intermediate depth seismicity and inferred high stretching rates Zhan & Kanamori, 2016). In this case, the origin of the foundering material could be Indian, Eurasian, or Tethyan lithosphere.
Although subduction is generally thought to initiate in old oceanic lithosphere (e.g., Carlson et al., 1983;McKenzie, 1977;Turcotte et al., 1977;Stein & Stein, 1996), which is typically more dense than asthenosphere, subduction of Eurasian continental lithosphere at the Pamir is now supported by many different observations and demonstrates at least the possibility of subduction initiation within continental lithosphere (Burtman & Molnar, 1993;Chatelain et al., 1980;Jay et al., 2017). Combined with the possibility of a Hindu Kush example of convective overturn of continental lithosphere like that invoked in the Carpathians (Houseman & Gemmer, 2007;Fillerup et al., 2010;Lorinczi & Houseman, 2009) or detachment of an Indian lithospheric slab, we infer greater exchange of mass between the continental lithosphere and the asthenosphere than has previously been assumed likely.

Conclusion
Most of the convergence between India and Eurasia is accommodated on the northern side of the Pamir and between the Fayzabad-area (NE Afghanistan) and the Tajik Depression. Based on an AIC, the most favorable regional tectonic boundary model of the options considered consists of a single tectonic boundary extending from the Main Pamir Fault to thrusts wrapping around the eastern and southern margins of the Tajik Depression then linking to the Chaman Fault system through Afghanistan (Lawrence et al., 1992;Szeliga et al., 2012). Additional tectonic boundaries introduce more free parameters to models without fitting the surface velocities much better. The paucity of localized shortening south of the Hindu Kush therefore favors one-sided subduction at the Pamir, at least in the present day, rather than two-sided subduction invoked in several recent dynamic models (e.g., Kufner et al., 2017;Liao et al., 2017;Schurr et al., 2014). Onesided models allow Hindu Kush seismicity to be hosted in the overturned Asian slab, a detached Indian slab, or convectively foundering lithosphere. Shortening in the Pamir-Hindu Kush-Karakorum is much slower and much more diffuse than at the Main Pamir Fault, precluding a present-day India-Hindu Kush subduction boundary in the crust.