Variations of the crustal thickness in Nepal Himalayas based on tomographic inversion of regional earthquake data

Introduction Conclusions References


Introduction
Collision processes are related to the convergence of continental blocks and lead to significant shortening and thickening of the crust.The collision zones with strong seismic activity often coincide with highly populated areas, leading to damage and destruction of human habitation and suffering of population.The Himalayas, which are the highest mountain chain on the Earth, has formed due to the collision of the Indian and Asian plates.The mechanisms of mountain building in Himalayas and Tibet have been extensively discussed by many authors for decades (e.g., Dewey and Bird, 1970;Seeber et al., 1981;Molnar and Tapponier, 1975;Allegre et al. 1984).As the Indian landmass moved northwards the sedimentary piles with their older crystalline foundation complexly folded, faulted and thrusted, which caused varied crustal structure all along the 2500 km long Himalayan arc from west to east.According to the most popular tectonic model of Himalayan collision (Seeber et al., 1981), the Indian plate underthrusts the Asian plate along a gentle north-dipping (4-10 • N) detachment plane, called the Main Himalayan Thrust (MHT) (Fig. 1).Most of the Himalayan earthquakes are shallow and occur at 15 to 20 km depth on the MHT.Several recent seismological studies, however, suggest that the tectonic model varies from west to east.For example, earthquakes in the eastern Himalaya tend to be much deeper than in the western part (Kayal, 2001(Kayal, , 2010;;Mukhopadhyay and Sharma, 2010).More definitive geodynamic concepts can only be constructed based on reliable information on the deep structures in the crust and the mantle.However, due to many political and natural reasons the Himalayas is a difficult region to make a detailed study with most of geophysical methods.
The Mohorovičić (Moho) discontinuity depth is one of the key types of information which is directly linked with the major geodynamical processes.For the Himalayas, the existing Moho depth models are either too generalized or too fragmentary.For example, an existing global model CRUST2.0(Bassin et al., 2000) provides an over-smoothed solution in Himalayas with extrapolation in some parts.The gravity modeling in the Himalayas also provides fairly smooth variations of the Moho depth (e.g., Tenzer and Chen, 2014).Another gravity study by Jin et al. (1996) reported that the Moho depth varies from 38 km below Indo-Gangetic Alluvial Plains (IGAP) to about 75 km below high Himalayas.The regional tomography models though depict reliable images of the lithospheric behavior beneath Himalayas and Tibet (e.g., Li et al., 2008;Koulakov, 2011), but they cannot provide many constraints on the crustal structures.On local scale, the existing receiver function sites, deep sounding profiles and local earthquake tomography results from east to west of the Himalayas and Tibet (e.g., Kind et al., 2002;Kumar et al., 2005;Galve et al., 2002;Hauck et al., 1998;Mitra et al., 2005;Ramesh et al., 2005;Schulte-Pelkum et al., 2005;Rai et al., 2006, Mukhopadhyay andSharma, 2010) provide reliable but too local and sparse information which is hard to use to build a generalized crustal model for the entire Himalayas.
In this paper we make an attempt to estimate the variations of the crustal thickness based on results of tomography inversion using the travel time data recorded by the networks of Nepal and northern India.In most cases, seismic tomography is used to derive smooth velocity distributions and appears to be not sensitive to sharp first-order interfaces.However, in some cases it can provide useful information to estimate the variations of the main interfaces.For example, Koulakov and Sobolev (2006) provided the map of the Moho depth beneath the Middle East area based on the inversion of the regional travel time data from the ISC catalogue.This model is fairly corroborated by later studies based on receiver functions and active seismic profiles (Mechie et al., 2013).Koulakov and Sobolev (2006), however, put forward some conditions which would make possible studying the Moho depth using travel time tomography: (1) stations in the study area should be distributed densely and uniformly as much as possible; (2) sufficient number of sources should be located inside the study region; (3) size of the area should be in the range of 150-500 km; (4) both travel times of crustal (Pg, Sg) and mantle (Pn, Sn) rays should be presented in the data set.To some extent, all these conditions are fulfilled in the Nepal Himalayas region.Thus, we claim that the tomographic results in this study provide new information on the variations of the Moho depth beneath the Nepal Himalayas.

Data analysis and tomography model
We have combined the data of regional networks in northern India (run by India Meteorological Department, IMD) and Nepal (run by the Department of Mines and Geology, Nepal, DMN) along with the global ISC catalogue for the years of 2004-2014.In total, we used the information from 78 seismic stations installed in India and Nepal.The data selection was based on three criteria: (1) the residuals for the P data and S data after location of sources in the 1-D model should not exceed 2 and 3 s, respectively; (2) the number of picks per event should not be less than 8; (3) the distance from an event to the nearest recording station should not be more than 250 km.In total, 10 864 P arrival times and 5293 S arrival times from 821 events in the study region were selected for this study (on average, almost 20 picks per event).The distributions of stations and selected events used for computations are shown in Fig. 1.Note that only in Nepal we have fairly dense distributions of both stations and earthquakes.In China to the north, there were many events, but no stations were available; in India, there were some stations available, but very few events were reported.The analysis of data is performed using the iterative tomographic algorithm LOTOS (Koulakov, 2009).Because of the large size of the area, we have modified the code by taking into account the sphericity of the Earth.All the calculations are performed in the Cartesian coordinates.However, the reference model is kept radially symmetric, and Z coordinates for the events and stations are corrected according to the spherical shape of the Earth.In other aspects, the workflow of the analysis was similar to that used in other studies based on this algorithm (e.g., Koulakov et al., 2010).The processing starts with preliminary source locations with the use of reference table containing travel times in the 1-D model.In the next step, the sources are re-located using 3-D algorithm of ray tracing based on bending method.The velocity distributions are parameterized with nodes distributed inside the study area according to the ray density.To avoid any bias of the model due to predefined parameters of the grid, we performed the inversions for four different grids with different basic orientations.Examples of node distributions for two grids with basic orientations of 0 • and 45 • are shown in Fig. 2. Note that, in map view, the node projections look regularly spaced.However, along the vertical lines, the number of nodes and spacing depend on the data distribu- tion (cases with denser node distributions are depicted with darker points in Fig. 2).The inversion was performed simultaneously for the 3-D P velocity distributions and S velocity distributions, source parameters and station correction.The matrix was inverted using the LSQR method (Paige and Saunders, 1982;Nolet, 1987).The inversion results obtained using differently oriented grids are averaged into one model which was then used to update the 3-D model for the next iteration.In total, for the analysis of both synthetic and observed data, we used three iterations, each of which included the steps of source locations in the updated velocity models, matrix calculations and inversions.
To avoid any predefinition for the Moho depth, we set the reference model without any sharp interfaces and even without high gradient levels.We defined a constant V p /V s ratio equal to 1.75 and set the P-velocity values at different depth levels: 5 km s−1 at −1 km, 6 km s−1 at 25 km, 7.2 km s−1 at 40 km, 7.7 km s−1 at 65 km, 8 km s−1 at 120 km and 8.2 km s−1 at 210 km depth.Between these levels, the velocity was linearly interpolated.Starting velocities for P and S models are shown in Fig. 3 with dotted lines.
Unlike the tomography algorithms used by Koulakov and Sobolev (2006) for studying the Moho depth in the Middle East, here we do not parameterize the Moho as a sharp first-order interface with variable depth.Instead, we derive the geometry of Moho based on consideration of velocity anomalies and absolute velocities.In the starting 1-D model the velocity around the Moho depths was faster than expected crustal velocities and slower than mantle velocities.As a result, the crust was revealed as low-velocity anomaly, whereas the uppermost mantle is associated with high-velocity anomaly.The variation of thickness of the crust-related low-velocity anomaly may represent the perturbations of the Moho depth.After inversion, absolute velocity forms a zone of higher gradient around the presumed location of the Moho interface, as seen in representation of average absolute velocities in Fig. 3.Note that in this case we present an average for a long profile, for which the crustal thickness may vary significantly.For local points, the Mohorelated high-gradient zone is seen more prominently.
To examine the adequacy of the detection of crustal thickness, we have performed a series of synthetic tests.In Figs.4a and 5 we present result of one of the tests.The synthetic model was defined as a superposition of a reference 1-D velocity model and a low-velocity anomaly with the amplitude of −15 % of variable thickness representing the crust.The lower limit of this anomaly is indicated in vertical sections in Fig. 5 with solid lines.In map view (Fig. 4a), a zone of variable thicker crust is highlighted with a dotted line.Outside this zone, the Moho is set at 30 km depth.Note that starting 1-D model for the synthetic reconstruction was different for the "true" reference models.This represents the realistic situation in the case of observed data analysis when the true reference model is unknown.
To compute the synthetic data, we have used the same source-receiver pairs as in the case of the real data analysis.The computed synthetic travel times were perturbed with random noise having an average deviation of 0.1 s, which enables approximately the same variance reduction as in the case of observed data inversion.After computing synthetic travel times using the 3-D ray tracer, we "forgot" all information on the velocity distributions and source locations.Then we performed the full data processing including the steps of source locations.The restored anomalies at 35 km depth are shown in Fig. 4a.In vertical sections in Fig. 5, we present the restored relative anomalies and absolute velocities.It is seen that, in section 1 along Himalayas, the thickness of the derived low P-velocity anomalies correctly represents the undulations of the Moho interface in the input model, especially in the eastern part of the profile.We see that the location of the "true" Moho better corresponds to the velocity of 7.4 km s−1 (yellow zone) in the restored absolute velocity model.
For section 2 to section 5, the low-velocity anomaly is visible only beneath Nepal.Neither in the Indian nor in the Tibetan side is the crust-related anomaly restored.This test shows that the robust reconstruction of the Moho depth using the tomographic reconstruction can be achieved only in case of coexistence of stations and events in a sufficiently large area.Just availability of only stations (like in India) or only seismicity (in Tibet) is not enough for this purpose.The resolved area can be estimated from this test that roughly coincides with the shape of the restored crust-related anomaly.
The horizontal resolution is examined with another synthetic test, which is presented in Fig. S1 of the Supplement.
As mentioned earlier, for the inversion of the observed data, we used three iterations.During the inversion procedure, the residuals reduced from 0.807 to 0.500 s for the P data (37.99 % of reduction) and from 1.57 s to 0.83 for the S data (47.14 % of variance reduction).Note that the similar remnant residuals were obtained in the final synthetic model discussed above.
The results of tomographic inversions for P velocity are presented in five vertical sections: one along and four sections across Himalayas.We present both relative deviations with respect to the starting model (Fig. 6) and absolute velocities (Fig. 7).In addition, one horizontal section of Pvelocity anomalies at 35 km depth is presented in Fig. 4b.More horizontal sections, as well as S-velocity anomalies, are presented in the supplementary material in Fig. S2.We have defined the resolved area according to the results of various synthetic tests, mainly based on the capacity to retrieve crust-related anomalies in the test shown in Fig. 4. Areas with lower resolution in sections 2 to 5 are shaded in the resulting plots in Figs. 6 and 7.
As we see from the results of the synthetic test, thickness of the low-velocity anomaly beneath Nepal may represent the depth variations of Moho.We can also identify Moho in a zone of generally higher gradient of the absolute velocity, which is observed in plots of absolute velocities in Fig. 7. Based on the results of synthetic modeling, we have identified the Moho depth approximately at the contour line of 7.4 km s−1 (yellow zones in Fig. 7)., 6, 207-216, 2015  We have manually traced the lower limit of this lowvelocity anomaly in 21 vertical sections passing across the Himalayan chain (see Fig. S3) and created the 2-D surface of this limit beneath Nepal (Fig. 8a).Projections of this surface to the vertical sections below Nepal, where a satisfactory resolution is achieved, are shown in Figs. 6 and 7 with solid lines.It should be noted that the unambiguous tracing of the Moho is not possible everywhere.For example in section 5 to 7 (see Fig. S3) below the main low-velocity anomaly, there is another low-velocity pattern which appears to be weaker and separated from the upper one.This transitional anomaly may or may not be included to the crust.In case it is included, the total thickness of the crust in the frontal zone beneath Himalayas may reach 80 km (see Figs. S3 and S4), which seems to be less plausible.Thus the latter model may be accepted as more realistic model; the former one is shown in Supplement.Thick transitional zone in this area may represent underthrusting of one continental block underneath another one, which results in doubling of the crustal thickness.However, as shown by synthetic tests, the resolution in this part of the area is not high; thus, we should be prudent and avoid interpretations that are too speculative.

Solid Earth
In sections 4 and 5, we can compare our results with Moho depth determinations obtained in other studies.In section 4, dotted line depicts the Moho depth derived from combined interpretation of receiver function, gravity and thermomechanical modeling by Hetényi et al. (2006).Same as in our results, they depict flat Moho in the area where our model is resolved.To the north, their model shows gradual thickening of the crust; however, this trend is observed outside the resolved area of our model.In section 5, we can compare with the Moho depth determined from receiver functions by Schulte-Pelkum et al. (2005).This model shows generally same dipping trend as in our result.In both profiles, Moho interface in our results appears to be deeper to 5-7 km than the one determined by other authors.We should emphasize that the absolute values of Moho depth derived from tomography should be considered with prudence because of unambiguity of the conversion of continuous seismic anomalies into the interface.The relative variations of the crustal thickness, however, appear to be correct.

Discussion
The variations of the crustal thickness in the frontal zone of the Himalayan thrust belt, as seen in our tomographic model (Fig. 8a), may be attributed to variable mechanical properties of the collided plates.The existence of weaker or more rigid segments in the underlying Indian plate may cause weaker or stronger folding in the Himalayan thrust zone.However, due to several reasons, it is not easy to quantify this correlation because the Indian plate is mostly covered by thick sediments of the Gangetic alluvium, which hide major tectonic features.
To identify hidden crustal structures, the observations of potential fields might be useful.In Fig. 8b we show the freeair gravity anomalies for the Nepal Himalayas region ex- tracted from the global model by Andersen et al. (2010); the smoothed anomaly is obtained using a pseudo-Gaussian weight function with the characteristic radius of 10 km.In the Nepal Himalayas, the gravity field demonstrates very strong variations.To the south of the Himalayas, there is a strong negative anomaly, which is partly caused by isostatic compensation related to the mountain growth, and it might also reflect thick sediments of the Gangetic alluvium brought from the Himalayas due to very fast erosion.The maximum value of the free-air gravity field is observed in the higher Himalayas along the Nepal-China border.It is clear that these strong variations across the Himalayan thrust zone are mostly due to abrupt Moho dipping from relatively thinner crust in the Indian Plate to almost doubled crust beneath the Himalayas and Tibet.Along the Himalayas we also observe strong variations of the gravity anomalies which might be associated with laterally inhomogeneous thickness of the crust.The lateral variations in gravity anomalies in Nepal correlate rather well with our estimates of the crustal thickness.For example, areas of thinner crust indicated with "1", "3", "6" and "8" correspond to lower-gravity anomaly patterns.On the contrary, thicker crust segments numbered with "2", "4", "5", "7" and "9" are associated with higher values of gravity anomalies.
Here we also examine the magnetic anomalies extracted from the global compilation by Maus et al. (2009).Besides the map for the Nepal Himalayas and adjacent areas in Fig. 8c, we present the map of magnetic anomalies for a much larger area in Fig. S5.In the Nepal Himalayas, the correlation of crustal thickness with the magnetic anomalies is not as clear as that found with the gravity map (partly due to non-availability of high-resolution magnetic data in Nepal).However, it is worth noting that the largest pattern of thinner crust "6" corresponds to the negative magnetic anomaly.On the contrary, the positive anomaly "4" in the frontal thrust zone is located close to the positive magnetic anomaly to the south.
Based on the observed correlation of the seismic model with gravity and magnetic anomalies, we propose a mechanism which may explain the variability of crustal thickness along the Nepal Himalayas.As observed in regional magnetic map of India (Fig. S5), the crust of the Indian plate appears to be very heterogeneous.One of the reasons for strong magnetic anomalies may be the presence of large provinces affected by relict igneous processes.Although these occurred in the geological past, the large magmatic intrusions might considerably strengthen the crust.The negative magnetic anomalies, on the other hand, may be explained by thicker sediments.In case of collision, these two types of the crust behave differently and cause different mechanical effects.The crust affected by igneous processes is stronger, and thus the compression of the overlying crust in the collision zone would be more prominent.This would explain the thicker crust in anomaly "4" close to the contact area, which can be explained by a stronger pushing effect of more rigid partition of the Indian plate.The presence of southward curve of the mountain limit line in front of the thicker crustal pattern "4" may represent a broader shortening zone produced by the more rigid incoming block.The segment of the anomalously thinner crust "6" may be explained by a smaller compression rate of the crust because of weaker incoming crust and lubricating effect of thicker sediments having a lower magnetic effect.In this case, the mountain front line is curved northward, indicating less intensive shortening.

Conclusions
In most tomographic studies, the main target is the smooth distribution of seismic properties which is not sensitive to the geometry of first-order interfaces.In this study we made an attempt to reconstruct variations of the Moho interface beneath the Nepal Himalayas where more or less uniform distribution of stations and sources takes place and travel times of both crustal (Pg, Sg) and mantle rays (Pn, Sn) are available.Based on synthetic modeling, we found that, for most of the Nepal Himalayas area, the crustal thickness variations can be robustly retrieved.For the surrounding areas, like northern India and Tibet, crustal structures cannot be resolved with the available data.
The obtained crustal thickness varies from 40 to 75 km along the Nepal Himalayas.There is a fair correlation of the derived crustal structures with the observed gravity and magnetic anomalies.The areas of thicker crust are associated with higher values of the free-air gravity field and vice versa.This correlation is a good argument to prove the reliability of our findings.The magnetic anomalies may provide important information on the mechanical properties of the crust.We see that different segments of the Indian crust behave differently leading to various collision rates.We expect that thicker crust in the frontal thrust zone can be associated with the more rigid incoming crust.Weaker crust segments may penetrate underneath an overlying plate with less resistance, and thus the weaker compression rate leads to thinner crust in the frontal thrust zone.The presence of thick sediments may have a lubricating effect and thus may also reduce the shortening of the crust.
This study gives us a fair understanding of the Moho configuration beneath the central Himalayas.However, western and eastern parts of the Himalayas have not yet been well studied.Such a comprehensive study based on joint consideration of seismic, gravity and magnetic data for the entire Himalayas will make a better understanding of the mechanisms of the India-Asia collision possible.
The Supplement related to this article is available online at doi:10.5194/se-6-207-2015-supplement.

Figure 1 .
Figure 1.Map of the study area and data distributions.Background is topography.Yellow dots are the earthquakes, and red triangles are the stations used in this study.Locations of three profiles used for visualization of the results are shown.Blue lines indicate the Main Central Thrust (MCT) and the Main Boundary Thrust (MBT).IGAP is the Indo-Gangetic Alluvial Plains.Inset shows the location of the region.

Figure 2 .
Figure 2. Distribution of P rays and two parameterization grids corresponding to basic orientations of 0 and 45 degrees.Intensity of the points represents the number of nodes in z direction corresponding to the current x-y coordinates.

Figure 3 .
Figure 3. P velocity and S velocity versus depth.Dotted lines depict starting 1-D velocity models; solid line represents average velocities along section 1 for the main result shown in Figs.6-7.

Figure 4 .
Figure 4. Velocity anomalies at 35 km depth after inversion of synthetic (a) and observed data (b).For the case of synthetic modeling, the limits of "thick crust" are marked with a dotted line.Triangles denote stations.Locations of profiles presented in Figs. 3 to 5 are shown.

Figure 5 .
Figure 5. Synthetic test with reconstruction of the "variable Moho" model in relative anomalies (upper part) and absolute P velocities (lower part).Locations of the profiles are shown in Fig. 2. The configuration of the synthetic "Moho" is indicated with a dotted line.Vertical lines with numbers mark locations where sections cross each other.

Figure 6 .
Figure 6.Vertical sections of the resulting P-velocity anomalies.Locations of sections are indicated in Fig. 4. Above each section, exaggerated topography is shown.Vertical lines indicate locations where sections cross each other.Areas with poorer resolution in sections 2 to 5 are shaded.Moho interface (black line) is traced on the bottom of the low-velocity anomaly.Dashed line in sections 1 and 2 indicates an alternative interpretation which is less plausible.Dotted lines in sections 4 and 5 indicate Moho depth determinations from Hetényi et al. (2006) and Schulte-Pelkum et al. (2005), respectively.More sections are shown in Supplement.

Figure 7 .
Figure 7. Same as Fig. 6, but for absolute P velocity.

Figure 8 .
Figure 8. Map of estimated Moho depth beneath the Nepal Himalayas (a) together with free-air gravity anomalies (Andersen et al., 2010) (b) and magnetic anomalies (Maus et al., 2009) (c).Red numbers indicate the locations discussed in the text.In (a), the locations of profiles used for presenting the main results are shown.