Early cretaceous bimodal volcanic rocks in Wuga Co area, central tibet: The first identification of direct products derived from slab sinking in the Bangong-Nujiang suture zone

Introduction: The Bangong-Nujiang Suture Zone (BNSZ) in central Tibet is a remnant of the Bangong-Nujiang Ocean that records its entire Wilson Cycle. The model of divergent double-sided subduction (DDS) is crucial for elucidating the evolution of tectonomagmatic activity on both sides of the BNSZ and for understanding why no high-pressure metamorphic rocks occur in the BNSZ. However, the DDS geodynamics remain poorly constrained. In particular, there is a lack of reports on magmatic rocks directly associated with slab sinking in the DDS terminal stage. Methods: This study presents new geochronological, geochemical, and isotopic data for the Early Cretaceous bimodal volcanic rocks around the Wuga Co area. Results: The bimodal volcanic rocks are divided into the Wuga Co rhyolites (SiO2 = 77.0–79.0 wt%) and the Wuga Co basaltic andesites (SiO2 = 53.9–55.5 wt%). The isotopic values of the Wuga Co basaltic andesites with low (87Sr/86Sr)i values (+0.7040 to +0.7044) and high εNd(t) values (+3.8 to +4.1) lie among three endmembers (the BNO sediments in accretionary wedge, depleted mantle and the BNO slab). Discussion: These values indicate the partial melting of a mantle peridotite that interacted with the subducted slab and sediment in the accretionary wedge, which was caused by the sinking of the Bangong-Nujiang oceanic slab. The Wuga Co rhyolites (108 Ma) have low (87Sr/86Sr)i values (+0.703 to +0.706), high εNd(t) values (+2.25 to +2.49), and high εHf(t) values (+5.6 to +10.0). These values indicate that the rhyolite formed by partial melting of juvenile basaltic crust. This study also collected Hf isotope data from both sides of the BNSZ to constrain its evolution. Our results show that the εHf(t) values of magma on both sides of the BNSZ were elevated simultaneously at 130 Ma, which may be caused by the Bangong-Nujiang oceanic slab rupture. Based on these new data, we propose that the Bangong-Nujiang oceanic slab ruptured from the two overlying terranes at approximately 130 Ma and subsequently sank into the mantle at approximately 108 Ma.


Introduction
Continent-continent hard collisions in single-sided subduction zones such as the Himalayan and Dabie-Sulu belts create suture zones (Luo A. B. et al., 2022a), resulting connection between the passive continental margins on one side and active continental margin on the other side (Yin and Harrison, 2000;Luo A-B. et al., 2022b). As oceanic lithosphere is consumed by subduction, collision leads to the thickening of the continental crust and the formation of high-pressure metamorphism, such as blue-schist facies. However, there are many suture zones globally where the blue-schist facies is not present, such as the Precambrian Jiangnan Orogen in South China (Zhao, 2015), the Early Paleozoic Lachlan fold belt of southeastern Australia (Soesoo et al., 1997), or the Late Paleozoic-Mesozoic Xing'an-Mongolia Suture Zone in the Central Asian Orogenic Belt (Pei et al., 2018). Soesoo et al. (1997) first elaborated the process of arc-arc soft collisional development caused by divergent double-sided subduction (DDS) to explain the Lachlan Belt in eastern Australia, where high-pressure metamorphism is absent. Therefore, different subduction models result in different collision forms. Single-sided subduction eventually leads to continent-continent hard collisions and DDS leads to arc-arc soft collisions. Analyzing the geochemical characteristics of these igneous rocks is essential to understand the subduction form and the collisional model it causes.
The Bangong-Nujiang Suture Zone (BNSZ) is a remnant of the Meso-Tethys Ocean and an important tectonic belt sandwiched between the southern Qiangtang Terrane and the Lhasa Terrane in the central Tibetan Plateau (Allégre et al., 1984;Girardeau et al., 1984;Yin and Harrison, 2000;Zhu et al., 2013). However, the subduction polarity of the Bangong-Nujiang Ocean (BNO) and the collisional model caused by it is still a matter of debate Guo et al., 2022;Liu et al., 2022;Zhang et al., 2022). Two main models are in dispute. The first is a hard collision caused by northward subduction (Han et al., 2020;Wang et al., 2020;Yan and Zhang, 2020) and the second is a soft collision caused by DDS Li S. M. et al., 2018;Li et al., 2019;Luo A. B. et al., 2022a;Luo A-B. et al., 2022b;Shi et al., 2022). The former is supported by the following: 1) A subduction complex (Mugagangri Formation) is distributed along the BNSZ and the southern Qiangtang Terrane (Hu W-L. et al., 2022). 2) continent-arc igneous rocks are distributed in the southern margin of the Qiangtang Terrane (Li S-M. et al., 2014;Tong et al., 2022). 3) The dip direction of thrust fault at the edge of the BNSZ is northward (Kapp et al., 2007). The latter is evidenced by the following points: 1) Early Cretaceous adakites and sodium-rich rocks occur in the northern margin of the Lhasa Terrane and the southern margin of the Qiangtang Terrane (Li X. K. et al., 2018;Liu et al., 2020). 2) High-resolution seismic reflection data reveals the Mesozoic DDS relict preserved within the crust. . 3) The comprehensive research on lava (High-K quartz-diorites, low-K tholeiitic andesites, bajaitic latites) in southern Qiangtang and sediments records the process from DDS to soft collision (Li et al., 2019).
To test the two disputed models, we selected bimodal volcanic rocks in the Wuga Co Lake for the following three reasons: 1) Felsic igneous rocks are easier to date than the mafic rocks; 2) The bimodal volcanic rocks are located in BNSZ, and such sample reports are scarce; and 3) Unveiling the petrogenesis of these samples has important implications for tracing the collisional orogenic processes between the northern Lhasa and southern Qiangtang terranes. Therefore, new geochemical and geochronological evidence for the bimodal volcanic rocks in the Wuga Co area is presented in this paper. These data, together with Hf isotopic data, allow us to investigate the evolution of the BNSZ.

Geological background and samples
The Tibetan-Himalayan Plateau is the highest plateau on Earth. It is a continuous accretion of the Songpan-Ganzi, Qiangtang, Lhasa, and Indian terranes. These terranes are separated by the Jinsha Suture Zone, Bangong-Nujiang Suture Zone, and Indus-Yarlung Zangbo Suture zones from north to south, respectively (Chang and Lo, 1973;Wang and Wright, 1987;Yin and Harrison, 2000) ( Figures  1A, B). Furthermore, the Lhasa Terrane is subdivided into northern, central, and southern terranes by the Shiquan River-Nam Tso Mélange Zone and the Luobadui-Milashan Fault. The Qiangtang Terrane is subdivided into northern and southern parts by the Longmu Co-Shuanghu Mélange Zone.
The Cretaceous volcanic rocks in the north-central Lhasa Terrane include the Zenong Group and the Duoni Formation, as well as the Qushenla Formation in the southern part of the BNSZ . The widespread Jurassic and Cretaceous magmatic activity in the north-central Lhasa Terrane is interpreted as the subduction of the BNO and continental collision (Li et al., 2015a;Hao et al., 2016a;Zhu et al., 2016;Yang and Wang, 2019). The Early Cretaceous Zenong Group contains high-K calc-alkaline felsic magma and pyroclastic rocks with an average thickness of 1000 m and an area of about 20000 km 2 (Chen et al., 2014). The eruption age of the Duoni Formation is 115-100 Ma, and the main lithology is basaltic andesite and rhyolite, with a thickness of 3000 m . The Qushenla Formation is sparsely distributed and consists of andesite and minor basaltic andesite .
Early studies in the sourthern Qiangtang Terrane found no evidences of arc magmatism. The lack of evidences for magmatic action in the BNSZ hinders its understanding (Allégre et al., 1984;Zhu et al., 2016). Later, identification of 185-170 Ma granitoids in the Amdo microplate, interpreted as a Jurassic continental arc along the sourthern Qiangtang Terrane, provides evidences of northward subduction (Guynn et al., 2006). Recent studies show that many Jurassic to early  arc-type igneous rocks are located within the southern Qiangtang Terrane. The lithology of these magmatic rocks consists of gabbroic diorites, diorites, granodiorites, and granites (Kapp et al., 2005;Ran et al., 2015).
The BNSZ spans more than 2,000 km in central Tibet, and is characterized by flysch, mélange, and ophiolite fragments, representing remnants of the BNO . In the BNSZ, the complete ophiolite sequence has been found in the areas of Ritu, Dongcuo, Dongqiao, Amdo, and Dingqing (Fan et al., 2014;Liu et al., 2014;Zhu et al., 2016). The Early Cretaceous Shamuluo Formation is unconformable on ophiolite rocks and flysch. The Jurassic Mugagangri Formation flysch consists of deformed deep-sea turbidites and ophiolite remnants. The lithology of the Shamuluo Formation is mainly sandstone and limestone deposited in a shallow marine environment .
The Qushenla Formation rocks have Wuga Co basaltic andesites (WGBA) and Wuga Co rhyolites (WGR) sandwiched within WGBA and have an extent of approximately 530 km 2 , located approximately 52 km northwest of Nima County in the Wuga Co area (Figures 1B,C). The Qushenla Formation unconformably contacts the Jingzhushan Formation and contacts the Late Triassic Wuga Formation along a fault ( Figure 1C).
The basaltic andesite is gray black, massive, and porphyritic (Figures 2A, B). The basaltic andesite contains plagioclase phenocrysts and a microcrystalline matrix ( Figure 2C). Plagioclase grains are euhedral and elongate and show polysynthetic twinning. These phenocrysts are about 0.5 mm long and 0.1 mm wide. Thin section observations indicate alteration. The occurrence of chlorite were altered from plagioclase and the carbonation were altered from matrix. The WGR is grayish brown and rhyotaxitic ( Figure 2D). The phenocrysts are mainly quartz and have a preferred orientation. The length of quartz grains is 0.2-0.4 mm, and the width is about 0.05 mm. No alteration is found under microscope. Detailed operating conditions for the laser ablation system, the ICP-MS instrument, and the data reduction are given in Zong et al. (2017). Laser sampling was performed using a GeolasPro laser ablation system that consists of a COMPexPro 102 ArF excimer laser (wavelength of 193 nm and maximum energy of 200 mJ) and a MicroLas optical system. An Agilent 7900 ICP-MS instrument was used to acquire ion-signal intensities. Helium was applied as a carrier gas. Argon was used as the make-up gas and mixed with the carrier gas via a T-connector before entering the ICP. A wire signal-smoothing device is included in this laser ablation system (Hu et al., 2015). The spot size and frequency of the laser were set to 30 µm and 6 Hz, respectively. Zircon 91500 and glass NIST610 were used as external standards for U-Pb dating and trace-element calibration, respectively. Each analysis incorporated a background acquisition of approximately 20-30 s followed by 50 s of data acquisition from the sample. An Excel-based software ICPMSDataCal was used to perform off-line selection and integration of the background and analyzed signals, time-drift correction, and quantitative calibration for trace-element analysis and U-Pb dating (Liu et al., 2008;Liu Y. et al., 2010). Concordia diagrams and weighted mean calculations were made using Isoplot/ Ex_ver3 (Ludwig, 2003).
In situ zircon Hf isotope analyses were performed on 36 previously dated zircon grains from the WGR samples. The points used for the Hf analyses were the same as those used during the LA-ICP-MS analyses, which were selected from cathodoluminescence (CL) images. The data were collected using a NEPTUNE Plus MC-ICP-MS. A single spot ablation mode with a spot size of 55 μm was used to acquire the data. Each measurement consisted of 20 s of background signal acquisition followed by 50 s of ablation signal acquisition. Detailed information on the operating conditions and analytical methods can be found in Hu et al. (2012). The analyzed 176 Hf/ 177 Hf ratios for the zircon standard (91500) were 0.282299 ± 31 (2σ n, n = 40), which are similar to the 176 Hf/ 177 Hf ratios of 0.282302 ± 8 and 0.282306 ± 8 (2σ) for the standard when determined by the solution method (Goolaerts et al., 2004;Jon et al., 2004). Off-line selection, integration of analyte signals, and mass bias calibrations were performed using the ICPMSDataCal program (Liu Y. S. et al., 2010b).

Whole-rock Nd isotopic compositions
Nd isotope analyses were performed on a Neptune Plus MC-ICP-MS (Thermo Fisher Scientific, Dreieich, Germany) at the Wuhan SampleSolution Analytical Technology Co., Ltd. (Wuhan, China). The exponential law, which initially was developed for thermal ionization mass spectrometry measurement (Russell et al., 1978) and remains the most widely accepted and used with MC-ICP-MS, was used to assess the instrumental mass discrimination in this study. The

Whole-rock major and trace elements
Whole rock elements compositions of samples were analyzed at Wuhan SampleSolution Analytical Technology Co., Ltd., Wuhan, China. Major element abundances were determined using Zsx Primus II wavelength dispersive X-ray fluorescence spectrometer (XRF). The agate crushed sample was dried first at 106°C for 2-4 h. After it, 25 g of the annealed sample was mixed with 5 g of sodium tetraborate in platinum crucible. This mixture was melted in high frequency furnace. After cooling, the melt was diluted to 100 ml by adding water. All major element analysis lines are kα and the standard curve uses the national standard materials, rock standard sample: GBW07101-14 The data were corrected by theoretical α coefficient method. The relative standard deviation (RSD) is less than 2%.
Trace element analysis of whole rock were conducted on Agilent 7700e ICP-MS at the Institute of Mineral Resources, Chinese Academy of Geological Sciences. The detailed sample-digesting procedure was as follows: 1) Sample powder (200 mesh) were placed in an oven at 105°C for drying of 12 h; 2) 50 mg sample powder was accurately weighed and placed in a Teflon bomb; 3) 1 ml HNO 3 and 1 ml HF were slowly added into the Teflon bomb; 4) Teflon bomb was putted in a stainless steel pressure jacket and heated to 190°C in an oven for >24 h; 5) After cooling, the Teflon bomb was opened and placed on a hotplate at 140°C and evaporated to incipient dryness, and then 1 ml HNO 3 was added and evaporated to dryness again; 6) 1 ml of HNO 3 , 1 ml of MQ water and 1 ml internal standard solution of 1ppm In were added, and the Teflon bomb was resealed and placed in the oven at 190°C for >12 h; 7) The final solution was transferred to a polyethylene bottle and diluted to 100 g by the addition of 2% HNO 3 .

In situ zircon Lu-Hf isotopes
Thirty-nine Hf isotopic analyses were carried out on 39 zircons separated from the WGR (nb-1). The initial 176 Hf/ 177 Hf ratios of the WGR samples vary from 0.282864 to 0.282987. The zircon εHf(t) values range from +5.6 to +10.0 and hafnium model ages range from 368.4 Ma to 539.9 Ma. The zircon Lu-Hf isotope data are listed in Supplementary Table S2.
Significant FC for these basaltic andesites can be ruled out, for the following four reasons: 1) An obvious linear relationship is not shown in the Harker diagram (Supplementary Figure S1). 2) Basaltic andesites produced by FC require basaltic magma as the parent magma; however, no contemporaneous basaltic magma has emerged around the study area . 3) The FC of plagioclase is often accompanied with negative Eu anomalies, which is not shown in our samples ( Figure 3C). 4) The WGBA samples show no FC trend in the La/Sm versus La diagrams or La/Yb versus La diagrams ( Figures 4A, B).
The mixture of mafic magma and felsic magma can produce basaltic andesite magma. However, the signature of homogeneous isotopic values ( 87 Sr/ 86 Sr = 0.704244-0.704494, 143 Nd/ 144 Nd = 0.512797-0.512806) suggests that no magma mixing occurred . Moreover, neither enclaves nor disequilibrium textures were observed in the field or microscope observations ( Figure 2C; Perugini and Poli, 2012;Liu et al., 2019). Further support for this conclusion is given by the consistent chondrite-normalized REE and primitive-mantle-normalized trace-element patterns (Figures 3C, D;Liu et al., 2019).
In the process of partial melting and fractional crystallization, the distribution coefficient of strongly incompatible elements relative to compatible elements can be ignored. According to this, Allègre and Minster, (1978) has carried out graphical and theoretical derivation for partial melting and fractional crystallization. When the concentration ratio of strongly incompatible elements to compatible elements is plotted against the concentration of strongly incompatible elements, the trajectory of partial melting is an inclined line, while the separation crystallization is a horizontal line. The values of WGBA have a strong positive correlation in the La/Yb versus La diagrams and that  Figures 4A, B), which is consistent with partial melting trends. Therefore, we suggest that the WGBA samples were generated by partial melting of their source rocks. Therefore, their source characteristics can be represented by their geochemical features.
Partial melting and fractional crystallization cannot change the isotope ratios. The distinct differences of the εNd(t) values between the WGBA (+3.8 to +4.0) and the WGR (+2.2 to +2.5) rule out the FC model. Furthermore, FC trends of the WGR were not observed in the Harker diagram (Supplementary Figure S1) and the low Mg # values (13.5-18.2) of the WGR are consistent with crustal melts (Mg # <40) (Rapp and Bruce Watson, 1995;Wei et al., 2017). All these geochemical characteristics suggest that the WGR resulted from crustal anatexis.

The WGBA
We suggest that the WGBA are derived from mantle peridotite because they have low 87 Sr/ 86 Sr values (0.703993-0.704366) and high ε(Nd) values (3.84-4.07). Nb, Ta, and Hf elements have the same geochemical behavior because of their similar incompatibility during magmatic evolution. The values of Nb/Ta and Hf/Ta can reflect the features of the source, because they do not change during partial melting or FC (Pearce and Peate, 1995;Wang et al., 2021;Zhao et al., 2021). The Nb/Ta (16.67-17.62) and Hf/Ta (5.82-6.18) values of our samples are similar to those of the primitive mantle (Nb/Ta = 17.39, Hf/Ta = 7.54), providing further evidences for their derivation from mantle peridotite. However, in the Th/Yb versus Nb/Yb diagram, the WGBA samples are distributed above the MORB-OIB array ( Figure 5). Therefore, we suggest that there are other endmembers involved in the source region of the WGBA (Pearce, 2008;Huang et al., 2017;Wei et al., 2017;Xu et al., 2017;Zhong et al., 2017;Tang et al., 2018).
In the Ni versus Cr diagram, the WGBA samples plot on the mixing line between slab melt and mantle peridotite ( Figure 6A), similar to the Biluoco magnesian andesite which is derived from the interaction between mantle peridotite and slab-derived materials  Figure 3B.

FIGURE 5
Nb/Yb vs. Th/Yb diagram (Pearce, 2008) of the Wuga Co basaltic andesites (WGBA) from the Wuga Co area. S, C, W, and F refer to enrichment by subduction environment, continental contamination, enrichment within plate, and fractional crystallization, respectively. Data and symbols are the same as in Figure 3B.

Frontiers in Earth Science
frontiersin.org (He et al., 2018). In addition, The WGBA also show the characteristics of low Y and heavy rare Earth elements (HREE; Figure 3C), which are generated by garnet residue in the source. These features, combined with high Na 2 O/K 2 O values (3.6-7.0), can be interpreted as the involvement of oceanic crust in the source region of the WGBA (He et al., 2018). In summary, the oceanic slab is part of the WGBA source materials. The WGBA samples show not only the geochemical characteristics of oceanic slab but also the geochemical characteristics of sedimentary rocks. The ratios of Th/Yb in volcanic rocks can be used to identify the involvement of sediments and/or sediment-derived melts in their generation (Woodhead et al., 2001). On plots of Th/Yb versus Th/Sm ( Figure 6B), the WGBA samples lie on the array between MORB and marine sediments, which indicates an obvious contribution from sediments in their origin (Li et al., 2015a). The sedimentary rock characteristics are also shown by the Ba/Th versus Th diagram ( Figure 6C).
Furthermore, we suggest that the endmember of the sediments in the WGBA source area is likely from the accretionary wedge because of the enrichment in large-ion lithophile elements (LILE) and depletion in HFSE ( Figure 3D). Previous studies have shown that there are three reasons for LILE enrichment and HFSE depletion: 1) Crustal contamination (Sun and McDonough, 1989;Li et al., 2017); 2) Involvement of aqueous fluids under a continental arc (Hawkesworth et al., 1993;Pearce and Parkinson, 1993;Chai et al., 2015); and 3) Involvement of mélange rocks (Marschall and Schumacher, 2012;Hao et al., 2016b). It has been explained in 5.1.1 that crustal contamination did not occur during petrogenesis. The involvement of aqueous fluids has been excluded by the Ba/Th versus Th diagram ( Figure 6C) and (Hf/Sm) PM versus (Ta/La) PM diagram ( Figure 6D) because the WGBA samples lie out the trend of aqueous fluids. Therefore, the conclusion that the mélange participates in the source rocks of the WGBA can be supported by the geochemical features of HFSE depletion and LILE enrichment. All these geochemical features indicate that sedimentary rocks in the accretionary wedge are involved in the source.
In summary, the geochemical data suggest that the WGBA is likely derived from partial melting of a mantle peridotite that interacted with subducted slab melts and melts from sediments in the accretionary wedge.  (Perez et al., 2018); (D) (Hf/Sm) N versus (Ta/La) N diagram (La Flèche et al., 1998). Literature data about the Rutog andesite and Oman andesite are from Zhao et al. (2021) and Kanayama et al. (2012), respectively. Related data about subducting sediments is from Zhao et al. (2020). Data and symbols are the same as in Figure 3B.
Frontiers in Earth Science frontiersin.org

The WGR
The zircon in which Hf isotopes have a high closure temperature has a very low 176 Lu/ 177 Hf ratio. There is no significant radiogenic Hf accumulation after formation (Wu et al., 2007). Therefore, the 176 Hf/ 177 Hf ratio in zircon can represent the source characteristics. The zircons in the WGR have high ε(Hf) values (+5.6 to +10.0) and young two-stage model Hf ages (TDM2) values (516.6 Ma to 788.0 Ma), implying that they are derived from a juvenile basaltic crust.
In summary, the major and trace element geochemistry, together with the Nd-Hf isotopic compositions, suggest that the WGR is the product of partial melting of juvenile basaltic crust.

Isotopic evidences
To test that the source of WGBA is involved in the BNO sediments in accretionary wedge, depleted mantle and the BNO slab, it is necessary to determine the geochemical characteristics of the three endmembers.
A mélange is a disordered rock mixture formed by a combination of sedimentary rocks scraped from a subducted oceanic plate, some oceanic crust fragments, continental fragments, in situ sediments, and turbidite (Li C. et al., 2020). In this definition, the geochemical features of sediments in accretionary complex can be represented by the geochemical features of subducting sediments. In existing research, the assumed subducting sediments in the BNSZ, as represented by the Yanshiping Group sandstone from Duobuza (Figure 7), can be interpreted as sediments in accretionary complex Li J-X. et al., 2016).
The ophiolite has been interpreted as the remnants of the oceanic slab. The geochemical characteristics of ophiolite in BNSZ reported by previous studies can represent BNO slab . On the other hand, the characteristics of depleted mantle under BNSZ have not been reported, so we quoted the global depleted mantle values (Figure 7).
The isotopic characteristics of WGBA that plotted among the BNO sediments in accretionary wedge, depleted mantle and the BNO slab (Figure 7) support that the WGBA samples are derived from partial melting of a mantle peridotite that interacted with subducted slab melts and melts from sediments in the accretionary wedge.

Divergent double subduction and oceanic crust rupture
There are two main viewpoints in dispute about the subduction polarity of the BNO. Some studies have shown that the BNO subducted northward under the southern Qiangtang Terrane, and others support the conclusion that the BNO subducted not only northward beneath the southern Qiangtang Terrane but also southward beneath the northern Lhasa Terrane Ma et al., 2020;Peng et al., 2020;Jiang et al., 2021). Recent studies have shown that arc-related magmatic rocks are found on both sides of the BNSZ . We consider that the BNO underwent DDS beneath the southern Qiangtang and the northern Lhasa terranes.
Because crust-mantle interaction is inevitable during mantle melting, the positive εHf(t) for the magmatic rocks in the Lhasa Terrane is likely underestimated, and should be higher (Zhu, 2011;Li Y. J. et al., 2020;Cao et al., 2022). Hence, the maximum value of εHf(t) collected is used as the participation degree of the mantle, which is reasonable. In a classical DDS system, the rupture from the overlying crust occurs at the initial stage of oceanic crust detachment. At this stage, mantle materials will fill the gap caused by the rupture, which would lead to the more mantle participate in the source of related magmatic rocks (Soesoo et al., 1997;Zhao, 2015;Zhu et al., 2016;Zhang et al., 2017). The maximum of collected ε(Hf) values in both the southern Qiangtang Terrane and the northern-central Lhasa Terrane increase at 130 Ma (Figure 8), simultaneously. Therefore, we suggest that the rupture of the BNO from both sides occurs at 130 Ma. This age has been endorsed by other studies Li Y. et al., 2018;Li et al., 2019). Soesoo et al. (1997) proposed that the magmatism caused by oceanic lithosphere completely sinking into the mantle has two characteristics. First, the related magmatic rocks are the products of three endmembers: sinking oceanic crust, mantle peridotite, and accretionary wedge sediments. Second, the related bimodal volcanic rocks are distributed in a hinge zone. The source of our sample similarly has three components: the BNO sediments in accretionary wedge, depleted mantle and the BNO slab (detailed explanation has Plots of ( 87 Sr/ 86 Sr) i versus εNd(t). Related data about the Baingoin A-type granites and Langla granites are from Hu et al. (2019) and references therein. The values for subducting sediments in the Bangong-Nujiang Suture Zone (BNSZ), depleted mantle, and Banggong-Nujiang oceanic slab are from Rehkamper and Hofmann, 1997;Wang et al., 2016), respectively. The εNd(t) and ( 87 Sr/ 86 Sr) i values of these data have been modified to 108 Ma. Data and symbols are the same as in Figure 3B.

Frontiers in Earth Science
frontiersin.org been given in Section 5.2.1). In addition, the WGBA (SiO 2 = 53.9-55.5wt%) and interbedded WGR (SiO 2 = 77.0%-79.0wt%) constitute a bimodal volcanic suite. In summary, our sample is the product of the BNO slab sinking into the mantle, which occurred at approximately 108 Ma.

Geodynamic model
The sustained DDS eventually leads to the closure of the BNO . The inverted U-shaped oceanic plate locked the flow of the asthenosphere. Therefore, magmatism was weak during this period Soesoo et al., 1997). However, closure of the oceanic basin leads to a locked but unstable situation. The continued sinking of the cold and dense BNO slab eventually led to its rupture from the two overlying terranes (southern Qiangtang and northern Lhasa terranes) at approximately 130 Ma (Figure 9). Subsequently, the rupture propagates from both sides to the hinge zone, resulting in the BNO slab sinking into the mantle at 108 Ma. Mantle materials fill the gap between the sinking BNO slab and the BNSZ mélange, which results in the interaction of mantle peridotite, BNO slab, and accretionary wedge sedimentary rocks. The interaction led to the generation of the WGBA. The partial melting of oceanic crust fragments in hybrid rocks resulted in the generation of the WGR, which was attributed to continuous heating by the mantle flow.
The "slab sinking modal" in DDS have been first elaborated by Soesoo et al. (1997) to explain the Lachlan Belt in eastern Australia. This model has been applied in many other orogenic belts, such as the Precambrian Jiangnan Orogen in South China (Zhao, 2015), the Late Paleozoic-Mesozoic Xing'an-Mongolia Suture Zone in the Central Asian Orogenic Belt (Li S. M. et al., 2018;Pei et al., 2018). We take WGBA as an example to verify the rationality of  the "slab sinking model" in BNSZ. Furthermore, we make the first identification of direct products derived from slab sinking in the Bangong-Nujiang Suture Zone.

Conclusion
(1) In the Wuga Co area, the bimodal volcanic rocks of the Qushenla Formation are distributed within the BNSZ and were active at approximately 108 Ma. Sr-Nd isotope data of the samples plot among sedimentary rocks in mélange, mantle peridotite, and the BNO slab, which suggests that WGBA is the direct product of slab sinking.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

Author contributions
WC: Investigation and writing; ZL: Investigation and editing and polishing; NW: Investigation and editing GW: Editing and polishing; MZ: Investigation; NH: English polishing; YH: English polishing; XY: Investigation.

Funding
Geological survey of mineral resources in key prospective areas of strategic mineral resources (DD20221684).