G Slip distribution beneath the Central and Western Himalaya inferred from GPS observations

of ≤ in this region, the elastic around the plate boundary. events can be to slip deﬁcits where the Indian and Eurasian Plates are locked during interseismic periods. Geodetic measurements can help discriminate the distribution of the interlocking areas and the steadily slipping areas beneath the Himalaya. To understand the deformation across the Central and Western Himalaya and the associated slip on thrust faults, campaign-mode GPS data were collected in the Garhwal–Kumaun region of the Western Himalaya. GPS sites velocities show that the deformation is currently concentrated between the Lesser and Higher Himalaya. Horizontal velocities are used to estimate the slip distribution. For the estimation, a model of the interseismic surface deformation caused by buried non-uniform creep dislocation (NUC) on a curved fault surface is used. The slip distribution shows that there might be structural discontinuity on the fault between the Kumaun and Garhwal regions of the Himalaya. The estimated slip rate at the depth around 20–40 km in the Central Himalaya and at the depth of ∼ 15 km in the Western Himalaya is 10 mm yr– 1 . The NUC model indicates that the shallow part ( < 20 km) of the thrust fault system along the plate boundary is almost locked. The locking depth appears to be deep in the Central Himalaya and shallower in the Western Himalaya. Further, most of the historical large seismic events are observed to have occurred in an area with a slip velocity less than 10 mm yr– 1 (i.e. in a locked zone) at the plate interface.


I N T RO D U C T I O N
The continuous process of continent-continent collision between the Indian and Eurasian Plates has led to the formation of the Himalayan range, since ∼50 Ma. This has contributed to the high elevation, fold-thrust formation and large crustal shortening along the ∼2500-km-long stretch of the Himalayan arc (Molnar & Tapponnier 1977;Srivastava & Mitra 1994). Three major south-verging thrust systems have developed along the length of the arc (Fig. 1). The northernmost fault is the Main Central Thrust (MCT), which emerges along the southern edge of the Higher Himalaya and separates the high-grade metamorphic rocks of the High Himalayan Crystalline Sequence from the weakly metamorphosed series of the Lesser Himalaya (Nakata 1989). The Main Boundary Thrust (MBT) marks the southern edge of the Lesser Himalaya and forms a series of north-dipping thrust faults, and it separates the lowgrade metasediments of the Lesser Himalaya from the sandstones of the sub-Himalaya. The southernmost fault is the Main Frontal Thrust (MFT), which is considered to be the most active of the three faults, and bounds the northern limit of the exposed Indian Plate.
The Indian Plate continues to move northeastward and to converge with the Eurasian Plate at the rate of ∼50 mm yr -1 (Patriat & Achache 1984;Bilham et al. 1997;Larson et al. 1999). Geological and seismological evidences suggest that the convergence rate in the Himalaya is ∼20 mm yr -1 (Molnar & Deng 1984;Armijo et al. 1986;Molnar & Lyon-Caen 1989;Bilham et al. 1997). Of the total convergence, ∼50 per cent occurs in the Himalaya, and the remaining involves distributed deformation (Demets et al. 1994;England & Molnar 1997;Holt et al. 2000;Wang et al. 2001;Chen et al. 2004a,b). The focal mechanism solution (Molnar & Lyon-Caen 1989) suggests that the convergence axis rotates counter-clockwise from west to east in the arc-perpendicular direction. This is consistent with the east-west hangingwall extension across southern Tibet (Armijo et al. 1986;Larson et al. 1999). Moderate to large earthquakes that have occurred along the Himalayan arc are shown in Fig. 2. The rupture zone between two great earthquakes, namely, the 1905 Kangra earthquake (M 7.8) and the 1934 Bihar-Nepal earthquake (M 8.1), is referred to as the seismic gap and is ∼800 km in length (Khattri 1987;Bilham 1995). This zone is considered to be a potential source of large earthquakes (Bilham et al. 2001). These earthquakes cause a slip on the MFT along the Himalayan arc and can be attributed to slip deficits on the locked plate boundary (Pandey & Molnar 1988;Ambraseys & Bilham 2000). It is important to understand the present-day deformation to asses the seismic hazard of this area. The Global Positioning System (GPS) technique is an effective and unique tool to precisely measure deformation along individual thrust faults; such measurements are valuable for understanding the convergence between the Indian and Eurasian Plates.

G P S DATA A N A LY S I S A N D R E S U LT S
A GPS network was established in the Kumaun region ( Fig. 3) of the Western Himalaya in 2005 May. At each site, a concrete pillar was constructed, and a non-magnetized steel plate was fixed over the pillar to settle a GPS antenna. GPS measurements commenced in 2005 December, and campaign measurements have been repeated every year since then. We tried to obtain GPS data during the winter season so as to avoid seasonal variations. We also collected GPS data from the Garhwal Himalaya (west of the Kumaun Himalaya) during 1999-2003, after the 1999 Chamoli earthquake (M 6.9). The locations of the GPS sites and the observation times are listed in Table 1. Trimble 4000SSE receivers were used for the observations; in addition, choke ring antennas with Dorne Margolin elements were used to reduce multipath effects on the signal. Also, a satellite elevation cut-off angle of 15 • was used to reduce multipath effects. All the stations were continuously operated for three to four days, and the sampling interval was 30 seconds. We analysed the GPS data together with data from some IGS fiducial sites by using the GAMIT/GLOBK (King & Bock 2005) post-processing software. Double-difference GPS phase observations were used to estimate the station coordinates, receiver clock parameters, atmospheric zenith delay and horizontal delay gradients and orbital and Earth rotation parameters. To reduce errors in the GPS-derived positions, the solid Earth tides and the ocean loading effect were taken into account by using the global tide model FES2004 (Lyard et al. 2006;Vergnolle et al. 2008).
To obtain the GPS site coordinates and velocities with reference to ITRF2000 (Altamimi et al. 2002), loosely constrained daily solutions were combined with the permanent tracking solutions of the IGS stations (Herring 2002). We determined the weights for each set of combined quasi-observations by adjusting the chi-square increments per degree of freedom through both forward and backward filtering of data to achieve the optimum value for the repeatability of the coordinates (McClusky et al. 2000).
Estimated ITRF2000 velocities of GPS stations are converted to velocities in the fixed Indian frame by using the Euler pole (52.970 • N, -0.297 • E with angular velocity of 0.499 • Myr −1 ) of relative motion defined by Banerjee et al. (2008). Velocity vectors with respect to the Indian Plate are shown in Fig. 3. Times-series of estimated site coordinates at some of the stations in the GPS network are shown in Fig. 4.
From the velocities with respect to the Indian frame, it can be seen that sites near the MCT show more variations than the stations near the MBT and MFT. This may suggest that the updip of the plate interface is locked beneath the study area. The zone with the maximum convergence rate, which accommodates much of the shortening, is well north of the frontal thrust system, and the convergence rate is comparable to the geological slip rate over the last few million years (Lav'e & Avouac 2000); further, the MFT is locked (Bilham et al. 1997;Larson et al. 1999;Jouanne et al. 1999Jouanne et al. , 2004. It is also observed that DHAR is moving southward while MUNS is moving northwestward. A cluster of microseismicity ( Fig. 3) can be seen to the southeast and northwest of MUNS. The activity of microseismicity is relatively low and non-homogeneously distributed in the region (Tsapanos 1990;Wason et al. 2002;Khan & Chakaraborthy 2007). The region MUNS situated may be high strength homogeneity.

I N V E R S I O N A P P ROA C H
Many authors (Vergne et al. 2001;Jouanne et al. 2004;Feldl & Bilham 2006;Bettinelli et al. 2006) have discussed the slip rate beneath the Himalaya by using creeping dislocation rather than the back slip model. These authors have assumed that the convergence rate is proportional to the creeping dislocation below the locked portion of the MFT. In the present study, the method devised by Yabuki & Matsu'ura (1992), which adopt Akaike's Bayesian Information Criterion (ABIC) to select a weighting coefficient for smoothness of slip distribution, is used to construct a model of the slip distribution on a curved fault surface. First, we try to understand the deformation in the Kumaun-Garhwal section, which lies close to the centre of the Himalayan fold and thrust belt. Then, we investigate the slip distribution beneath the Central and Western Himalaya by including GPS velocities for stations in Nepal (Central Himalaya) and the northwest Himalaya (Western Himalaya) obtained from the literature (Jade et al. 2004;Bettinelli et al. 2006;Banerjee et al. 2008). Velocity vectors with respect to the Indian frame are shown in Fig. 5.

Fault model
The fine subsurface structure was studied by deep crustal reflection seismic profiles of the region beneath the Greater and Tibetan Himalayan Zones in the INDEPTH project (Zhao et al. 1993;Nelson et al. 1996;Hauck et al. 1998). The fault geometry used in our analysis was obtained by considering their cross-sections of different observation lines presented in earlier studies (Jouanne et al. 2004;Yin 2006;Bettinelli et al. 2006). The assumed fault plane extends for 1700 km along strike (NE-SW direction) and 900 km along dip direction, and is meshed with 26 × 15 bicubic B-spline functions (Yabuki & Matsu'ura 1992). The plate interface was defined down to the depth of 120 km. We modelled the interseismic surface deformation in an elastic, homogeneous and isotropic half-space (Okada 1985) with a curved fault surface by using the buried non-uniform creep distribution (NUC) model (Fig. 6). We used ABIC approach, described by Yabuki & Matsu'ura (1992), to select the optimum values of the hyperparameter controlling the data fitness and the smoothness of the fault slip. For the best estimation of likelihood prior constraints by adjusting the hyperparameters were to get the minimum ABIC values, for achieving smooth distribution of fault slip model. Once the minimum ABIC values were obtained, the model parameters and covariance for the estimation errors were computed.

Resolution test
A resolution test was performed using slip patches of 30 mm (grey) and 0 mm (white) given on the fault plane shown in Fig. 7. The model is like Zebra test performed in the North Island subduction zone, New Zealand by Wallace et al. (2004). Slip patches were considered to be absent in the updip area of the fault model, since the plate interface of the MFT is locked. The random noise of standard deviation of 1.2 mm for horizontal component is added to the synthetic velocities at each site for the test. The result of the resolution test is shown in Fig. 7. The synthetic site velocities (Fig. 8) were reproduced fairly well. The locations of the areas of high slip matches the input model for the given slip distribution in the downdip portion, but it shows slightly less magnitude in the updip portion of the model. This shows that the observation network used in this study is capable of providing adequate resolution to constrain the slip distribution beneath the Himalaya.

728
M. Ponraj et al.  Past studies (Ni & Barazangi 1984;Srivastava & Mitra 1994) involving well log and earthquake data have indicated that this region is underlain by a basal detachment fault dipping gently to the north at an angle of ∼2 • and a depth of ∼20 km. The profile (b-b ) shows a lower slip rate in the plot of slip versus depth, probably because of less data in that section. The model reproducibility is shown in Fig. 9(b). Fig. 10(a) shows the optimum slip distribution obtained using the non-negative least-squares algorithms developed by Lawson & Hanson (1974); we assume a thrust-type dislocation on the grids. The model reproduces the overall site velocity pattern (Fig. 10b). Comparison of the model with the fault geometry (Fig. 10a) indicates a slip rate of 20 mm yr -1 in the fault depth range of 40-60 km. In the Central Himalaya and Western Himalaya, a creep motion of less than 10 mm is observed up to a depth of 30 km and 15 km, respectively. We infer that the Indian Plate is creeping in the deeper part of the frontal thrust fault, and this inference is consistent with that of earlier studies based on the dislocation model (Jouanne et al. 1999;Larson et al. 1999;Jouanne et al. 2004;Bettinelli et al. 2006). Bilham et al. (1997) have estimated the slipping rate of the Indian Plate to be ∼20 mm yr -1 at a depth of 20 km across the eastern Nepal Himalaya and at a depth of ∼25 km beneath western Nepal. The present model shows a low creeping rate (∼10 mm) at a depth of ∼30 km in the Nepal Himalaya. The MFT slip rate is estimated to be ∼14 mm yr -1 below a depth of 15 km in the Western Himalaya (Banerjee & Bürgmann 2002). It is also observed that the slip distribution in both the Central and Western Himalaya shows smooth variations, and a large slip rate (∼30 mm) is seen in the Central Himalaya at depths of 60-80 km.

D I S C U S S I O N
The velocity field obtained in this study indicates that the collision zone in the study area is not deforming rapidly. For the observation period, the GPS stations distributed around the MFT and MBT do not indicate any deformation, while those in the region around the MCT suggest the occurrence of deformation. Similar deformation has previously been observed along the Nepal Himalaya by using models of elastic strain accumulation for the region north of the locked Himalayan thrust system (Bilham et al. 1997;Larson et al. 1999). Focal mechanisms of earthquakes in the study region that occurred during the observation period are shown in Fig. 2; these earthquakes occurred to the north of the Himalaya and have been attributed to the thrust plane and the low-dip-angle strike plane (Bollinger et al. 2004;Jouanne et al. 2004). Slip vectors of thrust earthquakes that have occurred along the Higher Himalaya over the past few decades suggest that the earthquakes might have resulted from interseismic strains, because these events ruptured along shallow-dipping thrust faults (Pandey et al. 1995).
Profiles of the slip rate as a function of depth along 1-1 to 5-5 (Fig. 10c) are shown in Fig. 10(a) overall profile shows a smooth distribution of the slip rate with respect to the depth, except for profile 1, which lies across eastern Nepal. There may be possible trade-off between the magnitude of slip and depth in the lower crust of the Central Himalaya, since high value of the slip magnitude obtained in the model. It is also clear that the Central Himalaya (profile 2-2 ) shows differences in the structural anomaly in comparison to eastern and western Nepal. This suggests that the central Nepal segment (longitude: 82 • -84 • E) marks a discontinuity between eastern and western Nepal (Larson et al. 1999;Jouanne et al. 2004).

Comparison with seismicity
The Himalayan deformation is responsible for large earthquakes, and these earthquakes have ruptured the southern Himalayan arc (Seeber & Armbuster 1981;Pandey & Molnar 1988;Bilham et al. 2001), as shown in Fig. 2. A major portion of the current deformation is seen to lie between the Lesser and Higher Himalaya of the Kumaun-Garhwal region, which is coincides with the seismicity of the region (Fig. 3). This suggests that MCT is more active than MBT and basal detachment is present (Seeber et al. 1979;Ni & Barazangi 1984). A section of the plate boundary has not been ruptured over the last 300 yr (Bilham et al. 1997), and it is thought to be the location of great earthquakes. It is therefore necessary to evaluate the creep motion of the Indian Plate along the plate boundary. Recent earthquakes (Harvard CMT catalogue 2005-2007 along the plate boundary provide clues to the thrust and strike-slip mechanisms (Fig. 2). The seismicity in this region results mainly from the underthrusting of plates (Molnar et al. 1973;Ni & Barazangi 1984).
In the Nepal Himalaya, horizontal shortening due to ductile creep motion is observed along deep portions of the MFT, which is locked up to a depth of 20 km (Bettinelli et al. 2006). Along the strike of the Western Himalaya, the MFT is locked up to around 15 km (Banerjee & Bürgmann 2002). In our inversion model, we observe that the crust is creeping into the deep part of the MFT along the Nepal Himalaya and at shallower depths along the Western Himalaya. It is observed that the plate is creeping below 20 km and releasing elastic strain, which can cause minor to moderate earthquakes in the study region. Microseismic activities have been observed all along the Nepal Himalaya (Pandey et al. 1999), and it is being triggered by stress accumulation in the deep portion of the creeping zone (Cattin & Avouac 2000). Using our inversion technique, we have found that most of the historical large earthquake events have occurred in areas with slip velocities less than 10 mm yr -1 . Most of the earthquake events in this region are concentrated along the thrust zones and are confined to depths of about 20 km (Ni & Barazangi 1984).

C O N C L U S I O N
We have attempted to explain the observed crustal deformation in the Himalaya in terms of coupling between the Indian and Eurasian Plates. GPS observations near the MCT indicate the presence of deformation, which is significant when considering the seismic activity in the thrust areas of the MCT. The regional crustal deformation is examined by analysing campaign-mode GPS data recorded in this study, apart from published data for the Nepal Himalaya and the northwest Himalaya. It is possible that there is a discontinuity in the structural fault between the Garhwal and Kumaun regions of Western Himalaya. We propose a model for the interplate coupling in the updip portion of the plate boundary. The model (NUC model) indicates that the shallow part (<20 km) of the thrust fault system along the plate boundary is almost locked. The locking depth appears to be large in the Central Himalaya and shallower in the Western Himalaya. The inferred location of the locked area matches with the locations of large earthquakes in the Himalaya. We thank anonymous reviewers for constructive comments in improving the manuscripts. M. Ponraj (MP) conveys his sincere thanks to the Research Center for Predictions of Earthquakes and Volcanic Eruptions, Tohoku University, for providing him with excellent lab facilities. He is also grateful to the Japan Society for the Promotion of Science (JSPS), Japan, for providing financial support for his visits to Tohoku University, Sendai, Japan, and extends his sincere thanks to the Department of Science and Technology, New Delhi, for considering and recommending his application to the JSPS. We thank Dr. R. King for making the GAMIT/GLOBK GPS data analysis software available. We thank Prof. Mita Rajaram, Director, Indian Institute of Geomagnetism, Navi Mumbai, for encouragement and support. We also thank P. Wessel and W. H. F. Smith for the GMT software. We thank P. S. Sunil, S. K. Prajapati and K. Deenadayalan for their timely help.