Metamorphic P–T paths of Archean granulite facies metasedimentary lithologies from the eastern Beartooth Mountains of the northern Wyoming Province, Montana, USA: constraints from quartz-in-garnet (QuiG) Raman elastic barometry, geothermobarometry, and thermodynamic modeling

Metamorphic pressure and temperature (P–T) paths in late-Archean high-grade rocks of the eastern Beartooth Mountains of Montana (USA), a portion of the Wyoming Province, are established by a combination of imaging, analytical, and modeling approaches. Garnet inclusion mechanical and chemical thermobarometry, applied to several granulite-facies migmatites and an iron formation, is particularly useful in constraining the prograde P–T conditions. Quartz-in-garnet (QuiG) elastic Raman barometry was used on quartz inclusions in garnet for all samples studied. For a smaller subset of four representative samples, QuiG constraints were used in conjunction with Ti-in-quartz (TitaniQ) and Ti-in-biotite (TiB) thermometry to establish unique prograde inclusion entrapment P–T conditions. Ti measurements of garnet hosts and cathodoluminescence (CL) imagery of inclusion and matrix quartz grains were employed to check for Ti loss/diffusion. Lastly, inclusion studies were supplemented with thermodynamic modeling and matrix chemical thermobarometry to examine garnet nucleation temperatures and peak metamorphic conditions.

Disagreement between the volume strain and elastic tensor methods used to calculate quartz inclusion pressures implies that quartz inclusions studied are under strong differential strain. Prograde entrapment results from the two inclusion thermobarometry pairs used are distinct: 0.55–0.70 GPa and 475–580 ^∘C (QuiG–TitaniQ) versus 0.85–1.10 GPa and 665–780 ^∘C (QuiG–TiB). Garnet modal isopleth modeling indicates that the majority of garnet growth occurred at ∼ 450–600 ^∘C, implying that P–T conditions of garnet growth are interpreted to be most reliably represented by QuiG–TitaniQ inclusion thermobarometry. Normal distributions of calculated QuiG inclusion pressures and the concentration of mineral inclusions in garnet cores suggest that the majority of garnet inclusions were entrapped during a single stage of porphyroblast growth. A general lack of evidence from CL imagery for post-entrapment mechanical or chemical modifications to quartz inclusions suggests that quartz inclusions used to calculate entrapment P–T largely preserve their initial entrapment conditions. Biotite inclusions preserve higher temperatures than quartz inclusions in the same garnets, likely due to Fe–Mg exchange with garnet hosts that allowed Ti content of biotite to change after entrapment. Pseudosection modeling and matrix chemical thermobarometry of multiple, independent lithologies examined during inclusion studies suggest a range of peak granulite facies conditions of ∼ 0.50–0.70 GPa and 730–800 ^∘C. Peak metamorphic P–T modeling work from this study, together with evidence of regional amphibolite facies overprinting in immediately adjacent samples, indicates an overall clockwise metamorphic P–T path with nearly isobaric prograde heating to peak temperatures. Interpreted P–T path reconstructions are consistent with metamorphism developed in a more modern-style continental arc subduction zone and are observed in portions of the northern Wyoming Province as exemplified by metasupracrustal lithologies of the eastern Beartooth Mountains.

1 Introduction

Pressure and temperature (P–T) conditions during prograde metamorphism are commonly not well-constrained in metamorphic rocks due to a lack of thermobarically diagnostic inclusion mineral assemblages in porphyroblast hosts and/or re-equilibration of the matrix assemblage during peak metamorphism or later overprints. Prograde conditions are commonly obscure in granulite facies metamorphic rocks where simply discerning individual thermal events can be difficult due to resetting of mineral thermometers upon cooling, shrinking of sample equilibrium volumes at fluid-absent conditions, and/or partial to complete melt loss (e.g., Frost and Chacko, 1989; Pattison et al., 2003; Guevara and Caddick, 2016). Despite the elusive nature of prograde constraints, determining geologic conditions during prograde metamorphism provides for more accurate overall P–T path reconstructions that may be used to determine tectonometamorphic setting (e.g., Thompson and England, 1984).

A common approach to investigating prograde metamorphism is the use of mineral inclusions as recorders of the conditions at which they were entrapped. This approach relies on the assumption that mineral hosts act as chemical and/or mechanical insulators for inclusion phases upon entrapment, preventing the inclusions from re-equilibrating, either elastically or chemically, in response to evolving metamorphic conditions (e.g., Kohn, 2014). Inclusion studies have proven to be considerably useful by allowing workers to probe otherwise inaccessible portions of metamorphic P–T paths (e.g., Chopin, 1984). Several host–inclusion pairs have been utilized to constrain conditions of inclusion entrapment (e.g., olivine-in-diamond, Izraeli et al., 1999; coesite-in-diamond, Sobolev et al., 2000; garnet-in-diamond, Milani et al., 2015; feldspar-in-garnet, clinopyroxene, and olivine, Befus et al., 2018; apatite-in-garnet, Barkoff et al., 2019; zircon-in-garnet, Campomenosi et al., 2020; quartz-in-epidote, Cisneros et al., 2020; apatite-in-olivine, rutile-in-garnet, feldspar-in-epidote, Cisneros and Befus, 2020; quartz-in-zircon, Gonzalez et al., 2021; apatite-in-zircon, Guo et al., 2021; quartz-in-kyanite, Tomioka et al., 2022). Because of the prevalence of quartz inclusions in garnet porphyroblasts whose formation conditions can span a wide swath of P–T space, probably the most popular host–inclusion pair is quartz-in-garnet (e.g., Enami et al., 2007; Spear et al., 2014; Thomas and Spear, 2018; Alvaro et al., 2020; Spear and Wolfe, 2020; Gilio et al., 2021b; Mulligan et al., 2022). Improvements to elastic models and to the respective equations of state (EoS) of quartz and garnet have established quartz-in-garnet (QuiG) as a functional elastic barometer useful for exploring the conditions of quartz inclusion entrapment during garnet growth, commonly during prograde metamorphism (Angel et al., 2014, 2017; Campomenosi et al., 2018; Mazzucchelli et al., 2018; Murri et al., 2018; Angel et al., 2019, 2022). Ti-in-quartz (TitaniQ) thermometry (Thomas et al., 2010) has been paired with QuiG barometry to discriminate prograde temperatures of quartz inclusion entrapment (e.g., Gonzalez et al., 2019). Although there is uncertainty about co-entrapment of quartz and non-quartz mineral inclusions, non-quartz mineral inclusion thermometry (e.g., Ti-in-biotite (TiB) thermometry, Henry et al., 2005) may still provide complementary information regarding evolving P–T conditions during progressive metamorphism.

P–T paths in Archean high-grade metamorphic rocks are particularly enigmatic but important. Archean terranes represent an instructive avenue for the application of QuiG barometry as published work on inclusion elastic thermobarometry in Archean rocks is scarce (e.g., Parkinson and Katayama, 1999; Guo et al., 2021). The fragmentary and poorly preserved nature of the Archean rock record, the only indicator of the nature and timing of tectonic processes in the early Earth, has resulted in considerable debate regarding a number of aspects about the Archean earth (e.g., Korenaga, 2013; Brown and Johnson, 2018; Stern, 2018). Applying inclusion studies to Archean metamorphic rocks is an additional tool for resolving complex P–T histories characteristic of Archean tectonometamorphic processes.

The purpose of this study is to determine the prograde metamorphic conditions of multiple granulite-facies, garnet-bearing metasedimentary rocks from the eastern Beartooth Mountains, Montana, in the northern Wyoming Province. Specifically, this study seeks to establish prograde conditions of garnet growth using QuiG barometry combined with prograde temperature constraints from TitaniQ and TiB thermometry. Thermodynamic modeling of reconstructed sample bulk compositions (i.e., garnet modal isopleth diagrams) will be used to constrain the conditions of the garnet-in reaction and subsequent garnet growth as a supplement to inclusion-thermobarometry-based P–T estimates of garnet growth. Lastly, this study will utilize pseudosection modeling and chemical thermobarometry of the matrix assemblage to calculate peak conditions of metamorphism. The combination of P–T constraints from inclusion and matrix studies will allow for determination of more robust P–T paths. Rocks from the eastern Beartooth Mountains preserve a complex, polymetamorphic history whose geologic origins are tied to some of the earliest indications of modern-style plate tectonic processes such as compressional tectonics and arc magmatism (e.g., Mueller et al., 2014). Therefore, the implications of this work are geologically significant as they serve to further illuminate the metamorphic history of rocks that preserve evidence of nascent plate tectonic processes during the middle to late Archean. Additionally, the results of this work are significant as they demonstrate the applicability of inclusion elastic thermobarometry to rocks of great antiquity.

2 Geologic background

The Beartooth Mountains of Montana and Wyoming lie within the Wyoming Province (WP), an Archean craton that spans multiple states of the western USA (Fig. 1). The WP has three primary subdivisions going from north to south: the Montana metasedimentary terrane (MMT), the Beartooth–Bighorn magmatic zone (BBMZ), and the southern accreted terranes (SAT), as shown in Fig. 1 (Mogk et al., 1992; Mueller and Frost, 2006). While the Beartooth range largely falls within the BBMZ, the northwestern corner lies within the MMT, a zone of tectonic mixing of allochthonous units juxtaposed along ductile shear zones (Mogk et al., 1988; Wooden et al., 1988; Mueller et al., 1996). The BBMZ lithologies are predominantly ∼ 2.80 Ga (meta)plutonic rocks of the tonalite–trondhjemite–granodiorite (TTG) suite with older (up to 3.50 Ga) xenolithic metaigneous and metasupracrustal enclaves and pendants that are tectonically intercalated (e.g., Henry et al., 1982; Mueller et al., 1988; Mogk et al., 2020, 2022).

Figure 1 Map of the western United States modified after Mogk et al. (2022) showing the three principal divisions of the Wyoming Province noted by thick black lines (i.e., MMT, BBMZ, and SAT) relative to surrounding geologic blocks. Petrologic divisions of the Beartooth Mountains are labeled (i.e., NSB – North Snowy Block, SSB – South Snowy Block, SC – Stillwater Block, and BT – Beartooth Plateau block). State abbreviations are indicated where appropriate (e.g., MT – Montana, WY – Wyoming).

The WP is geologically significant because it records evolutionary trends of ancient continental crust from >4.00 Ga to ∼ 2.80 Ga (e.g., Mueller et al., 1992; Mueller and Frost, 2006; Mogk et al., 2022; Frost et al., 2023). Evidence of the oldest crust is inferred from Lu–Hf isotopic data of detrital or xenocrystic zircons found in orthoquartzites and gneisses. These data indicate formation of mafic proto-crust developing at ∼ 4.00–3.50 Ga via plume-dominated tectonic processes in a “stagnant-lid” regime from a primitive or slightly depleted mantle (Mueller and Wooden, 2012). The oldest rocks are TTG gneisses that formed starting at 3.60 Ga with major additional crust formation at 3.30–3.10 Ga. These gneisses have geochemical affinities indicative of a continental arc setting marking the “plume-to-plate” transition (Mueller et al., 2014; Mogk et al., 2022). Sediments such as arenitic quartzite, aluminous shale, psammite, and iron-rich sediments were deposited at 3.00–2.80 Ga (i.e., time between youngest quartzite detrital zircon age and age of enveloping TTG rocks), presumed to be in a passive margin setting. These sedimentary rocks were then metamorphosed to upper amphibolite to granulite facies conditions and tectonically mixed with metabasic and meta-ultramafic rocks as well as older TTG gneisses in pervasive and anastomosing high-grade ductile shear zones (Mogk et al., 2022). A second major crust-forming event occurred at 2.83–2.79 Ga with the production of voluminous calc–alkaline plutonic rocks, diorite to granite in composition, attributed to contemporary-style continental arc tectonism involving subduction and concurrent melting of the down-going slab, mantle wedge, mafic underplate, and deeply buried metasedimentary rocks (e.g., Mueller and Frost, 2006; Mueller et al., 2010, 2014).

The Beartooth Mountains are further divided into four distinct blocks based on petrologic character (Fig. 2) (Wilson, 1936). The North Snowy Block, part of the BBMZ, is a zone of tectonic mixing of kilometer-scale allochthonous units of metamorphosed igneous and/or sedimentary lithologies with the oldest being 3.50 Ga and tectonic juxtaposition occurring prior to 2.55 Ga (Mogk et al., 1988; Mueller et al., 1996). The South Snowy Block is principally comprised of metapelitic and metapsammitic rocks of the Jardine metasedimentary sequence (JMS) intruded by ∼ 2.80 Ga granitic plutons (Casella et al., 1982; Thurston, 1986; Mogk et al., 2012). The extensively studied Stillwater complex, a 2.71 Ga layered mafic–ultramafic intrusion, and its contact metamorphic aureole typify the Stillwater Block (e.g., Page, 1977). Of primary importance to this study is the Beartooth Plateau block. The Beartooth Plateau is characteristic of the BBMZ, i.e., volumetrically dominated by 2.80 Ga TTG rocks, but with notable enclaves of older metamorphosed igneous and sedimentary lithologies (e.g., Henry et al., 1982; Mueller et al., 2010).

Figure 2 Map of the Beartooth Mountains modified after Mueller et al. (1992). Sample localities for rocks studied here are marked in black along the eastern fringe of the Beartooth Plateau block (i.e., Hellroaring Plateau – HP and Quad Creek Plateau – QC). The inset shows the position of the range with respect to state boundaries whose abbreviations are noted.

Peak metamorphic conditions of the metasedimentary xenolithic rocks of the eastern portion of the Beartooth Plateau block have been well-constrained, but the prograde metamorphic path and retrograde overprints are less certain. Based on multiple thermobarometry models, petrogenetic grids, and pseudosection calculations, peak metamorphism (M1) reached upper amphibolite to granulite facies at ∼ 0.60–0.80 GPa and 775–800 ^∘C, representing pre-intrusion conditions at ∼ 2.80 Ga (e.g., Maas, 2004; Will, 2013). Subsequently, these lithologies were tectonically assembled and intensely folded at high grades at ∼ 740 ^∘C (Henry et al., 1982; Henry and Daigle, 2018; Mogk et al., 2022). The TTG plutonic suite rocks that dominate the BBMZ, typified by the Long Lake Magmatic Complex (LLMC), are derived from multiple, compositionally diverse mantle and crustal sources emplaced over a ∼ 40 Ma period (2.79–2.83 Ga) in a modern-style oceanic–continental convergent margin environment (Mueller et al., 2010). Intrusion of LLMC-like magmas over a ∼ 40 Ma period likely resulted in overprinting of the granulite facies peak assemblage by a regional amphibolite facies event (M2) ∼ 2.74–2.79 Ga at ∼ 0.50–0.60 GPa and 650–700 ^∘C with localized dehydrating fluids (i.e., fluids with a lowered activity of water) producing later granulite-facies mineral assemblages (Mueller et al., 1988, 1992, 1998; Henry et al., 2015; Henry and Daigle, 2018). A middle to upper greenschist facies hydrating overprint (M3) at ∼ 0.30–0.40 GPa and 450–500 ^∘C affects most rocks to variable degrees (e.g., serpentine, talc, and chlorite in ultramafic xenolithic lithologies and actinolite after hornblende in mafic xenolithic lithologies) (Henry et al., 2015; Henry and Daigle, 2018).

Select garnet-bearing xenolithic lithologies within the TTG (meta)igneous rocks from the eastern margin of the Beartooth Mountains are examined in this study. The polymetamorphic trends characteristic of these garnet-bearing rocks of interest commonly make geologic P–T reconstructions difficult and equivocal. For example, earlier potential constraints on the prograde metamorphic P–T paths in the eastern Beartooth xenolithic rocks have varied. Maas (2004) inferred a near-isobaric path to ∼ 750 ^∘C followed by a counterclockwise P–T path at higher T on the basis of mineral–chemical thermobarometry and sillimanite growth textures. In contrast, Guevara et al. (2017) argue for a clockwise P–T path as a result of near-isobaric heating to peak temperatures primarily using bulk-rock thermobarometry and garnet major element zoning models. The aim of this study is to distinguish between these contrasting options by analyzing the strain state and chemistry of inclusions in garnet as a proxy for prograde metamorphic conditions.

3 Methods

Samples analyzed in this investigation include garnet-bearing metasupracrustal rocks from the eastern Beartooth Plateau block, specifically migmatite and iron formation xenolithic lithologies. A total of 15 samples (i.e., 14 migmatites and 1 iron formation) were initially examined via Raman spectrometry using 100 µm thick polished thin sections to establish whether residual strains are preserved in quartz inclusions in garnet. Mineral chemical investigations were done on four select samples (i.e., three migmatites and the iron formation) with measurable quartz inclusion strains to provide additional constraints on P–T conditions.

3.1 Raman spectrometry

Inclusion and matrix quartz grains were measured on a Renishaw inVia Reflex Raman microscope at the Shared Instrumentation Facility (SIF) at Louisiana State University (LSU) for entrapment pressure determination using quartz-in-garnet (QuiG) barometry. Inclusion measurements primarily targeted the centers of rounded, fully enclosed, spherical to ovoid, non-faceted quartz inclusions located away from interfaces (e.g., host fractures, grain boundaries) in accordance with model assumptions (e.g., Zhang, 1998; Campomenosi et al., 2018; Mazzucchelli et al., 2018). Raman spectra were collected at ambient pressures and temperatures with a 532 nm incident laser in the center of each quartz inclusion using a 100× microscope objective in high-confocality mode. Experiments utilized an 1800 groove mm⁻¹ grating; single shot laser exposure times ranged from 20–30 s. The instrument was calibrated using a Ne wavelength calibration source and an internal Si metal standard. Data quality was monitored using the 520 cm⁻¹ peak of a Si metal wafer. Peak positions of the 128 cm⁻¹ (i.e., 127.5 cm⁻¹), 206 cm⁻¹ (i.e., 205.9 cm⁻¹), and 464 cm⁻¹ (i.e., 464.8 cm⁻¹) vibrational modes of quartz were assigned upon spectral acquisition using the Renishaw software; peak assignment error is approximately 0.2–0.3 cm⁻¹. Inclusion quartz Raman peak shifts were calculated by comparing the positions of the main spectral bands for inclusion quartz grains relative to their average positions in matrix quartz grains. The positions of the three main matrix quartz bands among measured representative samples, as well as the position of the 520 cm⁻¹ band of the Si metal wafer standard, can be found in the Supplement (Table S3).

Relative peak shifts were used to calculate inclusion strains with the Grüneisen tensor parameters for quartz and the stRAinMAN software package (Murri et al., 2018; Angel et al., 2019). Inclusion pressure (P_inc) and associated entrapment isomekes, lines in P–T space of possible inclusion entrapment, were calculated from inclusion strains using the online EntraPT platform (Mazzucchelli et al., 2021). P_inc values were calculated using a volume strain approach ($P_{\mathrm{inc}}^{\mathrm{vs}}$) that assumes quartz inclusions are under isotropic strains and an elastic tensor ($P_{\mathrm{inc}}^{\mathrm{et}}$) approach that accounts for any strain anisotropy in measured inclusions. While quartz, a non-isotropic hexagonal mineral, is elastically anisotropic (e.g., Murri et al., 2018), the former inclusion pressure calculation method was employed to monitor the effect of inclusion strain anisotropy on calculated inclusion pressures. EntraPT calculations utilized the elastic tensor for quartz at room pressure and temperature (Wang et al., 2015), an almandine third-order Birch–Murnaghan EoS (Angel et al., 2022) for host volume calculations, and a quartz Landau EoS with a curved α-β quartz boundary (Angel et al., 2017) for inclusion volume calculations.

3.2 Mineral chemistry

Mineral chemistry of inclusion and matrix grains was analyzed via wavelength-dispersive spectrometry (WDS) on a JEOL 8230 electron microprobe (EMP) in the Chevron Geomaterials Characterization Lab at LSU in the four representative samples. EMP analytical objectives included trace element chemistry (i.e., Ti content) of quartz and garnet and major element chemistry of garnet, biotite, plagioclase, and pyroxene. EMP data quality was ensured using well-characterized natural and synthetic standards. More information regarding EMP operating parameters, standards, and data quality can be found in the Supplement (File S1). Ti was measured in quartz inclusions exposed at the garnet–thin-section interface as well as in quartz grains in the matrix. Cathodoluminescence images observed with the electron microprobe (EM-CL) were collected in quartz grains prior to Ti measurements. EM-CL imagery was performed using a blue light filter as the blue region of the CL spectrum generally correlates to quartz Ti content (Spear and Wark, 2009). Mineral formulae were calculated from EMP chemistry data via normalization procedures based on the appropriate number of oxygens per formula unit (12 oxygens for garnet, 22 oxygens for biotite where OH + F + Cl = 4, 8 oxygens for plagioclase, and 6 oxygens for pyroxene). Fe³⁺ was calculated based on charge balance and stoichiometry.

Mineral chemistry information was used to calculate inclusion entrapment temperatures using two Ti-based thermometers: Ti-in-quartz (TitaniQ) thermometry and Ti-in-biotite (TiB) thermometry. TitaniQ and TiB inclusion temperatures were specifically utilized to isolate unique entrapment conditions along generated QuiG entrapment isomekes. While the latter thermometer assumes an equilibration pressure of ∼ 0.40–0.60 GPa and the presence of graphite, a typically absent phase in high-grade rocks, TiB temperature estimates can be used as limiting values (Henry et al., 2005). Both TitaniQ and TiB thermometers require a Ti-saturating phase, commonly rutile due to its capability as a buffer to fix the activity of titania (${\mathit{\alpha }}_{{\mathrm{TiO}}_{\mathrm{2}}}$) at unity. Metasupracrustal rocks in the study area generally contain ilmenite as the primary Ti-bearing phase; instances of rutile are typically the result of secondary alteration of ilmenite. In the absence of a Ti-buffering phase ensuring that ${\mathit{\alpha }}_{{\mathrm{TiO}}_{\mathrm{2}}}$ is equal to 1, ${\mathit{\alpha }}_{{\mathrm{TiO}}_{\mathrm{2}}}$ was assumed to be 0.9, a reasonable assumption for ilmenite-bearing metapelites (Ghent and Stout, 1984). Furthermore, ${\mathit{\alpha }}_{{\mathrm{TiO}}_{\mathrm{2}}}$ likely remained at this value throughout metamorphism given that ilmenite is found as both a matrix and inclusion phase in garnet for nearly all samples examined.

In addition to thermometers in inclusion studies, mineral chemistry was also used for matrix chemical thermobarometry to supplement calculated peak metamorphic P–T from pseudosection modeling. Specifically, the mineral thermobarometers utilized include the garnet–biotite Fe–Mg exchange thermometer (Holdaway, 2000), the garnet–biotite–plagioclase–quartz (GBPQ) net transfer barometer (Wu et al., 2004), the intracrystalline garnet Fe²⁺–Ca²⁺ exchange barometer (Wu, 2019), the net-transfer garnet–biotite–aluminosilicate–quartz (GBAQ) barometer (Wu, 2017), and the garnet–aluminosilicate–plagioclase (GASP) net-transfer barometer (Hodges and Crowley, 1985). Garnet-based barometers were applied with temperature estimates from TiB thermometry. Although the GBAQ and GASP barometers require modal aluminosilicate, a phase that is not present in all migmatites and a phase that is not present at all in the iron formation, calculated pressures may still provide geologically useful pressure constraints. In addition to these matrix thermobarometers, TitaniQ and TiB thermometers were also respectively applied to quartz and biotite grains in the matrix.

3.3 Thermodynamic modeling

Metamorphic phase equilibria were modeled in the samples of interest using the Theriak-Domino software package (De Capitani and Petrakakis, 2010). Thermodynamic modeling primarily utilized bulk sample chemistry reconstructed from thin section mineral modes and EMP mineral chemistry data using the Rock Maker program (Büttner, 2012). Mineral modes in samples of interest were constrained from both visual estimates and percentages calculated using a Keyence digital microscope. A representative migmatite (sample HP 82-88b) and the single iron formation sample (sample HP 81-82) examined were modeled in the MnNCKFMASHTO and MnCKFMASHTO systems, respectively, using the Holland and Powell (2011) dataset 6.2 formatted for Theriak-Domino by Doug Tinkham. Activity–composition models used included garnet, biotite, orthopyroxene, cordierite, staurolite, chlorite, chloritoid, white mica, and melt (White et al., 2014a, b); ternary feldspar (Holland et al., 2022); clinoamphibole (Green et al., 2016); and ilmenite and magnetite (White et al., 2000). Water content for thermodynamic models was principally inferred from the modes of hydrous peak assemblage minerals present. Inferred moles of H₂O were found to be broadly consistent with unpublished loss-on-ignition data from similar rocks in the study area. Inferring water content in this way is advantageous as it avoids adding low-temperature hydrous alteration products to the modeled bulk composition (e.g., pinite after cordierite) that would otherwise be included in bulk-rock measurements from X-ray fluorescence. Three additional moles of water were added to the iron formation bulk composition for pseudosection modeling; this is because Domino was unable to detect any reactions with the essentially anhydrous composition calculated by mineral modes. The incorporation of additional water into the modeled bulk composition is not unwarranted given evidence of high-grade fluids attendant during peak metamorphism of iron formation samples at an albeit reduced activity of water (Henry and Daigle, 2018). Migmatite models assumed a ${\mathrm{Fe}}^{\mathrm{3}+}/{\mathrm{Fe}}_{\mathrm{total}}$ ratio of 0.10 due to the presence of ilmenite without exsolved hematite as the primary oxide phase. More oxidizing conditions were inferred for iron formation samples due to the considerable modal abundance of magnetite; thus a higher ${\mathrm{Fe}}^{\mathrm{3}+}/{\mathrm{Fe}}_{\mathrm{total}}$ ratio of 0.25 was used in model calculations. Iron formation efforts also included modeling of the bulk composition with modal magnetite factored out. For these calculations, essentially all iron was assumed to be ferrous.

Thermodynamic models generated consisted of (1) garnet volume isopleth diagrams to identify the temperature range of the garnet-in reaction and subsequent garnet growth and (2) isochemical pseudosection diagrams to determine peak metamorphic conditions proceeding prograde entrapment conditions. Increasing amounts of H₂O (i.e., 12.5 mol followed by 25 and 50 mol) were increasingly added to the calculated bulk composition for construction of garnet volume isopleth diagrams to simulate the elevated water content likely present at lower metamorphic grades during garnet growth.

The representative samples chosen for modeling in Domino were migmatite sample HP 82-88b and iron formation sample HP 81-82. HP 82-88b contains the assemblage biotite–cordierite–quartz–plagioclase–garnet–ilmenite. On the other hand, HP 81-82 contains the assemblage orthopyroxene–garnet–quartz–magnetite–biotite–ilmenite. As a result of there being no major phases in sample HP 81-82 that contain significant Na, thermodynamic modeling was carried out in the MnCKFMASHTO compositional system. As is typically the case in Beartooth samples, the mineral constituents that comprise both samples do not exhibit major element zoning. Mineral compositions inputted into the RockMaker program are representative averages of matrix grains that were proximally located to one another and not in immediate contact. The modeled migmatite mineral modes and associated mineral composition represent an average of the leucosome and melanosome components of the rock as preserved at the thin section scale. Reconstructed bulk compositions calculated from mineral chemistry data for these two samples are shown in Table 1.

Analytical points from the outermost rims of garnet porphyroblasts were ignored for grain averaging due to retrograde Fe–Mg exchange with surrounding ferromagnesian phases like biotite. Cordierite grains in sample HP 82-88b are extensively pinitized; thus mineral chemistry information was not directly obtained. Instead, a cordierite composition with an X_Mg of 0.60 was assumed, which is consistent with the composition of non-pinitized cordierite grains in similar migmatitic rocks from the study area. Quartz (SiO₂), ilmenite (FeTiO₃), and magnetite (Fe₃O₄) were all treated as pure phases in bulk sample reconstructions.

Table 1 Reconstructed bulk-rock chemistry in weight percent oxide of the migmatite (HP 82-88b) and the iron formation (HP 81-82) modeled in Domino. Modal estimates for HP 82-88b are biotite (30.0 %), cordierite (27.0 %), quartz (19.5 %), plagioclase (13.5 %), garnet (8.0 %), and ilmenite (2.0 %). Original modal estimates for HP 81-82 are orthopyroxene (42.5 %), garnet (34.5 %), quartz (11.5 %), magnetite (9.0 %), biotite (2.0 %), and ilmenite (0.5 %).

4 Results

4.1 Mineralogy and host–inclusion relations in migmatite samples

Migmatitic rocks studied ranged from metatexites (i.e., display leucosomes and melanosomes as discrete layers) to diatexites (i.e., display intermixed leucosome and melanosome material). In both types of migmatites, melt is interpreted to have been derived locally via dehydration melting of hydrous phases as the pelitic to psammitic protoliths were heated to granulite facies conditions during M1 metamorphism (Maas, 2004). Minerals such as biotite and/or sillimanite typically define the main foliation as well as minor crenulations in samples. Melanosome minerals are typically comprised of garnet, biotite, and ± cordierite and/or sillimanite, whereas leucosomes are chiefly comprised of quartz, plagioclase, and ± alkali feldspar. Common trace minerals include ilmenite, apatite, zircon, monazite, magnetite, and pyrite. The most prevalent alteration products in migmatitic rocks are partial to complete pinitization of cordierite and partial sericitization of plagioclase. The mineralogy of the representative migmatite samples studied is given in Table 2.

Garnet hosts/porphyroblasts in migmatitic rocks vary widely in diameter (i.e., <1.0 to >10.0 mm), modal abundance, and overall quality (i.e., degree of fracturing and embayment). As is the case with the other major mineral constituents in migmatites, garnets studied do not typically display significant major element zoning. Porphyroblasts predominantly display flat interior compositional profiles with slight zoning along rims due to late Fe–Mg exchange (Fig. 3).

Table 2 Mineralogy of the four representative samples examined for mineral chemistry. × denotes minerals that are present above approximately 2 modal percent, whereas tr denotes trace minerals whose modes are below 2 modal percent. Alteration products are italicized. Mineral abbreviations used in the table and throughout this study are after Whitney and Evans (2010).

Figure 3 Compositional traverse of a garnet porphyroblast in migmatite sample HP 82-88b shown with respect to mole percent (mol%) of the four relevant garnet endmembers. Note the generally flat interior compositional profile and the Fe–Mg exchange at the rims marked by enrichment in almandine component and decrease in pyrope component.

In most samples, garnet hosts are commonly surrounded by varying amounts of biotite (± fibrolitic sillimanite). Additionally, garnet hosts sit atop rock fabrics in migmatite samples such as foliations and do not appear rotated. All porphyroblasts are typically inclusion-rich, and in many cases inclusions are concentrated within garnet cores and are generally not strongly aligned (e.g., Fig. 4). Quartz is the most ubiquitous inclusion phase along with lower amounts of biotite, ilmenite, magnetite, sillimanite, and various high-relief phases (e.g., apatite, monazite, and zircon). Regardless of size, quartz inclusions are most commonly spherical to ovoid in shape.

Figure 4 (a, b) Plane-polarized light photomicrographs of garnet porphyroblasts from migmatite sample HP 82-64 that show inclusion-rich cores and inclusion-poor rims typical of the samples examined. These samples exemplify the typical quartz-dominant nature of garnet inclusion mineralogy. Additionally, these photos show the generally spherical to ovoid morphology of quartz inclusions studied. The red scale bars represent 1 mm.

4.2 Mineralogy and host–inclusion relations in iron formation samples

Iron formation samples in the study area have been interpreted to be ferruginous sediments that experienced granulite facies metamorphism (Henry and Daigle, 2018). Samples may be either massive or banded; banding, if present, is expressed by alternating layers of quartz-rich and magnetite-rich bands on the centimeter to millimeter scale. The sole iron formation sample studied in this investigation (sample HP 81-82) is a non-banded sample that is fine-grained and granoblastic. Sample mineralogy is predominantly comprised of an anhydrous mineral paragenesis of orthopyroxene, garnet, quartz, and magnetite (Table 2). Biotite, ilmenite, apatite, zircon, pyrite, and pyrrhotite are common trace constituents. Clusters of actinolite–cummingtonite–grunerite mixed with iron carbonates ankerite and siderite along sample fractures are the only major alteration product observed.

In the iron formation sample, garnet hosts/porphyroblasts are all essentially uniformly round and ∼ 0.5 to 1.0 mm in diameter with minor fracturing. Additionally, garnet porphyroblasts do not show major element chemical zoning (Fig. 5). As observed in many migmatite samples, garnet hosts commonly display inclusion-rich cores and inclusion-deficient rims (Fig. 6). Also similar to migmatite garnets, garnets in iron formation samples also sit atop any foliations and do not appear rotated. Garnet inclusions are not aligned along a pre- or syn-growth foliation and predominantly consist of spherical to ovoid quartz with variable magnetite and lesser ilmenite and biotite.

Figure 5 Compositional traverse of a garnet porphyroblast in iron formation sample HP 81-82 shown with respect to mole percent (mol%) of the four relevant garnet endmembers. The compositional profile mirrors both the generally flat profile from sample HP 82-88b as well as the slight increase of almandine component and slight decrease of pyrope component along the rims.

Figure 6 (a, b) Plane-polarized light photomicrographs of garnet porphyroblasts from iron formation sample HP 81-82. The fine-grained nature of iron formation samples is apparent, given the smaller size of garnet porphyroblasts relative to migmatite garnets. As noted in migmatite samples, inclusions are generally dominated by quartz and concentrated within garnet interiors. The red scale bars represent 1 mm.

4.3 QuiG barometry

From the full 15 samples measured in this investigation, there were approximately 230 quartz inclusions whose main Raman spectral peaks were shifted to higher wavelengths relative to unstrained quartz. Among the four representative samples, the proportion of measured quartz inclusions that preserved positive inclusion pressures varied from 20 % to 38 % (Table 3). Raman peak shift magnitude varied irrespective of inclusion location in garnet hosts within the peak assignment error.

Table 3 Results of Raman investigations for the four representative samples, specifically the average peak positions of the three main quartz spectral peaks as preserved by matrix grains, the number of inclusions measured (n), the number of measured inclusions that preserve positive inclusion pressures (+P), and the percent of inclusions measured that displayed positive inclusion pressures.

Volume strain estimates of inclusion pressure ($P_{\mathrm{inc}}^{\mathrm{vs}}$) were found to be higher than elastic tensor estimates of inclusion pressure ($P_{\mathrm{inc}}^{\mathrm{et}}$) for nearly all inclusions measured. However, the difference between $P_{\mathrm{inc}}^{\mathrm{vs}}$ and $P_{\mathrm{inc}}^{\mathrm{et}}$ was consistently within uncertainty for all inclusions but a few outliers (Fig. 7). Both $P_{\mathrm{inc}}^{\mathrm{vs}}$ and $P_{\mathrm{inc}}^{\mathrm{et}}$ values follow approximately normal distributions despite minor peaks at higher P_inc values (Fig. 8).

Figure 7 Comparison of $P_{\mathrm{inc}}^{\mathrm{vs}}$ and $P_{\mathrm{inc}}^{\mathrm{et}}$ constraints on inclusion pressures from migmatite samples (a) and the iron formation sample (b). The number of inclusions for each sample is denoted by n values. In the majority of inclusions shown, $P_{\mathrm{inc}}^{\mathrm{vs}}$ estimates of inclusion pressure are higher than $P_{\mathrm{inc}}^{\mathrm{et}}$ estimates of inclusion pressure (i.e., fall above the 1 : 1 line that represents agreement between $P_{\mathrm{inc}}^{\mathrm{vs}}$ and $P_{\mathrm{inc}}^{\mathrm{et}}$). However, this difference is typically within P_inc uncertainty denoted by the red triangle and accompanying x–y error bars along the 1 : 1 line.

The lack of an observed statistical difference between $P_{\mathrm{inc}}^{\mathrm{vs}}$ and $P_{\mathrm{inc}}^{\mathrm{et}}$ has been previously observed in quartz inclusions that, like the quartz inclusions from this study, preserve inclusion pressures below 1 GPa (Gilio et al., 2021a). Therefore, use of either inclusion pressure constraint would be permissible in this case if preserved inclusion strains were approximately isotropic. However, it is noted that inclusions studied preserve anisotropic strains given by the linear relationship between differential strain (i.e., difference between strain calculated along the a and c crystallographic axes: ε₁–ε₃) and the difference between $P_{\mathrm{inc}}^{\mathrm{vs}}$ and $P_{\mathrm{inc}}^{\mathrm{et}}$ (Fig. 9).

Figure 8 Distribution of elastic tensor inclusion pressures ($P_{\mathrm{inc}}^{\mathrm{et}}$, solid line) and volume strain inclusion pressures ($P_{\mathrm{inc}}^{\mathrm{vs}}$, dashed line) in samples HP 82-88b (a), QC 81-20 (b), HP 82-64 (c), and HP 81-82 (d). Calculated inclusion pressures in all four samples follow approximately normal distributions, ignoring minor apparent peaks at higher inclusion pressures.

Figure 9 Differential strain (ε₁–ε₃) versus the difference between inclusion pressure (P_inc) from volume strain ($P_{\mathrm{inc}}^{\mathrm{vs}}$) and inclusion pressure from the elastic tensor ($P_{\mathrm{inc}}^{\mathrm{et}}$) for migmatite inclusions (a) and iron formation inclusions (b). The number of inclusions for each sample is denoted by n values. In both lithologies, there is a linear relationship between differential strain and the difference between $P_{\mathrm{inc}}^{\mathrm{vs}}$ and $P_{\mathrm{inc}}^{\mathrm{et}}$.

Table 4 Mean inclusion pressure values from the volume strain method (with and without differential strain criteria applied) and the elastic tensor method applied to representative samples. Lithology averages are bolded. A slight difference is noted between volume strain and elastic tensor values without application of the anisotropy criteria (i.e., column two and column four). This difference between volume strain and elastic tensor inclusion pressure increases when anisotropy criteria are applied (i.e., column three and column four). Note that $P_{\mathrm{inc}}^{\mathrm{vs}}$ and $P_{\mathrm{inc}}^{\mathrm{et}}$ approximate errors are ±0.036 and ±0.022 GPa, respectively.

Inclusion pressures calculated using the volume strain method assume that inclusion strains are isotropic, and thus $P_{\mathrm{inc}}^{\mathrm{vs}}$ does not account for quartz inclusion differential strain. In the absence of a published range of acceptable differential strain values that approximate isotropic strains (i.e., ε₁–ε₃=0), inclusions with differential strain values between −0.001 and 0.001 were deemed as having approximately isotropic strain values (e.g., Gonzalez et al., 2019). Isolating the quartz inclusions whose strains are approximately isotropic significantly reduces the mean $P_{\mathrm{inc}}^{\mathrm{vs}}$ such that the difference between the volume strain and elastic tensor residual pressures exceeds P_inc error (Table 4). Given the disagreement between $P_{\mathrm{inc}}^{\mathrm{vs}}$ and $P_{\mathrm{inc}}^{\mathrm{et}}$ estimates of inclusion pressure, $P_{\mathrm{inc}}^{\mathrm{et}}$ estimates of inclusion pressure in the samples of interest are retained for entrapment isomeke calculations. Average $P_{\mathrm{inc}}^{\mathrm{et}}$ values are 0.066 ± 0.022 GPa and 0.067 ± 0.022 GPa for migmatite and iron formation samples, respectively (e.g., Table 4).

4.4 Ti-based chemical thermobarometry

Quartz titanium values varied considerably between inclusion and matrix grains. Specifically, Ti content generally ranged from below the detection limit (i.e., 1–2 ppm) to 50 and to 120 ppm for matrix and inclusion quartz grains, respectively (Fig. 10). Broadly speaking, inclusion quartz Ti values are commonly higher than matrix quartz Ti values. However, this trend of quartz inclusion Ti concentrations being higher than quartz matrix Ti concentrations is not consistent.

Figure 10 Quartz Ti concentrations in parts per million among matrix and inclusion grains measured in migmatite samples (a) and the iron formation sample (b). In many cases, quartz inclusion Ti values exceed matrix quartz Ti values. However, this trend is not consistent, and, in many cases, matrix and inclusion Ti values are similar. Error in Ti concentration measurements is 1–2 ppm.

Biotite inclusion Ti concentrations and phlogopite content were generally similar to their matrix counterparts in two out of the three samples in which inclusion biotite chemistry was available (Fig. 11). In migmatite samples, HP 82-88b, inclusion biotite grains were generally found to be more Ti-rich and magnesian relative to matrix biotite grains.

Figure 11 Matrix and inclusion biotite chemistry for migmatite samples (a) and the iron formation sample (b) shown relative to TiB thermometry isotherms (dashed and solid lines).

In general, quartz grains are texturally featureless in EM-CL imagery. A common exception includes sub-grain boundaries that are typically restricted to matrix grains although sub-grain boundaries are apparent in larger quartz inclusions (Fig. 12). Inclusion and matrix quartz grains may show slight core-to-rim CL color zoning (Fig. 12c and f), but most quartz grains imaged do not display noticeable zoning in CL nor do corresponding intragrain EMP analytical traverses. Average quartz inclusion Ti content was largely found to vary independently with distance from the garnet core. Lastly, Ti content of the garnet hosts themselves did not vary systematically from core to rim within the detection limit regardless of whether inclusions were concentrated in porphyroblast interiors (Fig. 13).

Figure 12 Representative blue-light EM-CL images of matrix quartz grains (a–c) and inclusion quartz grains (d–f). Subgrain boundaries are the most common textural feature observed in EM-CL imagery, most prevalently in matrix grains although subgrain boundaries may be also found in larger inclusion grains. While zoning textures are typically absent, slight color zoning can be observed in select matrix and inclusion quartz grains as in (c) and (f).

Figure 13 Backscatter electron photomicrograph of migmatite sample HP 82-64 showing a garnet host with an inclusion-rich core and an inclusion-poor mantle and rim (a). Average quartz inclusion Ti content and garnet EMP analytical spots along the compositional traverse are marked. Garnet Ti content in parts per million along the compositional traverses indicated on the backscatter electron image: T₁, T₂, T₃, and the “rim-to-core” traverse (b–e). Distance along compositional traverses is shown in micrometers (µm). The grey box superimposed on each compositional traverse indicates the garnet Ti detection limit of ∼ 25 ppm. Ti content of garnet hosts does not show consistent spatial trends either with distance between individual inclusions as in T_1–3 (b–d) or from the inclusion-rich interior to the inclusion-poor exterior as in the compositional traverse (e).

Ti-based chemical thermometry of inclusions serves as the primary constraint on prograde entrapment temperatures. TitaniQ inclusion temperatures are identified by the region of P–T space where TitaniQ isopleths intersect QuiG isomekes, whereas TiB inclusion temperatures are calculated at an assumed pressure of ∼ 0.40–0.60 GPa. In each sample, TiB inclusion thermometry yields higher temperatures than TitaniQ inclusion thermometry (this is also the case in the matrix dataset addressed in the next section). Typical inclusion temperatures calculated from TiB thermometry in samples of interest range from ∼ 665–780 ^∘C as opposed to ∼ 480–580 ^∘C for TitaniQ (Fig. 14).

Figure 14 QuiG entrapment isomekes (blue lines) calculated from average elastic tensor P_inc values shown relative to the average inclusion TitaniQ isopleth and inclusion TiB temperature range calculated for each sample. Entrapment isomeke thickness is representative of error. In all samples, QuiG–TiB entrapment estimates are found to be predicted at consistently higher P–T than QuiG–TitaniQ entrapment estimates.

Garnet–biotite thermometry temperatures were also calculated from inclusion biotite chemistry, and these temperature values were generally found to fall between inclusion thermometry estimates from TitaniQ and TiB. However, garnet–biotite temperatures were not considered further for prograde entrapment estimates due the prevalence of partially reset temperatures owing to Fe–Mg exchange between biotite and garnet (e.g., Spear and Florence, 1992). While post-entrapment garnet–biotite Fe–Mg exchange can impact the Ti content of biotite (e.g., Ti decreases as biotites become more magnesian upon heating), the impact on temperature is comparatively small due to the shallow nature of the biotite Ti-saturation surface isotherms (Henry et al., 2005).

4.5 Thermodynamic modeling and matrix chemical thermobarometry

Results of pseudosection modeling for representative samples HP 82-88b (migmatitic diatexite) and HP 81-82 (non-banded iron formation), with and without modal magnetite, are shown in Figs. 15 and 16, respectively. While the migmatite sample has a relatively well-constrained peak assemblage P–T range of ∼ 0.60–0.68 GPa and 720–800 ^∘C, the iron formation peak assemblage has a considerably wider predicted P–T stability regardless of whether magnetite is included (i.e., ∼ 675–800 ^∘C over the entire pressure range modeled). The large stability range predicted for the iron formation peak assemblage is due to the predominantly anhydrous Fe–Mg mineral paragenesis characteristic of the sample.

Figure 15 Pseudosection for migmatite sample HP 82-88b calculated with Domino in the MnNCKFMASHTO system. The peak assemblage is marked in blue, and the solidus is noted in green. Mineral assemblages of the remaining stability fields are marked. The inferred peak assemblage is modeled to be stable at ∼ 0.60–0.68 GPa and 720–800 ^∘C.

Mineral modal estimates are not predicted by Domino to change significantly within the predicted P–T range of peak assemblage stability, prohibiting the use of modal isopleths of the constituent minerals to further narrow in on equilibration conditions in the iron formation sample. Therefore, in order to further constrain peak P–T estimates, results from select matrix chemical thermobarometers are superimposed on calculated iron formation pseudosections (full inclusion and matrix thermobarometry results for HP 81-82 in addition to results from the other three samples are given in Table 5). Doing so using garnet barometry in conjunction with TiB thermometry restricts the P–T range of peak metamorphism for sample HP 81-82 to ∼ 700–735 ^∘C and ∼ 0.50–0.52 GPa. A lack of plagioclase in sample HP 81-82 inhibits further chemistry-based metamorphic pressure estimates. However, it is noted that the ∼ 0.50–0.52 GPa pressure estimate from garnet barometry in this sample is similar to the calculated pressure range from garnet–amphibole–plagioclase barometry (i.e., 5.37–5.83 GPa from 700 to 800 ^∘C) in an adjacent plagioclase-bearing iron formation from the study area (Daigle, 2015). It is also noted that although ilmenite is a stable matrix and inclusion phase in HP 81-82, pseudosection calculations with modal magnetite do not predict ilmenite to be stable. This lack of predicted ilmenite is likely due to the high amount of ferric iron as well as the small amount of Ti in the reconstructed bulk composition.

Figure 16 Pseudosections for iron formation sample HP 81-82 calculated with Domino in the MnCKFMASHTO system, both with modal magnetite (a) and without modal magnetite (b). The peak assemblage is marked in blue, and the solidus is noted in green. Mineral assemblages of the remaining stability fields are marked. Regardless of whether magnetite is considered as an inert phase, pseudosection modeling calculates large stability fields for the inferred peak assemblage (i.e., ∼ 675–800 ^∘C over the full range of pressures modeled). Matrix-thermobarometry-based P–T constraints (i.e., garnet barometry and TiB thermometry) are boxed in darker blue and labeled to further narrow in on equilibration conditions within the modeled peak assemblage stability field.

Table 5 Inclusion and matrix thermobarometry results for representative samples examined. Inc denotes inclusion-chemistry-based values, Mat denotes matrix-chemistry-based values, and Trav denotes garnet traverse-chemistry-based values. TitaniQ inclusion temperatures correspond to the temperature range over which the average TitaniQ isopleth intersects QuiG isomekes, whereas TitaniQ matrix temperatures are calculated at an assumed pressure of 0.60 GPa. TiB thermometry error varies specifically because it steadily decreases going towards higher temperatures (i.e., the error decreases from ± 23–24 ^∘C at 480–700 ^∘C down to ± 12 ^∘C at 700–800 ^∘C). It is difficult to give discrete values for QuiG barometry error given that it is the combined effect of calculated isomeke error (i.e., <0.1 GPa) and the paired thermometer (i.e., ± 20 ^∘C for TitaniQ and ± 12–24 ^∘C for TiB).

While peak conditions are relatively well-constrained using pseudosection modeling for migmatite sample HP 82-88b, garnet barometry and TiB thermometry P–T constraints from garnet interiors and biotite inclusions, respectively, are in good agreement with pseudosection peak metamorphic P–T estimates. Matrix mineral barometry for sample HP 82-88b on the other hand is below pressure constraints from inclusion/garnet traverse chemistry and pseudosection modeling (i.e., ∼ 0.45–0.50 GPa as opposed to ∼ 0.55–0.70 GPa). Thermobarometry results from the two remaining migmatite samples are like those calculated for HP 82-88b and HP 81-82. However, select calculated chemical barometry pressures for sample QC 81-20, an aluminosilicate-absent diatexite, are estimated to have been notably higher than the ∼ 0.50–0.70 GPa range from the two samples modeled with Domino (i.e., up to ∼ 0.90–1.10 GPa).

Identifying the approximate geologic P–T conditions of garnet nucleation and subsequent growth is a potentially useful means for bolstering inclusion thermobarometry findings. Predicted garnet modes for the two samples modeled in P–T space were found to vary primarily as a function of temperature (Figs. 17 and 18). Given that the peak assemblages of both samples in their respective pseudosections were found to be at temperatures at which significant dehydration of hydrous phases had already occurred, the addition of extra H₂O to garnet modal isopleth modeling efforts was necessary. Prediction of garnet modes for a fixed ${\mathrm{Fe}}^{\mathrm{3}+}/{\mathrm{Fe}}_{\mathrm{total}}$ ratio did not vary with increasing water content above 12.5 mol of additional water; thus the results from an assumed 25.0 mol of additional water are shown. In aggregate, the garnet-in temperature is predicted by Domino to have initiated over a temperature range of ∼ 400–600 ^∘C. This temperature range is consistent with the range calculated from Domino-based MnNCKFMASHTO garnet modal isopleth models of similar bulk compositions in the study area from an unpublished XRF dataset.

Figure 17 Garnet modal isopleth diagram for migmatite sample HP 82-88b calculated with Domino in the MnNCKFMASHTO system shown relative to QuiG–TitaniQ and QuiG–TiB inclusion thermobarometry results. Numbers given next to each isopleth represent the modeled garnet mode; N represents the first appearance of garnet. QuiG–TitaniQ results fall along the modeled garnet-in reaction, whereas QuiG–TiB results are at P–T conditions at which significant garnet modes are already predicted to be stable.

Figure 18 Garnet modal isopleth diagrams for iron formation sample HP 81-82 calculated with Domino in the MnCKFMASHTO system, both with modal magnetite (a) and without modal magnetite (b) shown relative to QuiG–TitaniQ and QuiG–TiB inclusion thermobarometry results. Numbers given next to each isopleth represents the modeled garnet mode; N represents the first appearance of garnet. QuiG–TitaniQ results are above the modeled garnet-in reaction yet below the range of significant garnet growth predicted ∼ 560–600 ^∘C. In both cases, inclusion thermobarometry constraints fall above the garnet-in reaction but below a temperature range of significant garnet growth predicted around ∼ 550–600 ^∘C. As in HP 82-88b, QuiG–TiB results are at P–T conditions at which significant garnet modes are already predicted to be stable.

Modeled garnet-in temperatures for both samples agree well with the temperatures of intersection between the average $P_{\mathrm{inc}}^{\mathrm{et}}$ QuiG isomeke and the average quartz inclusion TitaniQ isopleth. In sample HP 82-88b, the predicted inclusion entrapment conditions are found to lie along the modeled garnet-in reaction, approximately 20–30 ^∘C below the temperature range of a predicted increase in the garnet mode from 5 modal percent to 15 modal percent. The garnet-in reaction for HP 81-82 is predicted to occur at ∼ 400 ^∘C with modal magnetite versus at ∼ 450 ^∘C without modal magnetite. In both cases, inclusion entrapment conditions are predicted to lie at temperatures above the modeled garnet-in reaction and about 50–100 ^∘C below a significant T range of garnet growth predicted to occur from ∼ 550 to 600 ^∘C. Entrapment conditions predicted by the range of calculated TiB temperatures are at significantly higher temperatures than the modeled temperatures of significant garnet growth. Moreover, all garnet modal isopleth diagrams, regardless of the amount of ferric iron or water content, predict a significant amount of garnet to be stable at the prograde conditions of inclusion entrapment predicted by QuiG–TiB inclusion thermobarometry.

5 Discussion

5.1 Prograde entrapment P–T conditions from inclusion thermobarometry

Inclusion thermobarometry offers two P–T ranges of inclusion entrapment during prograde garnet growth depending on which inclusion thermometer is paired with QuiG elastic barometry. Among the samples analyzed, QuiG–TitaniQ inclusion thermobarometry estimates prograde conditions of interest to be approximately 0.55–0.70 GPa and 475–580 ^∘C, whereas QuiG–TiB inclusion thermobarometry constraints on prograde P–T conditions are comparatively higher: 0.85–1.10 GPa and 665–780 ^∘C. The lack of an observed difference in calculated inclusion P_inc values from garnet cores to garnet rims as well as approximately normal distributions in both sets of P_inc values suggests that the majority of inclusions studied were entrapped during a single garnet growth event. Inclusion entrapment conditions during this garnet growth event were likely the same among samples studied given that average $P_{\mathrm{inc}}^{\mathrm{et}}$ values were essentially identical between the representative migmatite and iron formation samples (i.e., 0.066±0.022 and 0.067±0.022 GPa, respectively). A single stage of inclusion entrapment is texturally supported by observations of inclusion-rich garnet cores and inclusion-deficient garnet mantles and rims in both lithologies studied. These lines of evidence strongly suggest that the two seemingly disparate views of prograde inclusion entrapment conditions offered by the inclusion thermobarometry pairs utilized are not both geologically meaningful in regard to prograde entrapment conditions.

It is interpreted that prograde conditions predicted by QuiG barometry and TitaniQ thermometry more accurately represent P–T conditions of inclusion entrapment and associated garnet growth. This interpretation is made primarily based on thermodynamic modeling results for the two samples modeled with Domino. The garnet-in reaction and significant subsequent garnet growth is predicted to have occurred at ∼ 400–650 ^∘C, a temperature range that better aligns with QuiG–TitaniQ entrapment estimates. A further piece of evidence in support of QuiG–TitaniQ inclusion-thermobarometry-based constraints is that a significant amount of garnet is predicted to already be stable above ∼ 600–650 ^∘C for all samples modeled, conflicting with the idea that TiB inclusion thermometry temperatures mark the main stage of inclusion entrapment via garnet growth. However, most of the quartz and biotite inclusions were almost certainly entrapped at the same relative P–T conditions during garnet growth. Thus, barring significant post-entrapment mechanical and/or chemical modifications to quartz and biotite inclusions, their respective thermobarometry calculations should ideally yield similar entrapment P–T predictions. The difference in inclusion thermometry estimates for TitaniQ versus TiB is likely due to quartz and biotite equilibration not occurring contemporaneously. Higher TiB inclusion temperatures relative to TitaniQ inclusion temperatures are most likely an expression of post-entrapment modifications of biotite Ti content in response to Fe–Mg exchange between garnet and biotite during prograde heating and, in part, retrograde cooling. The higher phlogopite content (i.e., increased Mg) of inclusion biotite grains relative to matrix biotite grains, most notably observed in HP 82-88b, is evidence for this Fe–Mg exchange behavior.

Despite the agreement between QuiG–TitaniQ inclusion efforts and thermodynamic models, it is worth further evaluating whether calculated pressures and temperatures reliably reflect prograde conditions of garnet growth instead of an amalgamation of post-entrapment chemical and mechanical modifications. Beginning first with inclusion thermometry constraints, the main question to assess is whether current Ti contents of quartz inclusions reflect their initial values upon entrapment. While Ti contents of biotite inclusions probably changed after entrapment due to Fe–Mg exchange with garnet hosts, it is interpreted that measured Ti values of quartz inclusions have not significantly departed from their initial entrapment values. Evidence for this interpretation lies in the lack of observable Ti variations in garnet porphyroblasts among the samples in which Ti in garnet was measured, suggesting that the Ti is not significantly partitioned into garnet. An additional key observation mandating against widespread Ti diffusion out of quartz inclusions into garnet hosts is the overall lack of quartz intragrain Ti zoning and the generally featureless nature of quartz inclusions in EM-CL imagery. If extensive quartz inclusion Ti movement and/or loss had occurred, then one would expect to see increasing CL intensity from the core of quartz inclusions towards the rim (i.e., the interface with garnet) (e.g., Spear and Wark, 2009; Spear et al., 2012). Such zoning was only weakly identified in one of the inclusion grains imaged (Fig. 12f); thus quartz inclusion Ti partitioning by garnet could have occurred, but it does not appear that this diffusion behavior served as a major influence on Ti content of quartz inclusions in the samples of interest. Of course, quartz could lose Ti via recrystallization at lower temperatures (e.g., Nachlas and Hirth, 2015). Ti loss in this manner likely resulted in calculated matrix TitaniQ temperatures being well below granulite facies conditions. However, microprobe analysis generally targeted encased quartz inclusions for Ti measurements that did not display visible recrystallization textures in EM-CL imagery, precluding post-entrapment inclusion Ti loss via recrystallization. From Ti profiles in garnet and EM-CL imagery of quartz inclusions, it is interpreted that the use of average TitaniQ results as a proxy for garnet growth conditions is a reliable representation of prograde entrapment temperatures. Use of the average TitaniQ results helps to mitigate the effects of underestimating entrapment temperatures due to inclusion Ti loss or perhaps even inclusions whose Ti contents represent pre-garnet conditions that did not re-equilibrate upon entrapment. Similarly, use of the average TitaniQ inclusion temperatures helps to account for areas of elevated Ti due to any diffusion or apparently higher Ti contents as a result of another phase within the sample volume with higher Ti content (e.g., ilmenite).

Similar to scrutinizing inclusion thermometry findings to discern whether calculated temperatures are meaningful for prograde path reconstructions, it is also worthwhile to discuss the geologic meaning of measured inclusion strains utilized in QuiG elastic barometry calculations. The main concern is that quartz inclusions have experienced post-entrapment modifications that have altered the original strains developed because of entrapment conditions, especially given the present inclusions' great antiquity as well as their preservation in granulite facies rocks. There are several modes by which inclusion strains might have been modified after entrapment, and it is considerably more difficult to confirm or dispute their occurrence relative to post-entrapment inclusion chemical modifications. Likely, the most prevalent mode of post-entrapment mechanical modifications is partial to complete strain release of quartz inclusions during exhumation and cooling as well as during thin section preparation. Observations that support some influence of this mode of post-entrapment mechanical modification on quartz inclusions include (1) the varying percentage of measured inclusions that display rightward shifts of Raman spectral peaks in each sample (i.e., ∼ 20 %–38 %, Table 3) and (2) the presence of some degree of fracturing in the majority of garnet hosts analyzed. Although Raman measurements intentionally targeted fully encased, spherical quartz inclusion and deliberately avoided inclusions proximal to visible fractures or embayments in garnet hosts, the key word is “visible”. Garnet host features such as fractures below the thin section surface could certainly lead to partial to complete strain relaxation in seemingly encased inclusions that are visibly located away from fractures or edges of the garnet host. Additional post-entrapment inclusion strain modifications that may have affected samples in this investigation include viscous creep and/or plastic yield of host garnet as well as shape maturation of inclusion quartz. The extent to which these factors may have affected host–inclusion dynamics during heating to peak temperatures in this study is particularly difficult to determine. However, it has been shown that these processes are more prevalent in rocks metamorphosed to granulite facies conditions (i.e., peak T's in excess of ∼ 700–750 ^∘C) (Zhong et al., 2020; Cesare et al., 2021).

As it is, there are several processes by which quartz inclusion strains could have been modified since they were entrapped in their garnet hosts during the Archean. These processes would ultimately result in either inclusions that are no longer strained or inclusions whose current strain states no longer reflect initial entrapment conditions. While it is impossible to fully characterize how one or more of the above post-entrapment mechanical processes ultimately impacted calculated P_inc values, the data suggest that these ancient inclusion strains are still largely reflective of their original entrapment conditions. This assertion is likely valid given both the inter-sample consistency of mean $P_{\mathrm{inc}}^{\mathrm{et}}$ values and the agreement of QuiG–TitaniQ entrapment conditions with thermodynamic modeling results. Additionally, inclusion/garnet interior chemical barometry constraints generally agree with the 0.55–0.70 GPa range of inclusion entrapment pressures indicated by QuiG isomeke intersections with TitaniQ isopleths (Table 5). Agreement between inclusion thermobarometry and thermodynamic models is further illustrated by the fact that average QuiG–TitaniQ entrapment conditions plot above the garnet-in reaction calculated by Domino (Figs. 17 and 18). It is worth acknowledging that the average inclusion entrapment temperatures indicated by inclusion thermobarometry are below the predicted start of significant garnet growth by ∼ 25 and ∼ 50 ^∘C in garnet modal isopleth models for the migmatite and iron formation, respectively. While this slight temperature mismatch could reflect post-entrapment strain modifications, this cannot be determined, because this temperature difference may also be attributable to the combined effect of several other factors (e.g., post-entrapment chemical modification of inclusions, instrument or thermobarometer error, calculated bulk composition, solution models used in thermodynamic models).

5.2 Peak metamorphic P–T conditions from mineral chemistry

Constraints on peak metamorphism from both matrix chemical thermobarometry and pseudosection modeling paint a consistent picture of peak metamorphic conditions during M1 granulite facies metamorphism of ∼ 0.50–0.70 GPa and ∼ 700–800 ^∘C. This M1 P–T range is consistent with determinations of M1 granulite facies metamorphism by previous researchers (Henry et al., 1982; Maas, 2004; Will, 2013; Daigle, 2015; Guevara et al., 2017). It is probable that minimum M1 temperatures were on the higher end of the 700–800 ^∘C range calculated in this study, likely above ∼ 730–750 ^∘C. This preference for the higher end of the calculated M1 temperature range is due to the highest TiB thermometry results from among the four representative samples (Table 5). The majority of TiB temperatures that are above 730 ^∘C are from inclusion grains of biotite in garnet and not matrix biotite grains, further demonstrating that, in many cases, biotite inclusions in garnet preserve temperatures at or near the terminus of the prograde path as opposed to temperatures of initial inclusion entrapment.

As discussed previously, rocks in the study area were subjected to at least one lower-grade overprinting event, principally a regional amphibolite facies event (M2) related to the emplacement of LLMC-type magmatic rocks. Additionally, there is evidence of a second middle to upper greenschist facies overprint (M3) related to subsolidus infiltration of fluids. However, the metasupracrustal rocks studied in this investigation do not commonly preserve textural evidence of overprinting events. Pinitization of cordierite grains and sericitization of plagioclase grains noted in select migmatite samples and the presence amphibole–iron carbonate clusters along iron formation fractures likely occurred at lower temperatures. Rather, the main samples that most commonly preserve widespread evidence of M2 overprinting are the immediately adjacent xenoliths of metaigneous mafic granulite and ultramafic lithologies. Widespread evidence of M2 overprinting in mafic granulites, for example, is primarily indicated by partial to complete hydration of pyroxene grains to amphibole. Thermobarometry calculations from these transitional granulite–amphibolite samples suggest that M2 overprinting occurred at approximately 0.50–0.60 GPa and 650–700 ^∘C. The pressures and temperatures calculated from matrix mineral chemical thermobarometry shown in Table 5 that are lower than M1 conditions commonly fall within this M2 P–T range and may certainly reflect M2 (or even M3) overprinting. However, it is difficult to make this determination in the absence of M2 textural evidence. It is likely that some of the recrystallization textures in EM-CL imagery of matrix quartz grains may reflect M2 deformation. Preservation of M2 recrystallization temperatures may be indicated by select matrix TitaniQ temperatures approaching ∼ 580–600 ^∘C. However, the majority of matrix TitaniQ temperatures are ∼ 400–550 ^∘C, temperatures that are more relevant for deformation during the waning stages of the Archean and into the Paleoproterozoic (e.g., Rowan and Mueller, 1971; Mueller, 1979).

Figure 19 Inferred P–T paths (black arrows) from representative samples relative to (1) average inclusion entrapment conditions (maroon boxes) indicated by QuiG isomekes (blue solid lines) and TitaniQ isopleths (dashed black lines), (2) M1 granulite facies peak conditions (green boxes) indicated by matrix thermobarometry and/or pseudosection results, (3) M2 amphibolite facies overprint conditions (gold boxes), and (4) the approximate location of the “wet” granite solidus (black solid line) based on the work of Luth et al. (1964). In each sample, a clockwise P–T path is indicated with largely isobaric heating to peak conditions along the prograde portions of the inferred P–T path.

5.3 P–T path implications from inclusion studies

The principal objectives of the present study are to determine the prograde metamorphic conditions of garnet growth in the eastern Beartooth Mountains and to use these prograde conditions to place fundamental constraints on the overall nature of metamorphic P–T paths in the study area. QuiG–TitaniQ inclusion thermobarometry paired with thermodynamic modeling results indicates a general inclusion entrapment P–T range of 0.55-0.70 GPa and 475–580 ^∘C. Pairing inclusion findings with the estimated range of M1 granulite facies conditions from matrix chemical thermobarometry and pseudosection modeling (i.e., ∼ 0.50–0.70 GPa and ∼ 730–800 ^∘C), it is interpreted that lithologies of interest were heated to peak temperatures isobarically at approximately ∼ 0.50–0.70 GPa. Following the M1 event, samples of interest later proceeded along a retrograde P–T path that included at least one overprinting event, likely the regional M2 amphibolite facies overprint whose conditions were ∼ 0.50–0.60 GPa and 650–700 ^∘C. Placing each of these metamorphic steps into a continuous framework, an overall “flat”, clockwise P–T path emerges for the high-grade metasupracrustal rocks studied (Fig. 19), a P–T path similar to the path predicted by Guevara et al. (2017). Although this study does not advocate for the prograde constraints calculated from QuiG–TiB inclusion thermobarometry, it is worth noting that QuiG–TiB constraints also indicate an overall clockwise metamorphic P–T path. However, the higher entrapment temperature estimates from TiB inclusion thermometry require significant isothermal decompression near peak temperatures that is not supported by the data.

A clockwise P–T path that includes approximately isobaric heating to peak temperatures is consistent with the late Archean (i.e., ∼ 2.74–2.83 Ga) continental arc subduction zone setting inferred by previous researchers. The P–T path modeled here implies that the passive-margin pelitic to psammitic protoliths of the studied migmatites and iron formations were buried at ∼ 25–30 km depth and tectonically mixed via collisional tectonics and thickening where they were then heated to granulite facies temperatures, likely due to intruding arc plutons. Intrusion of subduction-derived magmatic bodies served as the primary means of overprinting the granulite facies peak assemblages. While the metasupracrustal rocks studied likely experienced regional amphibolite facies overprinting (M2) and perhaps even local granulite facies fluid conditions all due to intrusion of igneous arc sills and dikes, these later overprints are not easily resolvable in the present samples. Upon cessation of metamorphism, it is inferred that the rocks followed a path of cooling and decompression that produced many of the low-temperature alteration phases observed (e.g., pinite after cordierite, sericite after feldspar).

6 Summary and conclusions

Prograde metamorphic conditions of garnet growth in the eastern Beartooth Mountains of Montana are estimated to be 0.55–0.70 GPa and 475–580 ^∘C. Quartz inclusions were likely entrapped during a single stage of garnet growth along the prograde path of metamorphism. Matrix chemical thermobarometry and pseudosection modeling indicate that, after entrapment, samples were heated isobarically to peak granulite facies conditions of ∼ 0.50–0.70 GPa and 730–800 ^∘C, followed by at least one amphibolite facies overprint (∼ 0.50–0.60 GPa and 650–700 ^∘C) along the retrograde path. Overall, the samples of interest experienced a clockwise P–T path, consistent with metamorphism in a convergent margin setting such as a continental arc subduction zone previously demonstrated to be the operative tectonic regime in the late Archean when samples of interest were metamorphosed. Despite likely experiencing some form of post-entrapment mechanical and/or chemical modifications, it appears that the subset of quartz inclusions used to make prograde entrapment estimates still largely records their original entrapment conditions despite their great age (i.e., ∼ 2.80 Ga). Additionally, although there may have been slight deviations from the inferred P–T path, similarities between calculated entrapment and peak conditions from samples examined suggest that all metasupracrustal samples studied likely experienced the same general P–T path described above. This consistency in geologic history is noted despite the interruption of lithologic continuity among xenoliths throughout the study area by the surrounding volumetrically dominant ∼ 2.80 Ga TTG plutonic rocks.

The results of this investigation are geologically significant as they serve to clarify a portion of the metamorphic history in the eastern Beartooth Mountains that has not been previously well-constrained. Furthermore, inclusion constraints described here serve to demonstrate the clockwise nature of the overall metamorphic P–T path in the metasedimentary garnet-bearing rocks of interest. Proof of a clockwise metamorphic P–T path lends an additional line of evidence to the ever-growing body of work strongly indicating that late-Archean-aged rocks throughout the Beartooth range and the wider Wyoming Province were produced via modern-style collisional plate tectonic processes. The present work is important to the field of inclusion studies, specifically inclusion elastic thermobarometry. Finally, this work demonstrates that quartz-in-garnet barometry may be meaningfully applied to ancient rocks (among the oldest rocks used for inclusion elastic thermobarometry) if elastic thermobarometry model assumptions are met.