Crustal deformation rates in Kashmir valley and adjoining regions from continuous GPS measurements from 2008 to 2019

We present GPS velocities in Kashmir valley and adjoining regions from continuous Global Positioning System (cGPS) network during 2008 to 2019. Results indicate total arc normal shortening rates of ~ 14 mm/year across this transect of Himalaya that is comparable to the rates of ~ 10 to 20 mm/year reported else-where in the 2500 km Himalaya Arc. For the first time in Himalayas, arc-parallel extension rate of ~ 7 mm/year was recorded in the Kashmir valley, pointing to oblique deformation. Inverse modeling of the contemporary deformation rates in Kashmir valley indicate oblique slip of ~ 16 mm/year along the decollement with locking depth of ~ 15 km and width of ~ 145 km. This result is consistent with the recorded micro-seismicity and low velocity layer at a depth of 12 to 16 km beneath the Kashmir valley obtained from collocated broadband seismic network. Geodetic strain rates are consistent with the dislocation model and micro-seismic activity, with high strain accumulation (~ 7e−08 maximum compression) to the north of Kashmir valley and south of Zanskar ranges. Assuming the stored energy was fully released during 1555 earthquake, high geodetic strain rate since then and observed micro-seismicity point to probable future large earthquakes of Mw ~ 7.7 in Kashmir seismic gap.

These studies are inconclusive due to limited data duration from sparse network of GPS sites in Kashmir. We present results from total of 16 cGPS sites data (Table 1) out of which 10 cGPS sites are spatially spread across the length and breadth of Kashmir valley, 4 cGPS located northeast of the valley in Zanskar ranges and 2 cGPS sites further northeast at Leh and Hanle during 2008 to 2019. Our study gives well constrained, precise estimates of crustal deformation and strain rate in Kashmir Himalayan transect with more data and dense network of stations.

Results and discussions
GPS displacements. Velocities of all the cGPS sites of the present study are estimated using the methodology detailed in Data and Methods section. ITRF14 and Indian plate frame of reference 17 velocities and the associated uncertainties are given in Table 1. India fixed velocities of cGPS sites in Kashmir valley and adjoining regions are plotted in Fig. 1. Arc normal and parallel deformation rates for sites located in Kashmir valley and the adjoining regions are computed and plotted in Fig. 3a,b. These rates are determined by rotating the site velocities, perpendicular and parallel to the local arc geometry as defined by 18 . In addition, we have plotted in Fig. 3a,b arc normal and arc parallel rates from the published 14,16 velocities of four GPS sites (ARJA, DHUL, URII, RAUJ) located to south of valley and five sites (RANK, KULG, NISH, SUMB, NARA) located in the valley (Figs. 1, 2, 4). For the rest of northwest Himalaya, we plotted arc-normal and arc-parallel rates from the published GPS velocities of continuous and campaign sites 15,19 . All the published GPS velocities are transformed to ITRF14 reference frame and converted to India fixed reference frame using the angular velocity of Indian plate given by 17 [1][2][3]7,8,14,42,43,47,48 Fig. 4). Post seismic after slip dislocation models of 2005 Muzaffarabad earthquake in IKSZ zone 7 suggest post seismic displacements in this region. Afterslip of this earthquake released 56 ± 19% of the seismic moment of main event. GPS measurements 10,11,16,25 Fig. 2) with error bars determined by rotating the site velocities on to the direction locally orthogonal to the arc using arc geometry defined by 18 . Figure was created using qtiplot software version 0.9.8.9 (https ://sourc eforg e.net/proje cts/qtipl ot.berli os/files /qtipl ot-0.9.8.9.tar.bz2/downl oad) 49 . (b) Arc-parallel rates of GPS sites in northwest Himalaya (boxed region of Fig. 2) with error bars determined by rotating the site velocities onto the direction locally parallel to the arc using arc geometry defined by 18 . Figure was created using qtiplot software version 0.9.8.9 (https ://sourc eforg e.net/proje cts/qtipl ot.berli os/files /qtipl ot-0.9.8.9.tar.bz2/downl oad) 49     www.nature.com/scientificreports/ There is no consensus regarding the geological structure beneath Kashmir valley from previous studies in this region. Deep Seismic sounding recordings in Kashmir Himalaya 26,27 indicate a complex Moho structure beneath the NW edge of the basin. For example 28 , reported Moho depth of 45 km beneath the Sopore (~ 37 km SSE of Kupwara; ~ 15 km NNE of BARM) and few km to the north Moho dips to a depth of 54 km near BAND cGPS site. Further 27 , reports that Moho dips to ~ 65 km beneath the Kanzalwan (North of Bandipora town), and further up-warping to the depth of ~ 60 km beneath the Nanga Parbat. The Moho offsets beneath the Sopore and Kanzalwan are ascribed to 2 deep rooted faults beneath the region, apparently oriented in NW-SE direction. The processed broadband data from our network (Supplementary section A1) could not find any evidence of the Moho dip between Sopore and Bandipora. However, increase in Moho depth from ~ 54 ± 2 km beneath the Sopore to ~ 64 ± 2 km north of Kanzalwan is clearly observed, indicating ~ 11° dip of the Moho along with a similar dip of the decollement underneath (Supplementary section A2). This step in decollement and Moho may represent one of the deep rooted faults described by 27  Three cGPS sites (PAND, BHIM, and KARG) located in the Indus Suture Zone (ISZ) (Fig. 1, Table 1) record relative motion of ~ 3 mm/year indicating active regional deformation in this region.  (Table 1) provide buried dislocation parameters and the associated slip. Assuming that the surface deformation is caused due to slip along the single planar dislocation (MHT) of finite length and width in elastic half space, the equation for slip, say m following 12 is given as m = (G T ·W e ·G) −1 G T ·W e ·d, Where G is a function of dislocation parameters, W e is weight matrix, d are GPS velocities. The slip m consists of two orthogonal components, dip-slip and strike-slip. The weight matrix is taken as the inverse of covariance matrix of the observed GPS velocities. Input to the inverse program are the range of a-priori dislocation parameters (starting point, dip, strike, depth, length, width), and GPS velocities. The output is the slip estimated using inversion technique for all possible combination of dislocation parameters. Best fit is the solution for which the residual between the observed and modelled displacements is minimum (Fig. 4). Complete detail of the dislocation modeling theory is given in 15 .
We used the GPS velocities from our study along with the published velocities 14,25 of 9 cGPS sites in this region (Fig. 4). The best-fit dislocation model for Kashmir valley gives a dip-slip of ~ 15.5 mm/year and strike-slip of ~ 4 mm/year along the single buried dislocation with a depth of ~ 15 km, length and width of ~ 145 km. Dip and strike angle of the buried dislocation is 7° and 142° respectively (Fig. 4). Micro-seismic events recorded by collocated broadband network in Kashmir valley and adjoining regions are plotted in Figs. 1 and 4 along with seismic events from ISC catalogue (Fig. 1) and one year local events from 5 . Micro-seismicity in Kashmir valley from collocated broadband data (Supplementary Section A1) points towards probable location of the northern edge of locked decollement plotted by dashed black line in the Fig. 1. In addition, results from Kashmir broadband network (Supplementary Section A1) found a low velocity layer at a depth of 12-16 km representing the locked decollement 6 . Our dislocation model is consistent with the northern edge of locked decollement, as highlighted by location of micro-seismicity, and the depth of decollement reported from Kashmir broadband network data.
Previous studies with limited spatial spread and duration of GPS data 16,25 and using 33  Modified least square method 35 is used in which adjustment based on the effect of nearest sites is made on the least square covariance matrix using a scale factor. This gives the choice to use different scales of analysis and evaluate the scale-dependent behavior of the observed system and select a scale factor that best fits the observed system. Grid-Strain program 34 enables to include or exclude a GPS point to estimate the effect of the specific point and is advantageous to compute strain rates using unevenly distributed cluster of GPS points. The inputs to the program are GPS velocities and positions and the output is the strain tensor. We used two dimensional grid strain program in which computations are performed at each node of XY grid (Fig. 5) and the output given as strain field tensors consisting of eigen vectors, eigen values and the normalized shear. Strain at a point is considered highly significant if the distance between the reference point and at least one GPS point in that sector is less than the scale factor chosen for the observed system. High significant strain for Kashmir valley along with associated rigid rotations are plotted in Fig. 5, which indicates high compression (red) with a mean rate of 4.5e−08 year −1 approximately normal to the arc in Kashmir valley and adjoining regions. Maximum compression of 7e−08 year −1 occurs along the northern edge of decollement. High extension (blue) rate of 7.4e−08 year −1 is observed to the northwest of the valley which is due to the anomalous westward motion recorded at KUPW GPS site. In Zanskar range and ISZ the strain rates are higher compared to the valley.   (3) to better constrain the crustal deformation in ISZ. Inverse modeling of GPS velocities in Kashmir valley and adjoining regions give oblique slip of ~ 16 mm/year along the 145 km wide decollement at a depth of 15 km. Our study using good spatial spread and long duration cGPS data gave reliable constraints on the depth of decollement which is consistent with the recorded seismic activity and low velocity structure 6 mapped around same depth (12-16 km). Strain rates estimated using the GPS velocities indicate high strain accumulation (7e−08, compression) to the north of the valley and south of Zanskar consistent with the northern edge of the locked decollement mapped using the seismic activity and inverse models. Strain rates (extension) to the SE of Kupwara are due to the anomalous westward GPS velocity of KUPW cGPS site which needs further investigation.
Historical earthquake record suggests about 14 damaging earthquakes occurred during 1123 to 1885, with only one earthquake of Mw ~ 7.6 in 1555 with poorly estimated rupture parameters. Kashmir seismic gap of ~ 250 km length between the rupture zones of 2005 Muzaffarabad and 1905 Kangra earthquakes did not record any major seismic slip during the intervening ~ 465 years. Considering that all the accumulated strain is fully released during the 1555 earthquake and fresh strain accumulation started since then, high strain rates recorded in this region give a scalar geodetic moment accumulation rate 36 of ~ 8.91e24 dyne-cm/year (~ 4.14e27 dyne-cm for 465 years) which is capable of generating Mw 7.7 earthquake. Accumulated slip deficit since 1555 event is ~ 7.4 m (16 mm/year for 465 years) and assuming rupture of the probable earthquake to be similar to 2005 Muzaffarabad earthquake, seismic moment 37 is ~ 5.3e27 dyne-cm suggesting Mw 7.8 earthquake is long due in this region. If the seismic moment released due to the after slip of the 1555 event is assumed to be similar to the 2005 earthquake 10 , the reduction in the Mw for a single event would be ~ 0.14. Micro-seismicity recorded in Kashmir valley further corroborate that a high magnitude earthquake is long due in this region. Our hypothesis is broadly consistent with the conclusions drawn by 6,14 that large magnitude earthquake is anticipated in this region, may be with low probability of occurrence.

Methodology
GPS data and analysis. Complete details of cGPS sites and data span used for the study is given in Table 1.
In addition to 16 cGPS sites in Kashmir we have used cGPS data of 9 sites in the rest of Himalaya and 10 plate interior sites (Fig. 1). All the cGPS sites are located on hard rock exposures (except the Kashmir University (UOKA) cGPS site) with 3-10 years of data. cGPS data spanning 12 years (2008-2019) is analysed using GAMIT/GLOBK software 38 along with 22 IGS sites (Table 1) after performing quality check on the data using TEQC software 39 . Loosely constrained daily solutions were obtained using sampling interval of 30 s and elevation cut off angle of 15° after eliminating data with short duration (< 20 h), multipath and several cycle slips. Errors related to satellite clock, receiver clock, phase ambiguities, atmosphere, phase center and multipath are minimized in the daily solutions as detailed in 19 . Final daily solutions are combined using GLORG to determine the positions and velocities of the cGPS sites in ITRF2014 reference frame 40 by stabilizing the IGS sites to their ITRF2014 positions and velocities. ITRF2014 velocities and the associated uncertainties are given in Table 1 for all the cGPS sites used in the present study. India fixed velocities (Table 1) and associated uncertainties of the cGPS sites are determined using the angular velocity of India-ITRF14 as given by 17  www.nature.com/scientificreports/