The Mabja Dome Structure in Southern Tibet Revealed by Deep Seismic Reflection Data and Its Tectonic Implications

The ongoing India‐Asia collision has led to the formation of the northern Himalayan gneiss domes belt in southern Tibet. The domes are the result of the ongoing convergence and were formed by geological processes that may include crustal thickening, metamorphism, partial melting, and exhumation of middle crustal rocks to the surface. A combination of compressional, extensional, and diapiric processes has been invoked to explain the formation and evolution of these domes. Differentiating among these competing hypotheses requires well‐defined geophysical images of the internal structure of the domes. The Mabja dome is the largest dome within the northern Himalayan gneiss domes belt. A 70‐km long deep seismic reflection profile across Mabja dome was acquired in 2016. The seismic data could provide new information about the structural elements beneath the domes to the depth of 25 Km. We address the structure of the Mabja dome by conducting an integrated analysis of shallow crustal velocity structure and a normal‐incidence seismic reflection image of deeper crust. Our work suggests that the Mabja dome is underlain by shear zones at depths of 10–15 km and two high‐velocity bodies at depths of 3 km possibly representing the eclogitic‐facies rocks or mafic intrusions. We propose that the dome formation may have been controlled by collision‐induced north‐south shortening expressed by thrust stacking of middle crustal rocks, which led to the doming of the upper‐crustal rocks. The proposed mechanism inferred for Mabja dome can be applied to interpret the widespread domes throughout the southern Tibet and other related structures in orogenic mountain belts.

between India and Asia (Burg, Brunel, et al., 1984;Burchfiel et al., 1992;Diedesch et al., 2016;Hodges, 2000;Watts et al., 2005;Yin & Harrison, 2000;Zeng et al., 2014). Several mechanisms have been proposed to explain the origin and evolution of the domes (Figure 2), these include: diapirism by buoyancy, low viscosity of hot, lateral crustal flow of low-density, middle-crustal rocks and/or low-density magma (Calvert et al., 1999;Ramberg, 1980;Teyssier and Whitney, 2002;Whitney et al., 2004), or middle channel crustal flow (Beaumont et al., 2001;Grujic et al., 2002;Langille et al., 2010;Nelson et al., 1996). Some studies have suggested that east-west extension caused the exhumation of the middle-crust domal footwalls along low-angle normal faults (Brug & Van Den Driessche, 1994;Escuder Viruete et al., 1994). According to Burg, Brunel, et al. (1984), the compressional exhumation of middle-crustal rocks via thrust duplexing and rapid surface erosion played an important role in the development of such domal structures. Lee et al. (2004) suggested that a combination of thrusting and extension may have been responsible for the dome formation. The differing mechanisms that have been proposed for the development of the northern Himalayan gneiss domes ( Figure 2) have substantially different implications for the controlling tectonic regimes (Yin, 2004).
Deep seismic reflection profile is an effective method for probing crustal and upper-mantle structures with relatively high resolution compared with other geophysical techniques. It was successfully applied to the study of the deep structure of southern Tibet in the early 1990s by the INDEPTH project (Brown, 2013;Brown et al., 1996;Nelson et al., 1996;W. Zhao et al., 1993). A new deep seismic reflection profile (Figure 1c) was acquired in 2016 by the Institute of Geology, Chinese Academy of Geological Sciences. The deep seismic reflection profile was specifically acquired to determine the crustal structure and Moho topography below the northern Himalayan gneiss dome. The acquisition and processing parameters (e.g., source-receiver spacing, CDP fold, source energy and normal moveout stretch and mute) of this seismic reflection transect were designed to target the middle and low crust, and thus these were not optimized to image the shallow upper-crustal structures (below 3 km from the surface) which is the critical link between the surface geology and the deeper crustal structures. Nevertheless, the relatively high signal-to-noise ratio of the early arrivals recorded at long offsets can be used to characterize the shallow subsurface velocity structures when employing full waveform tomography (FWT). The combination of the new velocity model of the shallow subsurface structures with the deeper penetrating normal incidence seismic reflection image could provide new information on the detailed crustal structure related to dome formation.
Seismic reflection surveys contain a significant amount of information on physical properties of rocks commonly not used in conventional seismic reflection data processing and interpretations (Pratt, 1999). Here we analyze refracted waves that are often treated as systematic noise to be removed during basic reflection filtering. However, such waves are often useful for static corrections in reflection data processing. For example, ray-based travel-time tomography (TTT) uses refracted arrivals to estimate shallow subsurface velocity variations (Marti et al., 2019;Zelt et al., 2003). Velocity models derived by conventional travel-time tomography feature relatively limited resolution due to approximations used in the ray theory (Zelt & Smith, 1992). On the other hand, the FWT method uses all the information provided by the seismic recordings (e.g., travel-time, amplitude, and phase) rather than the travel-time alone to estimate the elastic parameters of the subsurface crustal materials. FWT improves the resolution and accuracy of the velocity models significantly (Fichtner et al., 2013;Operto et al., 2013;Pratt et al., 2002). The concept of FWT was originally developed in the 1980s for acoustic, isotropic elastic, and viscoelastic cases (Tarantola et al., 1988). The application of FWT was initially restricted to two dimensional (2D) and small-offset data sets because of limited computational resources. With the availability of supercomputing in recent years, the FWT approaches have become a powerful research tool. The FWT result is obtained by minimizing the misfit between modeled and observed data through a gradient optimization method (Fichtner et al., 2013;Operto et al., 2013;Pratt, 1999;Tarantola et al., 1988). The theoretical maximum resolution of FWT is on the order of a half wavelength (Operto et al., 2013;Virieux & Operto, 2009), the typical frequency from the seismic data is about 12 Hz, and the theoretical maximum resolution is about 21 m. Major difficulties in applying FWT include the large amount of computing resources needed, the limitations of recording equipment and survey design during data acquisition. FWT has been applied to marine data successfully (Prieux et al., 2011;Warner et al., 2013;Yelisetti, 2014). There are specific characteristics that facilitate the application of FWT to seismic reflection marine data sets, such as the homogeneous coupling of sensors and sources (i.e., both air-guns and hydrophones are in the water). Furthermore, acoustic wave propagation can be assumed in marine surveys. Source and receiver coupling, scattering, phase conversions, shear-wave propagation, along with multiples substantially hamper the application of FWT on land (Adamczyk et al., 2015;Plessix et al., 2012;Smithyman & Clowes, 2013;Stopin et al., 2014). The relative complexity of the shallow subsurface (strong lateral velocity variations, weathering, etc.) on land complicates even further the inversion of seismic reflection data (Adamczyk et al., 2015;Plessix et al., 2012;Smithyman et al., 2014;Stopin et al., 2014). FWT has largely been limited to 2D data sets with long offsets and, high S/N ratio (Jaiswal et al., 2009;Malinowski et al., 2011). To the best of our knowledge, this study is the first deep seismic reflection data to apply FWT to the geologically complex southern Himalayan region.
To obtain high-resolution crustal structure beneath the dome and study the processes of gneiss dome evolution in southern Tibet. This paper presents the results of new detailed structure of the Mabja dome from high-resolution seismic reflection data. First, we constrain the shallow upper-crustal structures (below 3 km from the surface) by FWT and study the resolution and reliability of the result. Second, conventional seismic reflection data processing was performed to obtain the normal incidence seismic reflection image. Third, the internal architecture of the Mabja dome is obtained by combining the shallow subsurface velocity structure with the normal incidence seismic reflection image. They provide new details, critical information on the evolution of the domal structures jointly in the southern Tibet.

Seismic Data Acquisition
A control-source 2D deep seismic reflection profile across Mabja dome was acquired in 2016. The profile was designed to image the lithospheric structure across the Himalaya (Figures 1 and 3). The survey consisted of a total of 288 shots: 226 small-size shots (96 kg), 60 medium-size shots (192 kg), and 2 big-size shots (1,000 kg). The latter are not used in this analysis ( Figure 3). The 96-kg sources (small size shots) were spaced 250 m apart buried in a single well at a depth of 30-50 m. The 192-kg sources (medium size shots) were spaced every 1,000 meters and consisted of two bore-holes with charges placed at a depth of 50 m. The receiver spacing is 50 m. Both were recorded with geophones spreads that had a minimum offset of 25 m and a maximum offset of 22.4 km. The two large 1,000-kg sources were also recorded for deeper penetration. However, these larger shots were excluded from the current study due to their location, they are too far away to provide coverage overlapping with the smaller shots. A Sercel 428XL seismic recording system was used. The acquisition parameters are listed in Table 1.
Shot records reveal few shallow reflections above 2 s. The lack of observable shallow reflection energy may be due to the lack of appropriate offsets. The relatively high amplitude direct, refracted arrivals were observed in the raw shot records. Several sets of high S/N reflections in the upper crust can be identified (Figure 4). High S/N shot records reveal early arrivals (above 5 s TWT) and long offset (up to 22.4 km). These can be used in FWT approaches to constrain the shallow seismic velocity structure. The seismic velocity LI ET AL.   structure of the shallow crust should provide a critical link between the surface geology and deep structures revealed by the normal incidence seismic reflection data.
Considering that the location of the large size shots ( Figure 1c) is far away from the survey line, this study primarily focuses on the shot records acquired using the small (96 kg) and medium-size charges (192 kg).

Data Pre-processing
The shot records were pre-process using the following: geometry specification and projection on to a 2D straight transect, time windowing, and amplitude scaling to adapt it to the acoustic forward modeling.

Geometry Specification and Projection
Seismic reflection data acquisition is a complex process, especially in an area with high topographic variations such as in the study area. Furthermore, the three-dimensionality needs to be carefully considered when analyzing these types of data sets. True three dimensional (3D) seismic reflection data acquisition is economically unfeasible in a Himalayan region like our study area. However, 2D seismic transects have been proven to be instrumental in unraveling deep-crustal processes in a wide range of tectonic settings (Brown et al., 1982;Carbonell et al., 1996;Hajnal et al., 2005;Pfiffner, 1986;Pfiffner et al., 1988).
The inversion scheme used in FWT involves numerical simulation of seismic waves propagating through a highly heterogeneous media. A 2D acoustic approximation to the solution of the wave equation was used in this study. Waveform inversion tomography requires that the survey geometry and topography characteristics of the survey can be reconstructed accurately when generating synthetic data (Smithyman & Clowes, 2012).
LI ET AL.
10.1029/2020JB020265 6 of 25 Simulating the seismic wavefield is challenging in the southern Tibet due to the irregular geometries that arise from the topography and offline locations of receivers and sources ( Figure 5). The irregular geometries cannot be modeled completely with a 2D inversion method. The exact limit quantifying how much deviation from 2D is acceptable depends on the selection of inversion methods. For 2D FWT, the travel-time errors should be much less than one-half cycle at the frequency of interest, or equivalently the path-length difference should be much less than onehalf wavelength (Smithyman & Clowes, 2012).
The shot and receiver locations were projected onto a straight line keeping the offset intact (Kashubin et al., 2009;Smithyman & Clowes, 2012;Zhang & Christopher, 2013). The source points were moved to a straight line by an orthogonal projection, and the receiver positions were moved to the straight-line profile according to their true source-receiver offset. This approach results in an irregular distribution of receivers along the projected line. Although the offsets are preserved, the distance between the sensors is modified in the process. Some receivers after the projection are very close to one another. In order to reduce the number of receivers, we decided that if the distance between two receivers was less than 50 m, they were treated as located at the same position (Kashubin et al., 2009;Zhang & Christopher, 2013). Figure 5 shows the good agreement between the original source, receiver locations, and the projected 2D locations of the profile used in the inversion procedure.

Time Windowing and Amplitude Scaling
The FWT inversion scheme requires only a limited window of data around the first arrivals that include the target wide-angle reflections and/or refractions. Because the algorithm aims to match the recorded seismic phases, bad traces and noisy seismic data need to be avoided and attenuated, respectively. Thus, it is necessary to remove bad traces, attenuate or eliminate ground roll, abnormal amplitudes, multiples, and scattered energy that probably originate from outside the imaging zone (Malinowski et al., 2011). Obviously, bad traces were removed by visual inspection. To eliminate surface waves and anomalous amplitude bursts (spikes), the data were transformed into the frequency domain and a spatial median filter was applied, frequency bands with amplitudes that deviated from the median amplitude by a specified threshold were re-scaled using neighboring traces. The waveforms were then time windowed to exclude late arrivals contaminated with surface multiples to ensure a more reliable convergence. Spectral whitening was applied to broaden the signal bandwidth. A 4-15 Hz bandpass filter was then applied. The amplitude of the data was further scaled by sqrt (t), where t is the twoway travel-time, to account for the fact that our models are strictly in a 2D space (Adamczyk et al., 2014) and partially correct for spherical divergence effects. The window size was increased in later stages of the inversion. A time window of 3.2 s was chosen to include as many arrivals at far offset as possible for the final inversion.

Travel-Time Tomography
The FWT used in this study is based on finite-difference acoustic wave equation solver in the frequency domain with a smooth topography. The topography is a critical parameter in the study area ( Figure 3) since the elevation varies from 3.5 to 5.0 km. A floating datum static correction has been applied to smooth the topography and minimize its effect on the inversion scheme. Inversion was achieved by minimizing the L2 norm and calculating the misfit between the synthetic and field data. A conjugate gradient algorithm guides the iteration toward the minima of the inversion function. The gradient is computed by multiplying the LI ET AL.  forward propagated wavefield with the back-propagated data residual at each time step. The initial starting model and the frequency selection are very important to stabilize the FWT process and obtain a reliable velocity structure (Koehn et al., 2012).
A relatively accurate initial model is required to avoid cycle-skipping, particularly for inverting higher frequencies since the FWT is a nonlinear process. A reliable initial model is usually obtained by first arrival TTT. The TTT scheme has been proven to provide relatively reliable models for shallow crustal surveys (Marti et al., 2019). About 175,712 first-arrival travel-time within the offset range of 50-24,000 m were manually picked from the raw data ( Figure 6a). Large reciprocal errors are often caused by geometry and/ or picking errors. Therefore, the reciprocal error is used as quality control of the geometry and travel-time picks. Computation of reciprocal errors leads to about 15% low confidence (with reciprocal errors beyond 100 m s) travel-time to be removed before inversion.
The frequency bandwidth present in the early arriving waveforms was approximately 4-16 Hz. For the travel-time tomography, a one dimensional (1D) velocity model (Figure 6c) was established based on the refracted velocity from fitting straight lines with different slopes to the first arrival travel-times ( Figure 6b). The TTT is based also on projected geometry with a grid spacing of 25 × 12.5 m. The travel-time residuals decreased from 160 to 12 m s after 12 iterations of the TTT. The synthetic travel-time from the TTT model matches the real picked first-break travel-time reasonably well (Figure 7a). Only a very minor decrease in the residuals takes place after 12 iterations (Figure 7b). Therefore, the inversion was stopped at the twelfth iteration. An average residual for all the shots of 12 m s was obtained. This was an acceptable accuracy, as the picking uncertainty is estimated to be on the order of 14 m s for this profile, and the data residuals are approximately evenly distributed. The ray density plots indicate that the TTT velocity model (Figure 8a

Full Waveform Tomography
The analysis of the shot gathers reveals that the lowest frequency of the first arrivals is about 5 Hz. This is most probably a result of the combination of the sensor hardware (10 Hz geophones natural frequency), the relatively large dynamite charges used as seismic sources and spectral whitening. As demonstrated by previous studies, the low frequencies are very important in FWT (Bleibinhaus et al., 2009;Pratt, 1999), because they are required to recover long-wavelength structure. Without the low frequencies, it is difficult to achieve convergence during inversion. The higher the starting frequency used in the FWT inversion, the closer the starting velocity model needs to be to the hypothetical solution. Bleibinhaus et al. (2009) developed a scale-independent method to measure the resolving power at the lowest frequency by relating the model size to the Fresnel zone width. According to this study, a small resolution number (N LR ) corresponds to a wide Fresnel zone width that indicates good resolving power at long wavelengths. A more accurate starting model is required for higher resolution numbers. In our case, the average P-wave velocity is about 6 km/s, which gives a resolution number of 6.9 at the lowest available frequency of 5 Hz. Previous studies suggested that N LR varies between 0 and 7 for 0.5-130 Hz (Bleibinhaus et al., 2009). The value obtained here is on the higher end of this range. This large value corresponds to a narrow Fresnel zone, indicating that non-linearity is relatively high and a relative sophisticated initial model is required for waveform inversion (Bleibinhaus et al., 2009). Hence, we conclude that the TTT model is a suitable starting model for waveform inversion.
The inversion frequencies were divided into two overlapping groups with increasing frequency. Each group consisted of four frequencies. The first group had frequencies of 5-10 Hz with an interval of 2 Hz and a grid size of 25 m, whereas the last group had frequencies within 9-15 Hz with an interval of 2 Hz and a grid size of 12.5 m. The adjoint gradient scheme for each frequency group, which consists of forward propagation, backpropagation, and a multiplication of the two fields was computed by using all the frequencies simultaneously. Inversion over five iterations was conducted for the first group, and the results were used for the initial model of the second group. A total of eight iterations were performed for the second group. The strategy of moving up through the frequency spectrum helps mitigate the nonlinearity of the problem (Table 2) This method is tolerant to velocity errors for lower frequencies, which are less likely to produce errors of more than a half cycle in the waveforms. As the inversion proceeded, the velocity model was gradually improved as the FWT progressively moved to higher frequencies. Figures 9a and 9b show the velocity models from FWT using frequencies of 5-10 and, 5-15 Hz, respectively. With the addition of higher frequencies and the increasing number of iterations, the resolution of the model was improved gradually. The final velocity model of the FWT is presented in Figure 9b. The details and resolution of the velocity model are significantly improved compared to the TTT model (Figure 8a).
The FWT results need to be validated owing to the non-linearity of the objective function and the risk of falling into local minima (Adamczyk et al., 2015). Consistency checks must be used to validate the obtained result since there are no boreholes with sonic-logs or rock samples available in the area.

Table 2
The Frequency Group Used in the FWT The 1D P-wave velocity profiles at various locations along the profile are shown in Figure 10. The TTT and FWT results have similar trends. However, the velocity changes from FWT are more abrupt and irregular than those associated with the TTT results ( Figure 10). Figure 11 shows the difference between TTT and FWT results. The FWT results contain more small-scale heterogeneities. A low-velocity anomaly between 3 and 4 km depth is resolved by both TTT and the FWT results. Velocity anomalies are better defined by the FWT, and velocity anomalies feature higher values (up to 250 m/s) when compared with the TTT velocity results (Figure 11). The most prominent features in the FWT results correspond to two high-velocity anomalies localized at depths of 2.5-3.5 km, and at horizontal distances between 30 and 40 km and between 55 and 64 km, respectively, along the profile. For the anomaly between 52 and 65 km, the FWT results indicate very high-velocity values approaching 7.2 km/s. Shot records reveal relatively high amplitude multi-cycle reflection event at the same location ( Figure 12). This high amplitude multi-cycle event correlates spatially with the velocity anomaly.
The pre-processed field data ( Figure 13) with the synthetic data (Figures 13b, 13d, and 13f) generated using the FWT model are displayed at the different locations along the profile. The field data generally match the synthetic data and the similarity coefficient is between 0.78 and 0.89, indicating that the FWT model is preferred by the seismic data, the FWT model is closer to the true subsurface P-wave velocity structure than the TTT model. The modeled first arrival times are within half a cycle with those observed, especially at shallow depths. Minor differences appeared between the modeled and raw data. These differences may be due to imaging 3D structure with a 2D profile, and the elastic effects which are not accounted for the acoustic wave equation solver.
LI ET AL.

Normal Incidence Seismic Reflection Image
Conventional seismic reflection data processing was performed using Focus and ProMAX commercial seismic processing packages (ProMax 5000 and Omega 2015 user guide). A straight line that is the same as the FWT geometry was defined with a 25 m CMP spacing. Bad/noisy traces were hand-picked and removed. Static corrections were estimated by using the FWT model (Figure 9b). A floating datum (the same datum used for the FWT inversion) was used to increase the quality of velocity analysis and the stacked image.
Noise attenuation was carried out by bandpass filtering and abnormal amplitude attenuation with different frequency bands. The time-domain data were transformed into the frequency domain, which were then processed using a spatial median filter with amplitudes that deviate from the median amplitude by a specified threshold. Amplitudes were then compensated for spherical-divergence-spreading. Surface-consistent-amplitude corrections were computed and applied to compensate for abrupt vertical and horizontal amplitude variations. Surface consistent deconvolution was also applied, using a 32 m s prediction length, 180 m s operator length, and 0.1% white noise to suppressed reverberations and improved the vertical resolutions. A raw and processed shot gather is shown in Figure 14.
LI ET AL.

10.1029/2020JB020265
12 of 25    were applied sequentially. Kirchhoff post-stack depth migration (0-6 s) was performed with a migration aperture of 6 km and an anti-aliasing frequency of 50 Hz, the smooth FWT and migration velocity analysis together provide an optimal imaging velocity model for migration. A 1,000 m s automatic gain control (AGC) and band pass filter (6-8 to 45-50 Hz) were applied before the final display. The migrated seismic profile showed many reflective events. Several sets of discontinuous arcuate events are at a depth of 1-7 km. There are some relatively weak reflections among the arcuate amplitude events. A sub-horizontal series of high amplitude reflections is resolved at a depth range of 10-20 km. The main seismic reflection data processing flows and parameters are listed in Table 3.

Results
Imaging 3D structure with a 2D surface survey usually results in complex images that can include out-ofplane high amplitude reflections that cannot be readily distinguished from those originating beneath a seismic line. This may increase the uncertainty of the P-velocity structure and seismic reflection structure (Calvert, 2016). Furthermore, the dips imaged by the stack section differ from the true dips of the geological strata. The dip of stratum on the reflection profile determined by the true dip needs to correct for the strike of subsurface structure. These issues may limit our ability to interpret the P-velocity structure and seismic reflection image.

P-Velocity Structure
When the final FWT velocity model is overlaid on the migrated seismic reflection image (Figure 15), several correlations can be observed. The velocity results can be considered in terms of four zones.
We interpret a relatively low-velocity (1.2-2.6 km/s) zone nearest to the surface as a layer of sedimentary cover. This relatively low-velocity layer can be further divided into two parts by a normal fault that correlates with the surface geology. The shallow low-velocity layer likely corresponds to the Quaternary sediments between 13 and 42 km (southern end of the profile  the profile), the shallow low-velocity layer likely corresponds to the Jurassic marine sediments and volcanic rocks (Figures 1c and 15b).
Underlying these inferred sedimentary units is the second zone, labeled THS, which is characterized by velocities up to 5.5 km/s ( Figure 15). The THS zone is relatively continuous from south to north, with a velocity ranging from 4.0 to 4.7 km/s. The relatively small velocity variation within this zone indicates a relatively homogeneous material. The velocity values would be consistent with those expected for the Tethyan Himalayan Sequence, which consists of sedimentary rocks that underwent low-grade metamorphism during the India-Asia collision.
A third zone extends from about 2.3 to 3.3 km (between the red and blue dotted lines, with the top labeled by STDS in Figure 15b). This zone is laterally heterogeneous with velocity increasing with depth from 5.0 km/s to values over 6.0 km/s. Two lensed shaped high-velocity anomalies (labeled H1 and H2 in Figure 15b) can be found within this zone. H1 is defined by a velocity of ∼6.6 km/s and H2 defined by velocities of LI ET AL. >7.2 km/s. The velocity is very high for crustal material, which may have important indicative significance for understanding the formation and evolution of the dome.
A fourth zone characterized by velocities starting at ∼5.9 km/s and reaching values up to 7.0 km/s extends from depths of 1.0-2.3 km. This relatively high-velocity layer correlates with the Greater Himalayan Crystalline Complex (GHC).

Seismic Reflection Structure
The migration of the seismic reflection data (Figure 16a) reveals distinct reflection patterns that provide insights on the structure of the upper crust to a depth of 20 km from the surface. Arcuate events R1, R2, and R3 are the highest amplitude reflections in the section. Amplitudes for the bright spots range from eight to 14 dB above the background in the INDEPTH seismic reflection profile , whereas amplitudes for the rest generally average about 4 dB higher than the background along the profile. The continuity and amplitude of R2 is weaker than R1 and R3. The reflections of R1 and R3 overlap each other at 15 km. The arcuate geometry R4 is consistent with interfaces that have been deformed by doming and, thus suggests that that the deformation that caused doming extends to depths of at least 12 km. The dip of R4 on the migrated reflection profile considering the true dip and strike of subsurface structure would be smaller than the true dip if the profile is not perpendicular to the strike of these structures.
The STDS is a down-to-the-north, low-angle detachment system that is traceable along the length of the Himalaya (Burchfiel et al., 1992;Burg, Guiraud, et al., 1984). The STDS represents an important physical boundary and outcrops at Kangmar dome  which feature places generally low-grade Tethyan metasediments against the Greater Himalayan gneisses (Burchfiel et al., 1992;Burg, Guiraud, et al., 1984;Edwards et al., 1996;Hodges et al., 1996). The STDS become shallow and exposed around NHGD and the STDS reflection event around Kangmar dome ( Figure 17f) has been imaged as relatively continuous and high amplitude reflection events at relatively shallow depth along the INDEPTH seismic reflection profile W. Zhao et al., 1993). Considering the geometry, the velocity jump from FWT results, and the depth of STDS reflection around the Kangmar dome, the relatively continuous high amplitude event (R5) that correlates with a velocity jump in the FWT results can be interpreted as the STDS reflection. It may also be seen that R5 has a similar spatial location and geometry compared with that of the Kangmar dome (Figure 16b). We interpret the R5 to represent the lower boundary of the STDS shear zone.

Discussions
The geological structure, identified from the seismic reflection image and the FWT velocity model, provides important constraints on the tectonic origin of the Mabja dome. The Mabja Dome, north of Kanchenjunga, is the largest dome of the NHGD belt. The NHGD is described as a discontinuous, east-west-striking antiform structure (Jessup et al., 2019) outcropping between the STDS to the south and the Indus-Tsangpo Suture Zone to the north (Figure 1), The STDS become shallow and exposed around NHGD which features a range of different degrees of metamorphism, from a gneiss core beneath high-grade metasedimentary rocks, to weakly metamorphosed to unmetamorphosed sedimentary rocks on top of the core (Burge, Guiraud, et al., 1984;Nelson et al., 1996) places the NHGD at the tip of an interpreted, major change in dip angle of the subducting Indian plate. The INDEPTH seismic reflection image reveals a complex compressional structure beneath the Kangmar dome ( Figure 17). Other geophysical and geological studies suggest a major lateral change across the NHGD from south to north (Jessup & Cottle, 2010;Jessup et al., 2006Jessup et al., , 2019Le Fort, 1975;Le Fort et al., 1987;Nelson et al., 1996;Unsworth et al., 2005;W. Zhao et al., 1993). LI ET AL.   Nelson el al. (1996). (a) Geology outcropping along the profile. (b) Velocity model derived from the FWT with the proposed interpretation of the high-velocity anomalies (in hot colors). (c) Normal incidence migrated seismic reflection revealing the major structural elements, the STDS which directly correlates with the velocity isoline identified in (b) and marked as STDS. Note that the STDS is at outcrop a high sheared structure. (d) The INDEPTH seismic profile . (e) Sketch drawing with the proposed interpretation. (f) Generalized crustal-scale cross-section revealing the major internal architecture of the Himalayas orogen from the INDEPTH profile . The blue square marks the location of the studied area.

The High-Velocity Body
The FWT velocity model reveals two shallow high-velocity anomalies H1 and H2 (Figure 15b), their velocities are higher than background average velocity ∼5.9 km/s. H1 is defined by a velocity of ∼6.6 km/s, and H2 is defined by velocities of >7.2 km/s. The velocity is very high for crustal material.
These velocities anomalies would result in relatively high reflection impedance contrasts If these velocity changes occur over a sufficiently small vertical distance, they would be consistent with the high amplitude reflections seen on the normal incidence seismic reflection image (Figure 15b). Shot records reveal relatively high amplitude multi-cycle events that correspond to the velocity anomalies (Figures 12 and 13). The final migrated image also reveals high amplitude multi-cycle events ( Figure 15) located at approximately 2-3 km in depth that correlate with the high-velocity anomalies.
The juxtaposition of the normal incidence seismic reflection image and the FWT velocity model has revealed a correlation between the high-velocity anomalies and the high amplitude multi-cyclic reflection fabrics ( Figure 15). A number of factors may contribute to generating high amplitudes reflection events in seismic reflection sections: (1) the variations of physical properties of the media (Seismic wave propagation velocities, densities); (2) the angle of incidence of the incoming wave; (3) the wavelength of the source (or incoming) signal at the target depth; (4) the vertical dimensions of the structure (thickness of layering); and (5) the horizontal dimension of the structure (horizontal extent of the layering). All these factors make it difficult to quantitatively infer lithology from observed reflection amplitudes and polarity. With these considerations, the dynamic characteristics of the reflected wavelets are difficult to estimate in crystalline and highly deformed terranes. The relation of the horizontal dimensions of the reflecting structure with respect to the Fresnel zone is a critical consideration (Hubral, 1983;Hubral et al., 1993;Eaton et al., 1991). It is also well known that the constructive interference effect due to layering can greatly affect reflection amplitudes (Poppeliers & Levander, 2004;Warner, 1990). This constructive interference effect is highly significant when imaging finely layered metamorphic structures (Holliger, 1997;Holliger et al., 1992), or mylonitic shearzones (e.g., the Lithoprobe transects across the Grenville province, Canada [Ji et al., 1997]).
Given the factors discussed above, there are several plausible interpretations for the high amplitude reflection associated with bulk velocity anomalies. For example, the observed seismic signatures would be consistent with a layered sequence of rocks characterized by alternating velocity values (high and low). The average value of the velocities within the layered structure needs to be consistent with the values of the velocity anomalies mapped by the FWT. A layered structure with a high velocity may represent mafic igneous rocks or ultrabasic rocks derived from partial mantle melting, and/or metamorphic crustal rocks that underwent different degrees of metamorphism. Lab studies of eclogitic-facies rocks have shown that they have a velocity in the range of 7.0-8.4 km/s. Eclogitic-facies rocks (Austrheim & Griffin, 1985;Fountain et al., 1994;Kern et al., 1999) would be consistent with the observations. Wang and Zhang (2015) documented granulite and eclogite-facies rocks in selected Himalayan gneiss domes, which were originated from Indian plate thrusting beneath the Asian continent followed by exhumation during doming to the surface (Cottle et al., 2009;Jessup et al., 2008;Wang et al., 2013;J. Zhang et al., 2012). However, mafic rocks in the form of dikes, sills, or relatively thick mafic intrusions cannot be ruled out as a cause of both the high bulk velocities and high amplitude reflections (Cook et al., 1999;Simancas et al., 2003). Laboratory measurements of physical properties and synthetic seismic studies have demonstrated that relatively high reflection coefficients and high velocities can be achieved by mafic layered intrusions alone (Deemer & Hurich, 1994;Simancas et al., 2003). We propose a mixed origin for the high velocities and high amplitude reflections, consisting of intrusion of mafic magmatic bodies of mantle origin into a structurally controlled level, probably a main decollement (Simancas et al., 2003).
The boundary between the broader low and high velocity (R5) indicated by a red dotted line in Figure 15b could plausibly correspond to the STDS. Many studies of the STDS infer that all the observed displacement occurred along a single fault zone during a deformation event that progressed along a low-angle shear zone that uplifted and exposed the domes at the surface (Yang et al., 2009;Yin, 2006;Zhang et al., 2007). The STDS separates low-grade metamorphic THS from high-grade GHC, marking one of the critical geological boundaries in the Himalayas (Figure 1) (Burge, Guiraud, et al., 1984;Zhang et al., 2007). The depth and geometry of the STDS may have been altered by the doming process. The R5 is consistent with the downward projection of the mapped STDS at the surface. The boundary (R5) is most likely the lower boundary of the STDS shear zone, which uplifted and exposed the domes at the surface.

The Properties of the High Amplitude Reflectors
Deep arc-shaped events are also observed at depths of −10 km to −15 km (R1, R2, and R3 in Figure 16). These seismic reflections have locally anomalous high amplitudes. The near horizontal and laterally continuous reflectivity may be explained in two ways. First, R1, R2, and R3 could represent laminated foliations in a shear zone. These reflections have been observed with similar characteristics in the Kangmar dome area in the INDEPTH deep seismic reflection image Makovsky et al., 1996;Nelson et al., 1996;W. Zhao et al., 1993). The arcuate reflector R4 above horizontal reflectors suggest that R1 and R3. R1, R2, and R3 are parts of a thrust decollement possibly representing the MCT (Figure 17d). The reflectors were offset by younger duplex thrusts leading the formation of the gneiss domes. Second, R1, R2, and R3 may be "bright spot" features interpreted as partially molten layers by Nelson et al. (1996), or a laminated structure with a high content of water (fluid) as suggested by Makovsky and Klemeprer (1999).

The Deeper Structure (3-25 km)
Our new seismic reflection data cover the transition zone of the NHGD (solid black line in Figure 1c). Seismic refraction recording by the reflection survey was used to estimate the shallow (e.g., depth < 3 km) velocity structure by FWT. The combined shallow velocity model and the deep seismic reflection image provide new details on the upper crustal structure of the transition zone. The seismic reflection image reveals a series of sub-horizontal reflections (R1, R2, and R3) at depths of 11-20 km ( Figure 16b). R1, R2, and R3 are parts of a thrust decollement possibly representing the MCT. MCT is defined by a thick shear zone of a few kilometers to 10 km thick and accommodated a minimum of 140 km of displacement (Yin & Harrison, 2000) which experienced multiple episodes of shortening (Kellett et al., 2009;Yin & Harrison, 2000;Z. Zhang et al., 2012). We interpret that the sub-horizontal reflections represent the MCT shear zone that has been deformed by thrust stacking from below (Z. . The stacked thrusts caused local crustal thickening expressed as doming. The proposed thrust stacks (R1, R2, and R3) at depths of 11-20 km favor discrete thrust shear-zone deformation as the cause for dome formation (Figure 16b).
The R5 (Figure 16b) is a relatively high amplitude reflection event in the depth migrated profile using the FWT velocity model (Figure 15b). The R5 reflection is most likely a marked shear zone -the STDS. The P-velocity model from FWT and seismic reflection image (R5) shows a similar pattern for the STDS. Both profiles could be combined by their common interface along the STDS (Figure 17). The shallow part is characterized by a P-velocity model with shallow dipping reflection structures (Figure 17b), whereas the deeper part is a structurally complex image in the depth migrated reflection profile (Figure 17c). The above results suggest that the deformation of the dome was caused by N-S compression, which resulted in thrust stacking along the MCT at a depth of ∼15 km from the surface (Figure 17e). The depth between 12 and 16 km is considered to correspond to the brittle to ductile transition according to the rheology of a typical continental crust (Aharonov & Scholz, 2019).

Evolution of the Dome and Its Implications
The above interpretation is consistent with the parallel geometry and possible relationships between the gneiss dome belt and the compressional Himalayan orogen. However, other processes such as metamorphic core complex-type extension, diapirism, or duplex formation have been proposed as the mechanism for the formation of the north Himalayan gneiss domes (Calvert et al., 1999;Ramberg, 1980;Teyssier and Whitney, 2002;Whitney et al., 2004).
The reflections in our seismic image reveal a complex compressional structure beneath the Mabja Dome (Figure 17c). The structure of the dome exhibits arcuate stacked thrusts. Burg, Brunel, et al. (1984) have interpreted the NHGD domal form as the result of a southward thrust duplex. The thrust-duplex could be the critical structures that controlled the formation and evolution of the dome (Figure 2d). Based on the new geophysical images of the internal crustal structures of the dome, we suggest that N-S compression created thrust structures that deformed the MCT at a depth of 15 km (Figure 17e). The compression associated with thrust stacking was accommodated by doming of the upper crustal strata (Figure 17e). Melts rising through fractures and shear zones were emplaced as relatively thin sills (H1 and H2) at shallow crustal levels. Doming, in turn, resulted in the extension of the upper crust (Figure 17b), the normal fault (Figure 17c) constitutes evidence for this shallow extension in THS. North-south compression at depths of 15-20 km and north-south extension above a depth of 10 km may have governed the dome formation. The N-S compression and shallow extension contributed jointly to the structural evolution and formation of the dome.
The Mabja Dome is the largest of the North Himalayan dome gneiss domes (NHGD) and has similar structural and metamorphic histories with other domes in the southern Tibet (Lee et al., 2004;Jessup et al., 2019). The proposed mechanism inferred for the Mabja dome can be applied to interpret the widespread dome development throughout the southern Tibet and other domes in orogenic mountain belts.

Conclusions
A detailed study has been carried out to constrain the tectonic evolution of the northern Himalayan gneiss dome belt. The dome belt is located at a tectonic boundary between continental subduction to the south and partial crustal melting to the north. The new seismic data provide constraints on the structural geometry of the Mabja in the belt. The shallow crustal velocity variations were constrained by a full waveform inversion of refracted seismic waves recorded in a recently acquired seismic reflection survey. A corresponding deep seismic reflection image provides a view of major structural features beneath the dome. Our results lead to the following conclusions: (1) Two high-velocity bodies, with an average velocity as high as 6.6 and 7.2 km/s, are constrained at a depth of 3 km from the surface. These velocities are consistent with the presence of eclogitic-facies rocks or mafic intrusions (2) A series of sub-horizontal reflections (R1, R2, and R3) are observed at depths of 11-20 km which are interpreted to represent a thrust system whose local development led to crustal thickening. R1, R2, and R3 are parts of a thrust decollement possibly the MCT. The arcuate geometry of R4 has been deformed by doming and suggests that the deformation that caused doming extends to depths of at least 12 km. The R5 reflection is most likely the lower boundary of the STDS shear zone, which uplifted and exposed the domes at the surface (3) The N-S compression created thrust structures that offset the Main Himalayan Decollement at a depth of ∼15 km. The compression was associated with thrust stacking in the middle crust and doming in the upper crust. Combined N-S compression and extension may have contributed to the formation and structural development of the seismically imaged dome. The thrust stacking in the middle crust is the main mechanism that controls the formation of the Mabja dome formation (4) The Mabja Dome is the largest dome within the NHGD belt and has similar structural and metamorphic histories as the other domes in the southern Tibet. The proposed mechanism inferred for the Mabja dome can be applied to interpret the widespread dome throughout southern Tibet and other dome structures in orogenic mountain belts

Data Availability Statement
The data in this study will be available at this site https://osf.io/snrp6/.