Haida Gwaii (British Columbia, Canada): a Phanerozoic analogue of a subduction-unrelated Archean greenstone belt

Understanding the formation and evolution of Precambrian greenstone belts is hampered by gaps in the rock record and the uncertainty of the tectonic regime that was operating at the time. Thus identifying a modern analogue of a Precambrian greenstone belt can be problematic. In this paper we present geological, geochemical and petrological evidence outlining the case for Haida Gwaii (British Columbia, Canada) as a modern example of a greenstone belt. Haida Gwaii is comprised of two rift-related volcano-sedimentary sequences. The older (Early Triassic) Karmutsen volcanic sequence consists of subaqueous ultramafic-mafic volcanic rocks that are capped by marine carbonate and siliciclastic rocks. The younger (Paleogene) Masset bimodal volcanic sequence consists of tholeiitic and calc-alkaline basalt along with calc-alkaline silicic volcanic and intrusive rocks that are capped by epiclastic sandstones. The Karmutsen and Masset volcanic rocks have indistinguishable Sr-Nd-Pb isotopes demonstrating they were derived from a similar mantle source. Some of the Masset calc-alkaline rocks are compositionally similar to magnesian andesites (SiO2 = 56–64 wt%; Mg# = 0.50–0.64) that are typical of subduction-related Archean greenstone belts. We show that the calc-alkaline signature observed in the bimodal sequence of the Masset Formation is likely due to fractional crystallization of a tholeiitic parental magma under relatively oxidizing (ΔFMQ + 0.7) conditions indicating that a calc-alkaline signature is not prima facie evidence of a subduction setting. Given the geological and geochemical evidence, Haida Gwaii represents one of the best analogues of a modern subduction-unrelated Archean greenstone belt.

Beyond the volcanic nature of Archean greenstone belts their origin remains one of the most debated issues in the solid earth sciences [1][2][3] . It is unclear what precisely a greenstone belt is and whether all are representative of the same geological feature 4 . For example, some greenstone belts are interpreted to be ophiolites or slivers of oceanic crust whereas others are thought to be oceanic plateaux, volcanic rift complexes, or island-arcs [4][5][6][7][8][9][10][11][12][13] . In other words, some greenstone belts may be subduction-related whereas others are subduction-unrelated 14,15 .
The principal contentious issue is whether Archean greenstone belts were formed by plate tectonic (subduction, sea-floor spreading, continental drift) processes that are similar to modern Earth or if they are the consequence of convective-force (mantle plumes, diapirism, foundering) or 'vertical' tectonics [16][17][18][19][20][21][22][23][24][25] . In addition to their belt-like nature and accretionary-like structure, amongst the most compelling evidence in favour of a subduction setting for some greenstone belts is the presence of boninitic rocks within the mafic-ultramafic volcanic series and magnesian andesite and calc-alkaline silicic rocks in the bimodal series [26][27][28][29][30] . However, the apparent absence of Atlantic-style passive margin sequences, deep-water sedimentary rocks, blueschists, ultra-high pressure rocks, paired metamorphic belts, and the extent of thrust-related structures imply that the thermal and tectonic processes that were operating during the Archean were different than today 3,[31][32][33][34] . Consequently, finding a modern analogue of an Archean greenstone belt is important so that it can help to address the knowledge gap in their structural, petrological and temporal evolution.
The islands of Haida Gwaii located ~70 km west from mainland British Columbia (Canada) are comprised of two rift-related volcanic series which are separated by and intercalated with sedimentary sequences of marine The volcanic sequences of an archetypical 1,2 greenstone belt are comprised of subaqueous and subaerial volcanic and volcaniclastic rocks whereas the upper unit is predominately comprised of sedimentary rocks (Fig. 2). The lower portion of the volcanic unit consists of subaqueous ultramafic (komatiite) to mafic (tholeiite, boninite) volcanic rocks with minor felsic tuffs. The upper portions of the lower unit consist of a bimodal sequence of tholeiitic flows and calc-alkaline silicic (andesite to rhyolite) volcanic rocks 1,49,50 . In some cases the mafic-ultramafic sequences are separated by thin layers of calc-alkaline rocks at time intervals ranging from 3 to 30 million years 33 . Furthermore, sedimentary depositional gaps also exist between volcanic episodes that range in duration from 2 to

stratigraphy and Geological structure of Haida Gwaii
The stratigraphy of Haida Gwaii matches that of the volcanic sequence of Archean greenstone belts (Fig. 2) 55,56 . The Upper Triassic Karmutsen Formation is the basal unit of Haida Gwaii and is comprised of lower subaqueous mafic tholeiitic rocks with subordinate picritic rocks and an upper subaerial mafic tholeiitic unit that is capped by limestones, calcareous siltstone, tuff and well bedded, fine grained detrital rocks 38,39,55 . Fossils constrain the age of the basalt between middle Ladinian (~230 Ma) and Carnian-Norian (~224 Ma) 57 . On Vancouver Island, the Karmutsen tholeiites are mainly submarine and 6.1-6.6 km thick. The formation consists of ca. 2900 m of basal submarine pillow lavas overlain by ~600-1100 metres of pillow breccia and aquagene tuff, followed upwards by ~2600 m of massive basalt flows and sills interbedded with shallow water and subaerial sedimentary rocks 38,58,59 . Dykes and sills are locally present however, sheeted dyke centres have not been found. The termination of volcanism was followed by submergence leading to the deposition of Late Triassic platformal carbonates.
The Karmutsen Formation is overlain by a Middle Triassic sedimentary sequence of argillite, siltstone and bivalve-bearing limestone 41 . Between the uppermost Karmutsen Formation and lowermost Masset Formation are sequences of Lower Jurassic to Upper Cretaceous conglomerates, feldspatholithic wackes, sandstones, and shales. The San Christoval and Burnaby Island plutonic suites were emplaced during the middle to late Jurassic just prior to accretion 60 . The Karmutsen plateau (Wrangellia) was accreted to the North American margin via east-ward dipping subduction 61 . During the Late Cretaceous to Early Paleogene the region became volcanically active as subaqueous volcanic debris, breccias and flows were deposited on Upper Cretaceous conglomerates, sandstones and mudstones of the Honna Formation 55,56 .
Volcanic rocks of the Early to Middle Paleogene Masset Formation occur throughout Haida Gwaii and along the western coast of British Columbia. They cover an area of ~5000 km 2 in a continental marginal basin (half-graben) that is ~150 km wide and about 500 km long. The formation has also been intersected in offshore wells and probably 62 covers an area of ~70 000 km 2 . In synvolcanic extensional grabens and down-faulted blocks, the volcanic sequence is up to 2 km thick but the geophysical data as well as the offshore wells indicate a thickness of ~3.5 km 63 . The volcanism is related to the development of an extensional transform pull-apart basin that experienced syntectonic Eocene extension of 2% to 15% 62,64 . The extension during Middle Tertiary led to significant crustal thinning 65 and heating 57,66,67 . It is possible that the thermal regime during the Eocene was elevated due to the subduction of the Farallon plate ridge further to the south (~250 km from southernmost tip of Haida Gwaii) as transtensional stress opened the Queen Charlotte Basin 36,62,64 . Subduction-related mantle melting was restricted to the wedge beneath the North American margin and did not affect the lithospheric or sublithospheric mantle beneath Haida Gwaii 36 . The magmatic activities of the Masset Formation lasted from about 46.2 Ma to 11 Ma with a peak between 25 and 20 Ma 62,64,68,69 . The formation is made up of lava flows with minor pyroclastic and subvolcanic (dykes and sills) rocks. Volcanic rocks include both subaerial and subaqueous varieties. Based on the land exposure, it was estimated that the proportion of subaerial to subaqueous eruptive rocks is about 60:40 35 . The Masset rocks are distinctly bimodal with a predominance of mafic over felsic types with a ratio of about 2:1 35,37 .
There are two types of mafic rocks within the Masset Formation: aphyric to sparsely porphyritic tholeiitic basalts and strongly porphyritic calc-alkaline basaltic andesites and rare basalts 35,36 . The tholeiitic basalts contain sparse phenocrysts of olivine or olivine-plagioclase-clinopyroxene-oxides. Calc-alkaline types include phenocrysts of Ca-poor pyroxene (hypersthene), plagioclase and clinopyroxene (augite). There are also occurrences of relict aluminous amphibole with rims of orthopyroxene 35 . Silicic rocks are typically younger than the mafic types. The volcanic sequences were also intruded by plugs and stocks of approximately coeval (46-27 Ma) 60 monzodiorites, diorites and granodiorites (SiO 2 = 55-70 wt%) 60,68 . These intrusive rocks form two belts extending along the coast of Haida Gwaii. The Masset Formation is overlain and partially intercalated with sedimentary rocks containing both continental and marine facies.

petrogenesis of the Volcanic Rock series
Karmutsen Formation. The Karmutsen basalts typically have SiO 2 within a narrow range of 49-52 wt%, the Mg# (Mg 2+ /Mg 2+ + Fe 2+ tot ) mostly between 0.6 and 0.4, and show tholeiitic fractionation trends such as an increase of Fe, Ti and V with increasing differentiation. In addition to the basalts with 6-8 wt% MgO, the formation also contains picrites with high MgO (9 wt% to 20 wt%) contents 38,58 . The tholeiitic basalts have chondrite-normalized patterns enriched in light REE with (La/Yb) N ratios ~2-2.5. On the other hand, the picrites have (La/Yb) N equal to ~0.5 38 . The shapes of mantle-normalized patterns of most tholeiitic basalts show a decrease from Nb-Ta to more compatible elements including the heavy REE (HREE) and Y. Such profiles are common among OIB and E-MORB. In contrast, the picrites have depleted-LREE patterns (Fig. 3) 38 . Both picrites and tholeiitic basalts have high positive but overlapping initial ε Nd (t) values (ε Nd (t) = +6 to +8) indicating a derivation from a common asthenospheric source 38,70 . The isotopic composition of the Karmutsen volcanic rocks is similar to a Pacific mantle-plume source (e.g. Ontong Java and Caribbean oceanic plateau) with whole rock compositions and an overall structure of an oceanic plateau (Fig. 4) 38 . The picrites were likely derived from a depleted spinel lherzolitic source by a high degree of partial melting 38,70 . In comparison, the parental magmas of the Masset Formation tholeiitic basalt are also considered to be derived from a spinel peridotite source and that their observed compositional variations are consistent with low-pressure fractional crystallization 35,36 .
The primary liquid compositions and mantle potential temperatures (T P ) can be calculated from volcanic rocks providing they are representative of a liquid (i.e. not cumulate) and have only experienced www.nature.com/scientificreports www.nature.com/scientificreports/ olivine fractionation 71,72 . The calculated primary liquid compositions and mantle potential temperatures of the Karmutsen volcanic rocks of Haida Gwaii and Vancouver Island are summarized in Table 1. The primary liquids are picritic to komatiitic and the T P are estimated to be between 1395 °C and 1670 °C with the higher temperatures (1570 °C to 1670 °C) indicative of an anomalously hot regime that is within the expected range of a mantle plume 71 38,117 . Data are from the Karmutsen mafic rocks whereas the Masset data include both mafic and felsic rocks [35][36][37][38]75,76 . NHRL is the Northern Hemisphere reference line 118 while S & K denotes the reference growth curve 119 . www.nature.com/scientificreports www.nature.com/scientificreports/ the Isua greenstone belt of West Greenland 29,73 . The identification of magnesian andesites and boninitic rocks is used to interpret a subduction zone setting for the Wawa and Isua greenstone belts but in the case of Haida Gwaii the rocks are definitively rift-related as their eruption was contemporaneous with the opening of the Queen Charlotte basin [62][63][64][65][66][67] .
The N-MORB-normalized trace element patterns of the intermediate rocks are broadly similar to the tholeiitic rocks but there is less variability in the Rb, Ba, and Th concentration but more variability with respect to Nb-Ta, Sr and Ti (Fig. 3) 35 . The patterns of the tholeiitic rocks resemble those of E-MORB, OIB, and many 'high-Ti' tholeiites of continental flood basalt provinces. Radiogenic Sr-Nd-Pb isotopes of the mafic rocks plot along the mantle array between the fields for MORB and OIB (Fig. 4). The tholeiitic and calc-alkaline mafic suites are isotopically indistinguishable. It is inferred that both suites were derived from similar parental tholeiitic magmas by polybaric fractional crystallization of different proportions of plagioclase, mafic minerals, and Fe-Ti oxide minerals 35,36 .
The highly silicic rocks are plagioclase-phyric, two pyroxene-bearing, mainly peraluminous types. They have SiO 2 ranging from 65 to 77 wt% (LOI-free basis). The N-MORB-normalized patterns show light REE enrichment with flat heavy patterns with depletions of Ti, and Sr, and variable Nb-Ta (Fig. 3). Furthermore, the normalized patterns of the felsic rocks are enriched in large ion lithophile elements. As expected, these rocks are more enriched in highly and moderately incompatible elements than the basalts. In contrast to the basalts, most rhyolites show depletion in Ti and Sr, which suggests fractionation of plagioclase and Fe-Ti oxide minerals (ilmenite and Ti-magnetite). The felsic rocks, unlike the mafic rocks, display variable to negative Nb-Ta normalized anomalies. The variability of Nb-Ta could be related to the effects of crustal contamination or influence from a subduction modified mantle source but it is more likely that ilmenite fractionation played the primary role. The partition coefficients (D) of Nb and Ta for ilmenite in intermediate liquids is quite high (D Nb = 2.3-4.6; D Ta = 2.7-6.6) 74 . Consequently, as magma becomes more silicic both Nb and Ta will be depleted at a faster rate thereby generating negative normalized anomalies. Moreover, the felsic rocks have Sr-Nd-Pb isotopic compositions overlapping those of basalts including high positive ε Nd (t) values (≥+ 6). In addition, the ε Nd (t) values of the Masset Formation are similar to the Triassic Karmutsen tholeiitic basalts 38,70,75,76 . It is suggested that the mafic rocks of the Masset Formation were generated by melting of a spinel peridotite asthenospheric mantle source 36 . the Nature of the Calc-alkaline signature. Calc-alkaline rocks are a subgroup of the subalkaline series that are defined by the relations between the alkali metals (Na 2 O + K 2 O), lime (CaO), the timing of Fe-Ti oxide crystallization, and water content of the magma system [77][78][79][80][81][82] . Although the term is widely used, there remain inconsistencies in its application. The calc-alkaline 'signature' is produced primarily due to water concentration, which affects the timing of plagioclase crystallization relative to ferromagnesian minerals, and relative oxidation state (timing of Fe-Ti oxide crystallization) but has little to do with CaO and alkalis 82 . Consequently, calc-alkaline rocks are commonly associated with volcanic-arc settings probably due to their greater likelihood of being water saturated coupled with auto-oxidation during degassing and crystallization in the crust 79,83,84 . However, calc-alkaline rocks are not restricted to arc settings and are also found at within-plate settings 81 . www.nature.com/scientificreports www.nature.com/scientificreports/ The calc-alkaline silicic rocks in the Masset Formation are likely derived by fractional crystallization of a tholeiitic parental magma under relatively oxidizing conditions 37 . Here we show models to demonstrate the effects of relative oxidation state (ΔFMQ −2 to +2) on the chemical evolution of the system using Rhyolite-MELTS (Table 2). Since the parental magma is not a primary liquid we assume a three-stage fractionation process with a starting water content of 0.75 wt% [85][86][87] . The first stage of fractionation that produced the model starting composition likely occurred in the uppermost mantle or lowermost crust and was probably dominated by olivine but may have included pyroxene and Cr-spinel 37 . The model begins at the second stage of fractionation (0.6 GPa) with an initial redox state at the FMQ buffer. The third stage of fractionation occurs at 0.1 GPa (~3.7 km) using the pre-spinel (magnetite) liquid from the 0.6 GPa model (at 1110 °C). The 'pre-spinel' liquid is used because the relative oxidation state does not affect the silicate mineral compositions of the crystallizing assemblage or residual liquid composition. The results show that liquids resembling the Masset silicic rocks can be generated at higher relative oxidation states (ΔFMQ > 0) as they yield magnesian, alkali-calcic to calc-alkalic compositions that are typical of volcanic-arc settings in spite of the fact that the volcanic suite formed at a rift-setting (Fig. 5) [35][36][37]80 . Furthermore, the shallow pressure (0.1 GPa) models show that ilmenite fractionates relatively early (1030 °C-1080 °C) in both the oxidizing and reducing models which encompasses a silica range of ~53.5 wt% to ~59.5 wt%.

Sample Composition
Amongst other geochemical criteria (trace element ratios, isotopes), the calc-alkaline signature of some silicic greenstone belt rocks is considered to be evidence in favour of a convergent margin setting for the bimodal volcanic suite 20,[88][89][90][91][92][93] . In order to test if the calc-alkaline silicic rocks from a greenstone belt can be derived by fractional crystallization, we apply the same modeling approach (H 2 O = 0.5 wt%, P initial = 0.5 GPa, ΔFMQ ± 2) as the Masset Formation to the tholeiitic sequence of the Archean Uchi-Confederation (Superior Craton) greenstone belt ( Table 2). The results, similar to the Masset Formation, show that the range of silicic rock compositions is reproduced under relative oxidizing (ΔFMQ > 0) conditions rather than reducing with ilmenite crystallizing only in the oxidizing (ΔFMQ +1, +2) models (Fig. 6). The implication is that the calc-alkaline signature of intermediate to silicic rocks from the Uchi-Confederation greenstone belt is likely related to the relative oxidation state of the magma system 81,82 . Although this does not preclude the possibility of a volcanic-arc setting for some bimodal sequences of greenstone belts, the models indicate that the presence of calc-alkaline rocks is not prima facie evidence of an arc setting. Moreover, it is possible that the depletion of Nb-Ta in some greenstone belt silicic rocks may be related to partitioning into ilmenite rather than as consequence of contamination by an 'enriched' source. Therefore, in order to fully utilize trace elements to constrain the tectonic setting of the intermediate to silicic volcanic rocks from a given greenstone belt, one must consider the influence that magmatic conditions (i.e., water content, relative oxidation state) impart on the melt during its evolution.

epithermal Gold Deposits of Haida Gwaii
Archean greenstone belts are major sources of gold 94,95 . For example ~80% of gold from Canada is mined from Archean greenstone belts of the Slave and Superior Provinces 96 . The types of gold deposits in greenstone belts are varied but typically fall under the terms 'orogenic' or 'lode' gold deposits 94,97,98 . Gold mineralization in greenstone belts is contextualized in terms of the relationship between the genesis of ore-bearing fluids and terrane accretion, thrusting, crustal shortening, metamorphism, and syn-orogenic extension. The style of gold mineralization on Haida Gwaii bears some resemblance to that of the Abitibi greenstone belt with respect to the influence of normal and listric faults as 'mineralization corridors' and to the emplacement of intermediate intrusions during extension 99 . However, Haida Gwaii has not experienced significant orogenic cycles or accretionary tectonics which could be a reason why there are no major gold deposits. www.nature.com/scientificreports www.nature.com/scientificreports/ There are a number of different types (skarn, epithermal Au, coal, perlite, alluvial sand) of mineral deposits identified on Haida Gwaii 100 . Iron, Cu and garnet skarns within the Kunga Formation limestone and the Karmutsen volcanics represent the largest (~34%) portion of deposits but there is an appreciable amount of epithermal gold deposits (38 known deposits). The gold deposits are often structurally associated with faults and/or Paleogene intrusive rocks (Kano intrusions). Higher grade mineralization is hosted by quartz veins and veinlets but there are prospects associated with hot springs as well.
The epithermal vein deposits are mostly located on Moresby Island (southern Haida Gwaii) and primarily hosted by Yakoun Group rocks and Karmutsen volcanic rocks with one occurrence within fault breccia of Masset and Karmutsen volcanic rocks. The formation of the epithermal deposits is associated with the emplacement of the Masset Formation intrusive rocks. The plutonic rocks were emplaced within older factures and faults that probably permitted remobilization and concentration of Au and Ag in silicic fluids that formed dense networks of veins during melt migration. Furthermore, the silicic magmas likely acted as the heat source for the hot spring systems.
The hot spring associated deposits are hosted by chalcedony and fine grained quartz veins and have Au-Ag ratio of 1:10. Many of the hot spring deposits exploited pre-existing faults or fractures and likely remobilized gold as they are often structurally above the older epithermal systems and hosted by Masset Formation and Skonun Formation. In some instances the sedimentary and volcanic rocks of the Middle Jurassic Yakoun Group also host deposits as well as the limestone and sedimentary rocks of the Kunga Group.

A Modern Analogue of an Archean Greenstone Belt
Many modern examples of greenstone belts are found along the convergent margins of the circum-Pacific region (Western Pacific and Cordilleran belts). This is, in part, due to their similarities to ophiolites but also due to the presence of boninites, calc-alkaline silicic rocks, and magnesian andesites in the arcs and accretionary complexes of Japan, Philippines and the islands of the Izu-Bonin-Mariana trench system 43,101-105 . Moreover, there are a significant www.nature.com/scientificreports www.nature.com/scientificreports/ number of large oceanic plateaux in the circum-Pacific region that are thought to be analogous to the lower volcanic series of greenstone belts due to their similarities in the stratigraphy, presence of subaqueous mafic-ultramafic volcanic rocks, possible association with mantle plumes, and absence of calc-alkaline rocks 8,9,[106][107][108] . The greater abundance of oceanic plateaux in the Pacific Ocean implies there is an increased likelihood that they could act as substrate for arc initiation or act as the nucleation site of the bimodal volcanic series [109][110][111][112] . However, given the uncertainty of Archean tectonics with respect to plate tectonics, and the geological complexity of greenstone belts in general, finding a modern analogue is problematic as it is likely there are a number of different 'types' of greenstone belts that were formed at markedly different settings 2,4,15,16,31,32,34 . At the moment, most modern examples of greenstone belts are associated with subduction zones and/or oceanic plateaux.
The rock units and stratigraphy of the Masset Formation along with the older volcanic basement sequence of the Karmutsen Formation bear a strong resemblance to the structure and lithologies of the volcanic portion of an Archean greenstone belt. The Karmutsen Formation is similar to the 'oceanic-plateau-like' basement of the ultramafic-mafic series of greenstone belts whereas the Masset Formation is similar to the bimodal volcanic series (Figs 3, 5 and 6). In the case of Haida Gwaii, both volcanic series were produced under tensional stress at a rift setting which may also be true for some greenstone belts. The magmatic rocks evolve from more primitive ultramafic compositions to evolved mafic and silicic compositions over a period of ~200 million years which is within the range (~50 Ma to 300 Ma) of many Archean greenstone belts 53,54 . The chemical evolution in the upper series is mostly controlled by crystal fractionation (Fig. 5). The calc-alkaline signature of the Masset Formation silicic rocks is due to crystallization under relatively oxidizing conditions (ΔFMQ + 0.7). Moreover, the calc-alkaline signature of the silicic rocks from the Uchi-Confederation greenstone belt may be due to the relative oxidation state of the magma during crystallization (Fig. 6). Thus, the presence of calc-alkaline rocks is not in and of itself definitive evidence for a subduction-related tectonic setting for Archean or younger greenstone belts. Finally, the presence of epithermal gold deposits within all volcanic units of Haida Gwaii is reminiscent of  Table 2. www.nature.com/scientificreports www.nature.com/scientificreports/ extension-related mineralization within the Abitibi greenstone belt 99 . Therefore, we propose that Haida Gwaii is an exemplar Phanerozoic analogue of subduction-unrelated Archean greenstone belts.

Methods
Mantle potential temperature calculations. The primary melt compositions and mantle potential temperature estimates were calculated using PRIMELT3 72 . The major elemental data of each sample was entered into PRIMELT3 and calculated using an Fe 2 O 3 /TiO 2 ratio of 0.5 and 1.0, pressure of 1 bar, H 2 O = 0 wt% and the lowest possible FeO content. The rock compositions and accumulated fractional melting (AFM) results are reported in Table 1.
Rhyolite-MELTS is optimization for silicic magma evolution. The software permits the user to vary magmatic parameters such as relative oxidation state (ƒO 2 ), pressure (GPa) and water (wt%) content so that multiple hypotheses can be tested. The starting compositions, initial magmatic conditions, and third stage compositions are reported in Table 1. We assumed initial water content of 0.75 wt% for the Masset models and 0.5 wt% for the Uchi-Confederation models and varied the relative oxidation state from ΔFMQ −2 to ΔFMQ +2. In order to simulate polybaric fractional crystallization the parental magmas are initially fractionated at high pressure (Masset = 0.6 GPa; Uchi-Confederation = 0.5 GPa) until 1110 °C. The 1110 °C composition was selected because it was the last liquid before Fe-Ti crystallization in both models. The volume of liquid remaining for the 1110 °C composition in the Masset models is 42.5% and 35.6% in the Uchi-Confederation models. The liquid composition at 1110 °C is then fractionated at a pressure 0.1 GPa (~3.7 km depth) until the model runs to completion (Stage-3). The output data files of each model are reported in the online supplementary material (Datasets S1 and S2).