Constraining the Moho Depth Below Bhutan With Global-Phase Seismic Interferometry

We use a novel technique named global-phase seismic interferometry (GloPSI) to image the lithospheric structure, and in particular the Moho, below two parallel north-south transects belonging to the GANSSER network (2013–2014). The profiles cross the Himalayan orogenic wedge in Bhutan, a tectonically important area within the largest continent-continent collision zone on Earth that is still undergoing crustal thickening and represents a challenging imaging target for the GloPSI approach. GloPSI makes use of direct waves from distant earthquakes and receiver-side reverberations with near vertical incidence. Reflections are isolated from earthquake recordings by solving a correlation integral and are turned into a reflectivity image of the lithosphere below the arrays. Our results compare favorably with first-order features observed from a previous receiver function (RF) study. We show that a combined interpretation of GloPSI and RF results allows for a more in-depth understanding of the lithospheric structure across the orogenic wedge in Bhutan.


INTRODUCTION
The temporary passive seismic experiment GANSSER (Geodynamics ANd Seismic Structure of the Eastern Himalaya Region) continuously recorded seismic data from a total of 37 broadband stations in Bhutan from January 2013 to November 2014. It was the first seismological project to address the lithospheric structure across the eastern third of the Himalayas. This is of interest because the segment beneath Bhutan is the only segment along the Himalayan collision belt where a basement uplift in the foreland, the Shillong Plateau, reaches high elevation, as a result of exhumation since the Pliocene (Biswas et al., 2007;Najman et al., 2016), and the region features mostly strike-slip mechanism earthquakes (Drukpa et al., 2006), unlike the Central Himalaya. Bhutan's structural segment boundaries within the Himalayas are highlighted both by gravity anomaly residuals and by orogen-crossing mid-to lower crustal seismicity of the underthrusting plate (Hammer et al., 2013;Hetenyi et al., 2016a;Diehl et al., 2017). Both the Bhutan Himalayas and the Shillong Plateau can accommodate great earthquakes, as witnessed by the 1714 and 1897 events, respectively (Bilham and England, 2001;Berthet et al., 2014;Hetenyi et al., 2016bHetenyi et al., , 2018Le Roux-Mallouf et al., 2016). However, as demonstrated by Coulomb stress modeling, the two regions are seismotectonically less linked than previously proposed (Grujic et al., 2018).
The structure of the orogenic wedge in Bhutan features the Main Himalayan Thrust (MHT), along which the underthrusting plate slides beneath the Himalaya. This is similar to the rest of the mountain belt, however, both the Moho and the MHT geometry present lateral variations between West and East Bhutan, over a distance of ca. 180 km (Singer et al., 2017a), as well as when compared to other cross-sections across the orogen (Subedi et al., 2018). On receiver function (RF) images by Singer et al. (2017a), in the West (along 89.5-90 • E) the Moho is at ca. 50 km depth at the front of the orogen in the south, and turns into a ramp ca. 70 km further north. The MHT is less clear, but is at roughly 15 km depth and turns into a ramp, in such a way that the underthrusting crust keeps a constant thickness. The trend at depth of the MHT, and the consequent separation between the over and underthrusted crusts is largely consistent with the thermokinematic modeling results of Coutand et al. (2014). In the East (around 91.5 • E), the Moho is at 50 km depth at the front of the orogeny as well, but appears sub-horizontal beneath Bhutan, and steepens gently farther north to connect to deeper Moho segments beneath southern Tibet (Shi et al., 2015). The MHT is somewhat shallower than in the West, being located at ca. 12 km depth. The intra-crustal microseismicity also varies laterally: while it is well focused along the MHT in the East, it is more diffuse and is partly "shifted" below the MHT in the West (Diehl et al., 2017), which corresponds well with the pattern of interseismic coupling obtained from GPS observations (Marechal et al., 2016). Both, the Western and the Eastern line RF images include intra-crustal converters beneath the MHT, which are not constrained by surface geology. These features can result from the deformational Himalayan orogeny (Singer et al., 2017a) or may be inherited elements of the underthrusting plate or terrain. In the upper part of the orogenic wedge, ambient noise tomography from the same seismic network data shows shear-wave anomalies that are well correlated with the tectonic structures known from structural geology at the surface (Singer et al., 2017b).
Here, we exploit data from the GANSSER network for an application of the global-phase seismic interferometry (GloPSI) method (Ruigrok and Wapenaar, 2012;Frank et al., 2014;Bianchi et al., 2021). The main aim is to image the Moho below the two parallel north-south transects crossing the Himalayan orogenic wedge in Bhutan, based on P-wave reflectivity, and to compare it with the RF observations (P-to S-wave conversions). GloPSI makes use of phase arrivals and their receiver-side reverberations from earthquakes at epicental distances of 120-180 • (Figure 2). These phases impinge on the lithosphere from below with small angles of incidence. The method hereby assumes horizontal and slightly dipping reflectors that extend at least as far as the diameter of the Fresnel volume. The maximum dip that can be imaged with GloPSI, is heavily dependent on the illumination. If there is "full" illumination with global phases, that is good coverage from ray-parameter p = 0 until p = 0.04 s/km, then the maximum dip that can be imaged is around 14 • for reflectors in the lower crust [phi = arcsin(pv), where v is the velocity in the lower crust]. The reflections are isolated from the earthquake recordings with auto-correlation and stacking. For more details on the SI theory illustrated with modeling, please refer to, e.g., Frank et al. (2014), Draganov et al. (2006); Ruigrok et al. (2008), and Ruigrok et al. (2010). With a migration step, the resulting reflection responses are turned into a reflectivity image of the lithosphere below the arrays. In this way, we use a different, independent waveform set to study the Moho. The challenge is twofold: how does GloPSI perform in an environment of relatively thick crust, with traces and imprints of orogenic deformation, and how feasible is it to retrieve useful signals from a relatively small-aperture profile with a small number of, but densely spaced stations, compared to typical GloPSI processing.

Data
We use 16 months of seismic data records collected from the temporary GANSSER network in Bhutan (data 1 ) between January 2013 and November 2014. The stations operated in continuous recording mode, at 100 samples per second, and data recovery rate exceeded 90%. To work along densely spaced profiles, we selected data from the 28 three-component broadband seismometers that were set up along two almost straight, northsouth arrays across the Himalayan orogenic wedge in Western (16 stations along a 116 km deployment) and Eastern (12 stations along a 98 km deployment) Bhutan, respectively ( Figure 1A). Station spacing was as equidistant and linear as access and local conditions allowed. All stations operated for 16 months, and about a third of the stations for the full 22 months.
In the raw seismic records, we looked for direct waves that reach the GANSSER array while traveling through the Earth's core. We are therefore interested in large magnitude M > 6 seismic events occurring between 120 and 180 • epicentral distance, at any depth range (Figures 1B, 2). For this purpose, we created a list of events occurring within the region of interest from January 2013 to November 2014 from the EMSC website 2 . We then visually selected waveforms with clear P-wave arrivals. For this purpose, we plotted each event for all stations and discarded unclear traces or entire events. We selected 59 earthquakes with M > 6 (Supplementary Table 1), in the desired epicentral distance range ( Figure 1B). For each earthquake we selected a time window of 50 s in length, starting at the first arrival (PKiKP, PKP, or PKIKP, Stein and Wysession, 2009) and containing subsequent reverberations in the lithosphere. To limit the occurrence of cross-terms in the autocorrelations that might appear if PKIKP and PKP phase responses are simultaneously present, we use only time windows with one single dominant phase response.

Method
The selected phase responses have a good signal-to-noise ratio and at most, only four traces are missing per event. We can detect PKP over a large distance range, which should enable a stack of sufficient illumination angles to correctly estimate the timing of the main zero-offset reflections. The different data processing  Hauck et al., 1998). (B) Selected earthquakes (yellow stars) at epicentral distances of 120-180 • (red dashed lines) from Bhutan (red triangle) with a sufficiently high signal-to-noise ratio. We found 59 earthquakes with M > 6.0 that satisfy these criteria.
steps are subsequently described and illustrated below in Figure 3 for Western (left) and Eastern Bhutan (right).
The first step is the application of seismic interferometry, which for GloPSI, amounts to the autocorrelation and stacking of selected time windows, and the flipping of the amplitudes (Figures 3 A1, B1). The Z-component records are used as a proxy for the P-wave transmission responses. For each phase, we apply a series of processing steps before stacking the autocorrelations. After loading an instrumentally deconvolved phase response, a mild spectral whitening is applied. This whitening function is implemented with a sliding frequency window. The large earthquakes that we typically use have a very low corner frequency. At larger frequencies, there is a significant drop in amplitudes. The whitening is thus necessary to boost the higher FIGURE 2 | (A) Illustration of raypaths (black lines) between a source (black star) and a receiver (green triangle) for the most prominent P phases. The typical epicentral distances for GloPSI and RF are marked with shading. (B) Close-up of the medium below the receiver, with a global phase wavefront just before it impinges on the lithosphere. (C) Illustration of exemplary zero-offset reflection responses used in GloPSI (modified from Ruigrok and Wapenaar, 2012).
frequencies and to construct a reflection image with sufficient resolution. The upper limit can be selected at a frequency where there remains a sufficiently high signal-to-noise-ratio for all events. In our case, we applied a spectral whitening in the range of 0.01-1.4 Hz, with a sliding frequency window length of 0.12 Hz. The whitening from very low frequencies onward, boosts the direct P-wave pulse, which can be easily removed with a taper (Figures 3 A2, B2). Since this part of the signal is of a much higher amplitude than the following wiggles, the signal of interest with relatively lower amplitudes can be better seen after removal. Subsequently, a bandpass filter is applied between 0.2 and 1 Hz in order to remove the low frequency sidelobes and consequently discard the red-to-blue difference from crust to mantle that could be observed in Figures 3 A1, B1. Similar to the HiCLIMB dataset across Tibet (Ruigrok and Wapenaar, 2012), there were not enough phase responses to stack out the source-side reverberations completely. They were suppressed by substracting the scaled average trace over the array from each individual trace. This processing step introduces collateral damage which suppresses actual horizontal reflectors, if such reflector exists over a large part of the array. Hence, this step is not to be applied to the data for the interior of the continents. Average removal was not needed in Bianchi et al. (2021) who applied GloPSI to an array that was deployed for a similar duration.
We then apply a static correction and eliminate free-surface multiples at the receiver side (step3, Figures 3 A3, B3) following the method from Verschuur and Berkhout (1997). There could also be triplications of the PKP phase from the source side, or possibly over the D discontinuity. As long as the different triplications have a similar ray-parameter, their detrimental effect is also suppressed with this step. However, if there were multiple strong arrivals in one phase response with different angles of incidence, such a panel would need to be excluded from further analysis. Since no strong reflectors are observed beyond a 15 s delay time, in Figures 3 A4, B4, we show the traces with accurate spatial distribution for the first 15 s delay time only, which corresponds to a depth of ca. 100 km. This is also the depth range to which we migrate the reflectors (Figure 4). The limitation to the top 100 km is not inherent to the GloPSI-technique, but a consequence of our spatially sparse data set. We would like to point out that during the processing, no further normalization is applied, except for the whitening. Over the array, and over various depths, the relative amplitudes should be meaningful. However, imaging is done without loss correction and hence the amplitude distribution is somewhat biased in the direction of the middle and upper crust. In the frequency band we selected, there is a strong focus on the Moho and upper mantle. As a last step, the station-wise stacks are time-to-depth converted using the P-wave velocities of the same simplified 1D velocity model that was used by Singer et al. (2017a) for the RF depth migration. The model assumes a constant crustal Vp velocity gradient of 0.0118 s −1 starting from 5.9 km/s at the surface. The average Moho depth is set at around 50 km depth, with a Vp jump from 6.5 to 8 km/s. The mantle Vp gradient is then 0.0016 s −1 . These gradients are based on a local minimum 1D Vp velocity model obtained from regional earthquakes (Diehl et al., 2017). The resultant first order lithospheric structure across the orogenic wedge in Western and Eastern Bhutan is shown in Figure 4. A depth of 100 km is reached after 14.3 s with the 1D Vp model.

RESULTS AND INTERPRETATION
From the obtained final GloPSI images alone (Figure 4), it would be difficult to confidently interpret the depth and geometry of the reflectors. We recall that the methods' strength lies in laterally tracing a physical parameter produced by material differences in the crust that produces a similar signal at depth. We cannot resolve fine structures, such as thin layering, but we can resolve larger structures that show continuity in space which can give us important information about the seismic facies. Moreover, because of the large station spacing and the inability to image Frontiers in Earth Science | www.frontiersin.org FIGURE 3 | Data processing steps for the profiles in Western (A) and Eastern (B) Bhutan. After autocorrelation and stacking (step 1), the delta pulse is removed (step 2), and the multiples are suppressed (step 3). In step 4, the traces are shown with accurate spatial distribution for the first 15 s, which correspond to a depth of around 100 km. The amplitudes are normalized to the maximum value after each processing step, which causes the shift in colors after the delta pulse removal.
steep dips, the lateral continuity is not accurately reproduced in the images. This could be potentially improved with more data. From the reflection retrieval, shown in Figure 3, step 1, we see a clear change in reflectivity at around 10 s which is an indicator of the boundary between the crust (mostly positive amplitudes) and the mantle (mostly negative amplitudes), for both Eastern and Western Bhutan. For the interpretation of structures, we show our migrated results together with the location of converter points identified in the RF study by Singer et al. (2017a) (shown as filled circles) in Figure 4. Thanks to these earlier observations, we can point to the two most prominent reflectors and interpret them, due to known interfaces, as being the Moho and the MHT, and also point to some reflectivity patterns within the crust. In Western Bhutan, the Moho can clearly be identified at around 50 km depth and is slightly dipping toward the North, in agreement with the flexural deepening of the downgoing plate.
The shallowest clear horizon is seen at approximately 20 km depth, and is identified as the MHT, with a flat geometry in the East, and northward deepening in the West. In the East, the Moho reflector dips at a lower angle toward the North. The shallowest clear horizon interpreted as MHT is seen at approximately 15 km depth. Within the crust, there is a near parallel layering in the East, while such a pattern is not detected in the West. As a result of the comparison between the two techniques, with RF being especially sensitive to S-wave structure, and GloPSI being especially sensitive to P-wave structure, we can say: 1. The Moho is well visible in Eastern and Western Bhutan with both methods, with matching results. The two techniques corroborate each other for both the flat geometry of the Moho (southern parts of both profiles) and the slightly dipping parts (northern parts). The Moho kink is not so obvious from GloPSI; the Moho signal is indeed not visible between approximately 60 and 80 km profile distance (on the Western profile), while it is visible elsewhere. In a similar way, the Moho signal is fading between approximately 60 and 80 km distance along the Eastern profile.
2. The shallower reflectors (less than 20 km deep) are well visible in the GloPSI images, better than from the RF. These reflectors are likely linked to the MHT, the predominant intracrustal collision structure in the Himalaya orogen (Coutand et al., 2014). To the images (Figure 4, violet line), we added the MHT depth proposed by Coutand et al. (2014). The authors then collected multithermochronologic data at various points along similar profiles. In particular, in Eastern Bhutan, the depth and geometry of the MHT correlate well with our results. In Western Bhutan, we added the extrapolated reflectors obtained from the wide-angle seismic reflection profile INDEPTH Tib-1 ( Figure 4A, black arrows) that was acquired in the Higher Himalayas (Hauck et al., 1998), northwest of Bhutan (Figure 1).
In Western Bhutan, the dip of the MHT and the presence of anisotropy associated to it, are structural reasons for preventing a good resolution of that feature by GloPSI, and might cause the discrepancy between the location of the MHT from different studies (yellow and purple lines in Figure 4A display some consistent offset). Schulte-Pelkum et al. (2005) interpreted the RF signals from the MHT in eastern Nepal as an anisotropic shear zone representing the upper crustal detachment zone. Singer et al. (2017a) instead imaged the MHT as a negative RF arrival and a ramp structure.
3. Intra-crustal reflectivity in Eastern Bhutan is stronger than in Western Bhutan, hinting at a difference in physical properties or deformational history of the crust penetrating beneath the orogenic wedge (Singer et al., 2017a). Although the flexural shape of the downgoing plate's lithosphere seems to be similar between West and East Bhutan (Hammer et al., 2013;Hetenyi et al., 2016a), physical characteristics of the crust reflect spatially smaller scale properties, for which there is no clear information available, concerning their origins. We speculate that the reflectivity is inherited, but note that further research in the footsteps of e.g., Talwani et al. (2016) and Mallick et al. (2020) is required.

DISCUSSION
With the use of temporary seismic network data and as few as 59 teleseismic events, the GloPSI approach allowed us to outline the main lithospheric structures of the Bhutan Himalaya along two transects. The interpretation was guided by earlier RF observations, however the GloPSI results provided confirmation and useful results. The main conclusion is that the incidence of the Moho and MHT geometry in the two sets of images suggests that P and S waves both yield coherent information about the structure. A benefit of using GloPSI, is that it only exploits the vertical component of seismic records. Therefore, there is potential to use GloPSI on seismic networks that rely on one-component sensors, as long as these can record data in the required frequency band (1 s to at least 5 s). Larger aperture networks and a comparatively simple structural environment, complemented with permanent network data, likely increase the chances for GloPSI results to be interpreted alone, without RF or other prior data.
Among the drawbacks of GloPSI, one can note the (partial) loss of relative amplitude information during successive renormalization of waveforms. Moreover, the structure that can be imaged with GloPSI is more limited in space than with RF. With GloPSI, the main contribution comes from structure directly below the stations, whereas with RF, conversion points may be selected that are various distances away from the stations. Hence, for GloPSI, denser array deployments are needed in order to obtain laterally continuous images.
Orogenic regions tend to have strong lateral variability in elastic properties, and even if the incidence angles for GloPSI are closer to vertical than RF studies, the 3D effects should still be visible on the resulting images. The existence of accurate 3D velocity models and denser field arrays would overcome this limitation and allow for more advanced migration techniques. In most areas on Earth, the distribution of sources is more favorable for GloPSI than for RF methods to be applied. Both methods use transmission responses with different incidence angles below the stations, and different sensitivities to S and P-wave velocity structures, which results in imaging overlap, but also partly exclusive, subsurface information. Combining GloPSI and RF with their individual weaknesses and strengths allows, in general, one to interpret lithospheric structure with a variety of network configurations. With two, rather short but densely spaced profiles across the Himalayan orogenic wedge in Bhutan, we observe a matching pattern of Moho and MHT geometry, as well as east-west differences in shapes and reflectivity patterns.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
AO processed the data with codes from ER and wrote the initial version of the manuscript. IB and GH contributed to the interpretation and manuscript writing. All authors contributed to the article and approved the submitted version.

FUNDING
The deployment of the temporary seismic GANSSER network in Bhutan was funded by the Swiss National Science Foundation (SNF), Grant 200021-143467 (GH). The GANSSER acronym was chosen to commemorate the famous geologist Augusto Gansser . The instrumental pool of the SED and SEG groups at ETH Zurich provided the station hardware.