The bedrock exposed in the Nashoba terrane of east-central Massachusetts records a complex history of deformation and metamorphism associated with the collision of Avalonia with the Laurentian margin during the Devonian Acadian orogeny. Although the structural history of the terrane has been well studied, its pressure-temperature (P–T) history is less well constrained, and the mechanisms by which the upper amphibolite facies Nashoba terrane was juxtaposed between greenschist facies rocks to the southeast and northwest have not been established. Here, we apply phase equilibria modeling, geothermobarometry, and petrographic analysis to three garnet-bearing migmatitic rocks from the Nashoba terrane to determine their P–T histories and provide key metamorphic constraints relevant to Acadian orogenic processes. All three samples are from the Nashoba Formation, a unit within the Nashoba terrane consisting of volcaniclastic rocks likely deposited in an arc/back-arc setting near the trailing edge of the Ganderia terrane. Peak subsolidus conditions are determined via the integration of petrographic analysis and thermodynamic modeling. Peak anatectic conditions are constrained with garnet-biotite thermometry + GASP barometry, garnet intersecting rim isopleths, and thermodynamic modeling.

Phase equilibria constraints suggest peak subsolidus conditions ranging from ~550°C to 700°C and ~6–12 kbar. Geothermobarometry, phase equilibria modeling, and garnet rim isopleths constrain biotite-out anatexis to ~700°C–715°C and ~5–8.5 kbar, up to ~2–4 kbar deeper than previously suggested. This synthesis of phase equilibria modeling and petrographic analysis suggests that all three samples record a clockwise P–T path with peak pressures achieved before anatexis associated with the Acadian orogeny at ~394 Ma. These results are inconsistent with a previously interpreted pre-Acadian period of low P/T metamorphism (<5 kbar peak pressures). Instead, this implies a previously unrecognized phase of intermediate P/T crustal thickening preceding anatexis, which we interpret as a result of Acadian orogenesis. We suggest that the Nashoba terrane exposes mid-crustal levels of the Acadian hinterland.

Central and southern New England are composed of several peri-Laurentian and Gondwana-derived terranes that together record a long history of continuous accretion and collision onto the Laurentian margin throughout much of the Paleozoic, resulting in a record of complex metamorphism and deformation (Figure 1) [1-5]. Constraining the metamorphic evolution of these terranes is key to developing tectonic models of Laurentia and evaluating paleotectonic reconstructions. The Nashoba-Putnam terrane in east-central Massachusetts and eastern Connecticut consists of high-grade metamorphosed volcanic, volcanoclastic, and pelitic rocks thought to record metamorphism and deformation associated with ~100 Myr of collision on the eastern margin of Laurentia [3, 5]. The Nashoba-Putnam terrane has been interpreted as a remnant volcanic arc complex at the trailing edge of the Peri-Gondwanan arc-back-arc terrane Ganderia [3-11]. Ganderia experienced a protracted Paleozoic tectonic history including (1) the development of an arc on the Gondwanan margin and subsequent rifting of the arc from Ganderia from ~540 to 500 Ma, (2) collision of the Ganderian leading-edge with Laurentia instigating the ~430 Ma Salinic orogeny, (3) ~420 Ma collision of the Avalonian terrane with the Ganderian trailing edge leading to the ~400 Ma Acadian orogeny, and (4) ending with the ~355 Ma collision of the Meguma terrane onto the composite Laurentian-Ganderian margin and subsequent exhumation [3-5]. Thus, the rocks of the Nashoba terrane may record metamorphism associated with a substantial portion of Paleozoic tectonism.

Although the geochemical affinities [4, 5, 10], structural geology [11-13], and geochronology [3, 5, 7, 10, 14, 15] of the Nashoba terrane have been previously documented, data on the complex pressure-temperature-time-deformation (P–T–t–d) history of these rocks are limited [7]. Hepburn et al. [3] and Walsh et al. [5] have synthesized geochronologic, structural, stratigraphic, geochemical, and petrologic data to develop tectonic models for the formation of the Nashoba-Putnam terrane. Both studies suggest that the bulk of metamorphism recorded in the Nashoba terrane is low-P/T, amphibolite facies arc-related regional metamorphism as a result of late Silurian westward-dipping subduction of the Acadian Seaway under the trailing edge of Ganderia at ~425 Ma, and widespread anatexis during subsequent collision of Avalon during the Devonian Acadian orogeny at ~400 Ma. The Acadian collision was followed by greenschist facies deformation and fluid infiltration as a result of the docking of Meguma, exposed in Nova Scotia, Canada, on the eastern margin of Avalonia during the Neo-Acadian orogeny at ~360–320 Ma [7, 8]. There is also some evidence (explored below) for metamorphism predating the Acadian orogeny, potentially linked to the collision of the Ganderian leading edge with Laurentia and associated closure of basins within Ganderia at ~450–430 Ma preceding the Salinic orogeny [3]. The Nashoba terrane experienced higher temperatures than the adjacent rocks of the Ganderian-derived Merrimack belt to the northwest and Avalonia to the southeast, but the peak pressures are not well constrained [11, 16, 17]. The inferred presence of kyanite based on the reported presence of pseudomorphs of sillimanite after kyanite in Nashoba rocks [16, 18] suggests an early period of intermediate-P/T metamorphism [19] that necessitates crustal thickening. This style of metamorphism is typical of continent-continent collision following the subduction of a major ocean basin [20] but can also be produced via back-arc closure and associated crustal thickening [21-25]. Although the metamorphic record in the Nashoba-Putnam terrane can potentially help elucidate the dynamics of the collision of the Avalon terrane with the composite Laurentian margin in Central New England, quantitative constraints on the P–T evolution of these rocks are lacking.

In this study, we combine thermodynamic modeling, geothermobarometry, and petrographic observations to produce a well-constrained metamorphic P–T history for three samples collected along a N-S transect in the Massachusetts portion of the Nashoba-Putnam terrane. These samples were selected roughly along strike in order to investigate samples believed to represent comparable metamorphic conditions, assuming limited post-emplacement displacement. These P–T histories are discussed within the context of the Nashoba terrane’s tectonic history and recently proposed tectonic models encompassing prograde metamorphism and exhumation.

Our field area within the Nashoba terrane occupies the ancestral, traditional, and contemporary unceded lands of the Nipmuc and Agawam tribes. We recognize that we operate on stolen land, and we extend our gratitude to the many Indigenous peoples who have rich histories here for their ongoing stewardship of the land. We honor and respect the enduring relationship between these peoples and this land, as well as the strength of Indigenous culture and knowledge. We commit to recognizing, supporting, and advocating for the many Indigenous peoples who live, work, and study at our respective institutions. By offering this Land Acknowledgment, we commit to holding ourselves and our institutions more accountable to the needs of Indigenous peoples.

3.1. Stratigraphic History

The portion of the Nashoba-Putnam terrane outcropping in east-central Massachusetts is known as the “Nashoba terrane” [2-7]. It consists of complexly deformed amphibolite-facies metamorphosed sedimentary, volcanic, and volcaniclastic rock [2, 12, 26, 27]. It is separated from the Merrimack belt in the west by the Clinton-Newbury fault zone (Figure 2). The Merrimack belt consists of metasedimentary rocks sourced largely from ~460 to 443 Ma Ganderian basement but with some Laurentian input [3, 28, 29]. The Merrimack belt was deposited at ~435–420 Ma, during or following suturing of the leading edge of Ganderia with the post-Taconic composite Laurentian margin [3, 4]. The Nashoba terrane is separated from Avalonia in the east by the Bloody Bluff fault zone (Figure 2; [2, 3, 5, 12, 30]). The structurally lowest units of the Nashoba terrane are the metavolcanic and metasedimentary rocks of the Marlboro Formation, Shawsheen (para)Gneiss, and Fishbrook (ortho)Gneiss. In this study, we investigate samples from the Nashoba Formation, which consists of high-grade metamorphic volcanogenic metasedimentary schists, gneisses, and migmatites that occupy about a third of the exposed terrane west of the Marlboro Formation (Figure 2) [2, 3, 27, 31, 32]. The Nashoba Formation lies stratigraphically above the Marlboro Formation, but the precise structural relationships between the formations are complex and overprinted by a history of ductile and late brittle deformation [5, 11-13].

3.2. Tectonic and Metamorphic History of the Nashoba Terrane

Kay et al. [4] and Hepburn et al. [3] argue that εNd values from metaigneous and U-Pb detrital zircon ages from metasedimentary rocks in the Nashoba terrane suggest a Ganderian affinity. Ganderia is a Gondwana-derived crustal fragment thought to have experienced three distinct cycles of arc/back-arc opening and closure, the first two of which are pertinent to the Nashoba terrane [3]. The first cycle resulted in the Penobscot arc/back-arc complex, which formed due to the subduction of the Iapetus Ocean under Ganderia as the Penobscot complex rifted from Gondwana with associated volcanism beginning at ~650 Ma and possibly continuing until the Cambrian [8, 33, 34]. The igneous protoliths of the Marlboro Formation metavolcanic rocks have been dated to ~500 Ma and are interpreted as rift volcanism products [3, 4]. Fragments of this rifted arc collided with each other and consolidated at ~486–478 Ma during the Penobscot orogeny [3, 8, 29, 35-38].

The 478–453 Ma Popelogan-Victoria arc/back-arc complex formed on the older Penobscot basement and associated rifting generated the Tetagouche-Exploits basin [8, 36-39]. Few magmatic zircons of Popelogan-Victoria age have been found in the Nashoba Formation. Instead, the bulk of Nashoba Formation detrital zircon ages represent the older Penobscot arc/back-arc volcanism at ~500–540 Ma, though the youngest detrital zircon age determined by CA-TIMS is ~466 Ma [3]. These ages suggest that while deposition of the Nashoba Formation metasedimentary protoliths into the Tetagouche-Exploits basin occurred primarily during the Ordovician Popelogan-Victoria cycle, the sediments were largely sourced from weathering of the Penobscot basement. The collision of the Ganderian leading edge with the eastern margin of Laurentia initiated the closure of the Tetagouche-Exploits basin within Ganderia (~450–435 Ma) generating what is called the Salinic orogeny or Salinic disturbance at ~430 Ma [5, 8, 36]. Basin closure and convergence likely persisted until the ~425–420 Ma initiation of the Acadian orogeny [8, 36, 39].

Volcanic activity associated with a west-dipping subduction zone between Ganderia and Avalonia generated dioritic and granodioritic intrusions and the earliest major amphibolite-grade metamorphic event in the Nashoba terrane (M1) [3, 5, 7]. Monazite U-Pb-Th microprobe dating from Stroud et al. [7] suggests that M1 occurred at 423 Ma with a ~30 Myr duration. Stroud et al. [7] argue that concentric overgrowths on metamorphic M1 monazite suggest this event was the result of deformation-free, static thermal metamorphism associated with subduction and the emplacement of I-type granitoids. M1 mineral assemblages are poorly preserved and represented by sillimanite-muscovite-biotite and sillimanite-garnet-biotite-muscovite assemblages in the Nashoba Formation, with peak conditions of greater than ~600°C and 3.2–4.8 kbar [7].

Stroud et al. [7] argue that peak metamorphic conditions in the Nashoba terrane occurred at ~394 Ma as a result of the accretion of the Avalonia terrane onto Laurentia during the Acadian orogeny (~400–385). The Acadian metamorphic event (M2) was responsible for widespread migmatization across the terrane. Metamorphic zircon in the Nashoba terrane indicates metamorphism lasting from ~440 to ~335 Ma [3, 5, 7]. M2 is thought to be characterized by sillimanite-K-feldspar-biotite and sillimanite-K-feldspar-biotite-garnet assemblages with widespread anatexis at conditions greater than 650°C and 2.5–4.5 kbar.

The Late Devonian accretion of the Gondwana-derived Meguma terrane with Laurentia in present-day Canada (the Neoacadian orogeny) imposed a minor, local metamorphic overprint on the Massachusetts portion of the Nashoba terrane, with recrystallized monazite rims indicating fluid infiltration from 372 to 305 Ma (M3 in [7]). Active convergent tectonics ended with the collision of Gondwana with Laurentia during the Late Mississippian Alleghanian orogeny, finishing the assembly of the Pangean supercontinent between ~325 and 260 Ma ([6] and references therein). Ganderian rocks exposed in Massachusetts experienced minimal low P–T overprinting from the Alleghanian orogeny, though there is evidence for more substantial overprinting in Southern New England [40-43].

Stroud et al. [7] investigated the petrochronology of the Nashoba terrane through the synthesis of monazite petrochronology and P–T estimates based on a comparison of regionally observed mineral assemblages with petrogenetic grids [44]. These results suggest a relatively low-P/T Buchan-style metamorphic history taking place entirely below ~5 kbar [45]. Abu-Moustafa and Skehan [46] present the earliest study of metamorphism within the Nashoba Formation and argue for peak conditions of ~650°C and 6 kbar based on pelitic mineral assemblages and early petrogenetic grids. Studies based on petrogenetic grids present broad constraints on metamorphic histories that can be refined with the addition of bulk composition-specific petrologic techniques such as phase equilibria modeling, the method of intersecting isopleths, and geothermobarometry. Bulk composition-specific techniques are also able to account for local variations in the effects of bulk composition on metamorphic phase equilibria [47], which can have large effects on predicted aluminosilicate stability [48] and the P–T conditions of melt reactions [49, 50]. As a result, the incorporation of sample-specific bulk compositions and thermodynamic modeling with petrographic analysis and other tools can uncover subtle aspects of P–T histories obscured by broader approaches. In this study, we investigate the presence of intermediate pressure assemblages within the overall framework of the Nashoba metamorphic history by presenting integrated petrographic analysis, geothermobarometry, and phase equilibria modeling for two kyanite or kyanite-pseudomorph bearing metapelitic samples and one sillimanite bearing sample along a N-S transect in the Nashoba Formation. We then relate this P–T history to the previously established regional tectonic framework in an effort to constrain the prograde evolution of the Nashoba terrane, as well as speculate on the mechanisms that emplaced the hot and high-grade Nashoba terrane between the cold and low-grade Avalonia and Merrimack belt.

4.1. Analytical Methods

Major element X-ray maps and mineral compositions for sample JWB13-8 were collected on the Cameca SX5 at the Syracuse University Electron Microprobe Laboratory. X-ray maps were collected with both wavelength-dispersive and energy-dispersive spectrometers. Maps were collected with a focused 100 nA beam current at a range of step sizes from 2 to 10 µm, and a count time of 60 seconds. Quantitative silicate analyses were conducted using a 20 nA beam with a 5 µm spot size for sheet silicates and 1 µm for other minerals. These analyses were calibrated using a variety of synthetic and natural standards loosely matched to the minerals analyzed, including well-characterized garnet and feldspar standards. This is the standard approach for the field, and readers are directed to the Syracuse University Electron Microprobe Laboratory for further details (Jay Thomas, Personal Communication). X-ray maps and mineral compositions for samples 21-AC-04 and 21-AC-06 were collected on the Cameca SX5-Tactis at the American Natural History Museum electron microprobe facility. The same analytical settings described above were employed.

Accurately constraining the bulk composition for subsolidus modeling of a migmatitic rock can be complicated, as the original bulk composition may be modified by melt segregation or migration syn- to post-anatexis [51]. In order to mitigate these effects, the suggested methods of Guevara and Caddick [51] were followed. Chiefly, representative hand samples were chosen to most accurately reflect the modal compositions of leucosome and melanosome at the outcrop scale. Additionally, the presence of muscovite and biotite in the thin section and patchy migmatite outcrop texture suggests minimal melt loss (online Supplementary Figure S1). The bulk compositions of the representative hand samples were collected via fused glass XRF on the Thermo ARL Peform’X XRF at the Hamilton Analytical Lab. The precision and accuracy of this method are on the order of <0.2 wt % or lower for major elements. Please see the official Hamilton Analytical Lab materials for further details.

4.2. Thermodynamic Modeling and Geothermobarometry

Mineral assemblage diagrams (MADs) are isochemical phase diagrams that display the stable mineral assemblages for a specific bulk composition over a range of P–T space. MADs constructed for this study were calculated using a module in the FORTRAN program GIBBS [52-54]. Program GIBBS is capable of a range of thermodynamic calculations, including a Gibbs energy minimization routine for the calculation of MADs and modules for modeling progressive metamorphism including phase fractionation, at equilibrium or overstepped conditions. Assuming the same bulk composition and thermodynamic dataset, Program GIBBS produces a MAD equivalent to what is produced by Theriak-Domino [54]. See Spear and Wolfe [54] for more details on the methodology of the program.

In this study, MADs were calculated in the MnNCKFMASHTi chemical system employing both the SPaC18 dataset and the Holland and Powell [55] ds6.2 (referred to as HP11 herein). Our HP11 modeling employed the thermodynamic data presented in White et al. [56], with the exception of plagioclase/K-feldspar, epidote, and ilmenite. These data are from Holland and Powell [57], Holland and Powell [55], and White et al., [58], respectively. The SPaC dataset is based on Berman [59] and employs ideal ionic solution models for all minerals except for garnet, which is nonideal [60]. Full details of the thermodynamic data and solution models employed can be found in Spear and Wolfe [54].

P–T histories were constructed by synthesizing textural observations with thermodynamic modeling and classical geothermobarometry. Geothermobarometric estimates were calculated with the garnet-biotite thermometer of Hodges and Spear [61] updated with the Berman [60] garnet model. Pressure estimates were obtained using the garnet-aluminosilicate-plagioclase barometer from Hodges and Crowley [62]. We adopt the minimum errors reported for each calibration, ~50°C and ~1 kbar.

5.1. Overview

Samples were collected at roadcuts defining a ~10 km N-S transect through the central portions of the Nashoba Formation (Figure 2). Sample JWB13-8 is the northernmost sample and is an example of the Nashoba Formation gneissic unit, Snsg (Table 1) [27]. 21-AC-06 is approximately ~7 km south of JWB13-8 and from the same unit (Table 1). 21-AC-04 is ~4 km south of 21-AC-06 and is from the biotite gneiss and amphibolite (Sngba) unit [27]. While all three samples are not mapped as migmatite, they all show evidence for variable degrees of anatexis and/or melt loss on either the outcrop or thin section scale. This includes patchy to stromatic leucocratic melt segregations at the outcrop scale (21-AC-04 and 21-AC-06) and peritectic garnet growth at the thin section scale (JWB13-5). We employ the migmatite nomenclature of Sawyer [63] when describing these samples. Additionally, to our knowledge, 21-AC-04 preserves the first published relict prograde kyanite in the Nashoba Formation. A detailed description of the petrography of these samples follows, and a summary of sample stratigraphy, mineralogy, and location is provided in Table 1.

5.2. Petrography

5.2.1. JWB13-8

JWB13-8 is a garnet-plagioclase-quartz-kyanite-sillimanite-biotite gneiss, with accessory apatite, zircon, monazite, and ilmenite. The dominant foliation is defined by aligned biotite and alternating layers of intergrown biotite + sillimanite and quartz + plagioclase and quartz leucosomes (Figures 3(a) and 3(b)). Garnet grains are spatially aligned within the foliation, and continuous seams of sillimanite needles bend around garnet grains. Rare mm-scale sillimanite-after-kyanite pseudomorphs within the matrix are aligned with the fabric as noted by previous researchers [18]. Shear bands associated with an incipient S-C structure crosscut all fabric elements, including the sillimanite seams. Tabular biotite grains range from <0.5 to over 1 mm in length. Quartz, plagioclase, and quartz + plagioclase pods appear continuous in outcrop scale but are somewhat discontinuous in thin section scale and are occasionally micro-boudinaged. There is no muscovite present in any of the thin sections prepared from this sample.

Garnet grains in JWB13-8 occur in two habits. The dominant habit (G1) occurs as 1.5–2 mm diameter subidioblastic porphyroblasts. These grains typically contain inclusions of biotite, quartz, zircon, plagioclase, and ilmenite. Subidioblastic garnet grains appear pretectonic with respect to the matrix foliation and appear rotated (Figures 3(e) and 3(f)). Rare <10 µm diameter (sub)idioblastic inclusion-poor garnet grains (G2) occur in quartz + plagioclase leucosomes (Figures 3(c) and 3(d)). G1 grains are interpreted as grains that nucleated and at least partially grew at subsolidus conditions, while G2 is interpreted to be peritectic in origin.

5.2.2. 21-AC-06

21-AC-06 is a garnet-plagioclase-quartz-sillimanite-biotite gneiss. It also contains the accessory phases zircon, apatite, ilmenite, and monazite. The gneissic banding is defined by alternating layers of foliated biotite + prismatic sillimanite and quartz + plagioclase ribbons. Biotite grains define the planar foliation and are often wrapped around large plagioclase microboudinaged xenoblasts (Figures 4(a) and 4(b)). Sillimanite occurs exclusively as small prismatic intergrowths within matrix biotite.

Garnet occurs as ~1–2 mm in diameter subidioblastic to xenoblastic grains (Figures 4(c) and 4(d)). Larger garnet grains typically preserve polycrystalline inclusions of quartz + biotite with rare plagioclase. Garnet grains also contain inclusions of zircon, apatite, ilmenite, and aluminosilicate. Garnet grains are typically rimmed or embayed by biotite + plagioclase.

5.2.3. 21-AC-04

21-AC-04 is a garnet-plagioclase-quartz-kyanite-sillimanite-biotite-muscovite gneiss. It contains the accessory phases apatite, zircon, monazite, and ilmenite. The matrix foliation is defined by intergrown biotite + prismatic sillimanite + fibrolite with rare muscovite (Figures 5(a) and 5(b)). Rare kyanite occurs as <0.25 mm long xenoblastic grains rimmed by muscovite + chlorite (Figures 5(c) and 5(d)). Lenses of quartz, plagioclase, and quartz + plagioclase are discontinuous on the thin section scale but continuous at the outcrop scale. ~0.25 mm long sillimanite porphyroblasts randomly overgrow the foliation.

Garnet in 21-AC-04 is rare and occurs as xenoblastic ~0.8 mm diameter grains and <0.5 mm diameter subidioblastic grains (Figures 5(c) and 5(d)). Garnet grains in this sample are relatively inclusion-free compared with JWB13-8 and 21-AC-06 but contain rare inclusions of quartz and biotite. Cracks are filled with chlorite + biotite. Both garnet habits appear wrapped by the matrix foliation.

6.1. JWB13-8

Subsolidus G1 garnet in JWB13-8 is variably zoned with respect to major elements (Figures 6(a)–6(d)). Fe zoning is essentially flat, while Mg content decreases from core to rim (Figures 6(a) and 6(b)). Ca zoning displays a subtle decrease from the core until ~50 µm from the rim where there are sharp, discontinuous high-Ca rims where the grain is in contact with biotite (Figure 6(c)). Mn is unzoned, except for an uptick in Mn at the rims (Figure 6(d)). Zoning in all cases is subtle (<0.01 mol% differences core to rim; Table 2). These patterns are interpreted to reflect at least partial diffusive relaxation/resetting of subsolidus growth zoning at peak or near peak temperature conditions. High Mn rims are interpreted to reflect retrograde net-transfer reactions. Peritectic G2 garnet is unzoned, with a composition broadly similar to the rims of G1.

Matrix biotite grains are unzoned and are typically XAnn = 0.57. Biotite replacing garnet rims is slightly more magnesium-rich, with XAnn = 0.56 (Table 3). Plagioclase grains are also unzoned and are typically XAn = 0.34 (Table 4).

6.2. 21-AC-06

Garnet in 21-AC-06 is unzoned with respect to Fe, while Mg records a subtle core-to-rim decrease, with more dramatic decreases in embayments around biotite (Figures 6(e) and 6(f)). Ca is largely unzoned with the exception of a subtle increase in the concentration at the rim, along cracks, and around polymineralic quartz + biotite ± plagioclase inclusions (Figure 6g). Mn displays reverse zoning, with a particularly high concentration on the outer ~10–20 µm (Figure 6h). Like JWB13-8, the preserved zoning is interpreted to result from a combination of diffusional relaxation at peak/high T (Fe, Mg, and Ca) and retrograde consumption of garnet by biotite + plagioclase + quartz (Mn).

Biotite replacing garnet is typically ~XAnn = 0.54, while matrix biotite is typically ~XAnn = 0.56, and unzoned (Table 3). Plagioclase grains in the matrix and replacing garnet are of similar composition, typically XAn = 0.27 (Table 4).

6.3. 21-AC-04

Garnet in 21-AC-04 is largely unzoned in major elements (Figures 6(i)–6(l)), with the exception being high Mn zones along cracks filled with biotite + chlorite. Biotite along garnet rims and in the matrix is all the same composition of ~XAnn = 0.60 (Table 3). Plagioclase in 21-AC-04 is typically nearly pure albite, at XAn = 0.05 (Table 4). Taken together, the lack of garnet zoning and homogeneous biotite compositions suggests diffusive relaxation of garnet growth zoning and homogenization of biotite compositions at (near-) peak temperature conditions.

Here, we present MADs and geothermobarometry results for samples JWB13-8, 21-AC-06, and 21-AC-04. The bulk compositions employed for thermodynamic modeling can be found in Table 5. Approximately 5 wt% water was added to each bulk composition to ensure fluid-saturated conditions. The application of geothermobarometry and thermodynamic modeling of high-grade rocks can be tricky as retrograde net-transfer and exchange reactions can modify mineral compositions post-peak conditions [64]. To account for this complication, we employed the method outlined by Storm and Spear [65]. Chiefly, garnet-core-matrix biotite pairs were used to calculate temperature maxima, while garnet rim-adjacent biotite pairs were used to calculate temperature minima. Additionally, melt segregation, migration, and loss can modify the original bulk composition of a rock enough to render modeling of pre-anatectic phase equilibria impossible without melt reintegration [49-51]. To help alleviate these effects, we selected samples that best represent outcrop modal abundances of leucocratic and melanocratic material. The general agreement between geothermobarometry, garnet rim isopleths, thermodynamic modeling, and petrography (explored below) suggests that retrograde and melt segregation effects are fairly minimal.

Figure 7 presents a comparison of MADs calculated with the SPaC (a, c, and e) and the HP11 (b, d, and f) thermodynamic datasets for all three samples. Phase assemblages of peak subsolidus and anatectic conditions are also highlighted in the HP11 MADs. For each sample, a minimum 50°C and ~1 kbar error for garnet core and rim barometry is plotted as a gray parallelogram representing a minimum uncertainty envelope, after [66]. In each case, core and rim thermobarometry overlapped within the minimum analytical error. Individual core and rim barometry with associated uncertainties can be found in online Supplementary Figure S2. Garnet near-rim compositional isopleth intersections are indicated by the star. Highlighted fields show P–T ranges and phase assemblages of peak subsolidus conditions in yellow, conditions of anatexis in red, and conditions of retrograde metamorphism in blue.

7.1. HP11 vs SPaC

Although both HP11 and SPaC calculate similar assemblages for the peak subsolidus and anatexis fields, the pressure of these assemblages in HP11 is generally lower than the same assemblages in SPaC, though they still fall within a conservative ~±1 kbar error of each other [67]. Both the garnet core and rim isopleths do not intersect on any of the HP11 MADs. The proposed P–T path in Figure 8 is based on the SPaC thermodynamic dataset results due to the consistency between the SPaC dataset with the thermobarometry results and the petrographic analyses, though our general conclusions are consistent regardless of the dataset used.

7.1.1. JWB13-8

In sample JWB13-8, the lack of staurolite and the presence of subsolidus garnet and kyanite constrains the peak subsolidus assemblage to quartz + H2O + plagioclase + muscovite + biotite + garnet + ilmenite + kyanite, which is stable between 650°C and 700°C and 7–11 kbar. The melt-in and the muscovite-out reaction constrains the anatectic assemblage to quartz + H2O + plagioclase + biotite + garnet + melt + ilmenite + sillimanite, which is stable between 700°C and 750°C and 3.5–7 kbar. The final equilibrium assemblage of quartz + H2O + plagioclase + biotite + K-feldspar + ilmenite + sillimanite is stable between ~625°C and 700°C and ~3–5 kbar. Thermobarometry results for JWB13-8 range from ~600°C to 750 °C and ~2–7 kbar. This is broadly consistent with the rim isopleth conditions of ~715°C and 5.5 kbar.

7.1.2. 21-AC-06

In sample 21-AC-06, the absence of K-feldspar and the presence of subsolidus garnet constrains the peak subsolidus assemblage to quartz + H2O + plagioclase + muscovite + biotite + garnet + ilmenite, which is stable from 575°C to 675°C and 6–9 kbar. The melt-in and sillimanite-in and muscovite-out reactions constrain the anatectic assemblage to quartz + H2O + plagioclase + biotite + garnet + melt + ilmenite + sillimanite, which is stable from 675°C to 750°C and 3.5–7.5 kbar. The final assemblage of quartz + H2O + plagioclase + biotite + K-feldspar + garnet + ilmenite + sillimanite is stable from 625°C to 700°C and ~3–5 kbar. Garnet near-rim geothermobarometry for 21-AC-06 predicts a P–T range from ~700°C to 800°C and ~3–7 kbar. These results are broadly consistent with the rim isopleths but hotter by up to 100°C.

7.1.3. 21-AC-04

In sample 21-AC-04, the presence of kyanite and subsolidus garnet constrains the peak subsolidus assemblage to quartz + H2O + plagioclase + muscovite + biotite + garnet + rutile + kyanite, which is stable from 625°C to 700°C and 7.5–12 kbar. The peak anatectic assemblage of quartz + H2O + plagioclase + melt + muscovite + biotite + garnet + rutile + kyanite is stable from 675°C to 725°C and 5–8 kbar. The kyanite-out and sillimanite-in reaction constrains the final metamorphic assemblage of quartz + H2O + plagioclase + muscovite + biotite + garnet + ilmenite + sillimanite to ~600°C–700°C and ~3.5–8 kbar.

The results of geothermobarometry range from ~680°C to 740°C and ~9–13 kbar. These temperatures are consistent with the ~700°C peak temperature predicted by garnet rim isopleths. GASP barometry predicts pressures of several kbar in excess of rim isopleths. This discrepancy is explored in more detail below.

8.1. Anatectic Reaction Textures

All three samples analyzed record evidence that they experienced some degree of anatexis. Although the Nashoba Formation has been recognized as a high-grade anatectic unit for some time [46], the mapping of migmatite locations has focused on K-feldspar bearing stromatic migmatites. While this approach may be a necessary concession to mapping such a large geographic area, it excludes migmatization dominated by biotite-melting or early K-rich melt loss and restitic compositions, both of which may result in rocks that do not have textbook migmatite textures in the field. Below, we summarize the evidence for biotite-out anatexis and potential melt loss in each sample.

8.1.1. JWB13-8

JWB13-8 is muscovite + K-feldspar absent and preserves leucosomes of plagioclase + quartz, which are interpreted to be recrystallized melt from the reactions described below. The presence of the sillimanite intergrown with matrix biotite and the aluminous nature of the bulk composition suggests muscovite was present during prograde metamorphism. Additionally, the results of geothermobarometry, the intersecting isopleth method, and phase equilibria modeling place peak temperature conditions above the solidus. This sample also contains rare, small, chemically unzoned peritectic garnet (G2) enclosed in leucosome. Core and rim isopleths from G2 garnet retrieve the same conditions as near rim isopleths from subsolidus garnet.

The two peritectic reactions are:

  1. Biotite + Plagioclase + Kyanite = Garnet + Melt

  2. Biotite + Plagioclase + Sillimanite = Garnet + Melt

Reaction (1) involves the dehydration melting of biotite to produce melt and a peritectic garnet phase [49, 67-70]. All reasonable P–T trajectories from the peak subsolidus assemblage require some amount of melting at kyanite-grade. This reaction is superseded during exhumation into the peak anatectic field by reaction (2) or biotite melting in the sillimanite stability field. These reactions are both responsible for the generation of rare peritectic garnet, and back-reaction with melt is potentially responsible for the biotite + plagioclase ± sillimanite rims around some subsolidus garnet. The lack of retrograde muscovite in this specific thin section, though late muscovite is present at the outcrop scale, suggests some loss of the earliest, highest-K melts, or sufficiently rapid exhumation or cooling to impede substantial muscovite recrystallization.

8.1.2. 21-AC-06

Sample 21-AC-06 belongs to the same map unit as JWB13-8 (Snsg) and records essentially the same textures with minor differences. In contrast to JWB13-8, 21-AC-06 records better developed gneissic banding on the thin section scale, though in outcrop it is best described as patch migmatite [63]. While the thin section explored here does not contain peritectic garnet, peritectic grains are clearly visible at the outcrop scale (online Supplementary Figure S2). Garnet in this thin section, however, does contain polyphase biotite + quartz + plagioclase inclusions that are likely melt-inclusions trapped during the growth of garnet as a result of reaction (2). Additionally, phase equilibria constraints and geothermobarometry indicate peak temperature conditions above the solidus. The lack of muscovite in this sample, coupled with the patchy outcrop texture suggests that there has been some segregation and potential migration of early melt.

8.1.3. 21-AC-04

The anatectic textures preserved in 21-AC-04 are broadly similar to those in the previous samples. The sample is a patch to stromatic migmatite gneiss at the outcrop scale (online Supplementary Figure S1; [63]) and preserves lenses of dominantly biotite + sillimanite ± relict kyanite restite. Elsewhere, the sample is muscovite poor (<5 mode%), with a matrix foliation defined by intergrown biotite + sillimanite melanosome. Phase equilibria modeling predicts initial melting via reaction (1) at ~700°C and ~8 kbar which is then followed by reaction (2) at <8 kbar. The presence of rare (<5 mode%) muscovite in recrystallized leucosomes suggests 21-AC-04 has experienced the least melt segregation or migration of the samples explored.

All three samples explored record textural evidence for anatexis at the thin section scale. Outcrop scale evidence is more difficult to discern due to both the limited exposure of these outcrops and the potential restitic nature of these units. In general, if quartz + plagioclase layers are considered leucocratic, the Snsa and Sngba units can be described as a patch to stromatic quartz-plagioclase migmatites, interbedded with layers of biotite + sillimanite restite, though we argue it is both simpler and more prescriptive to describe these units as “variably restitic.” The more subdued field appearance of these rocks relative to classic Nashoba Formation K-feldspar migmatites is likely due to the leucocratic phases present (plagioclase + quartz) and variable melt loss and probably the main reason for their mapped designation. Our results suggest that other samples of the same or similar units should be investigated to determine the geographic extent of anatexis preserved by the Nashoba terrane.

8.2. Reconciling Thermodynamic Methods

Assuming metamorphism at near-equilibrium conditions and the absence of extensive retrogression, the results of rim geothermobarometry, thermodynamic modeling, and the intersecting rim isopleth method should be consistent. The error parallelogram defined by rim + core pairs in JWB13-5 and 21-AC-06 overlap with the core isopleths and the peak field (Figure 7). The weak or flat major element zoning in garnet from both samples suggests that garnet and biotite compositions were modified via diffusion at high T and that any retrograde reactions acting in these samples must have occurred at conditions near enough to peak temperature as to result in indistinguishable results. This is consistent with the consumption of garnet during initial (nearly) isothermal decompression as predicted from thermodynamic modeling (see online Supplementary Figures S4–S6).

Core + rim garnet-biotite thermometry and GASP barometry for 21-AC-04 present a more complicated picture. Garnet near-rim isopleths plot within the melt + sillimanite field but are likely within error of the melt + kyanite field. Additionally, the results of GASP barometry are slightly higher pressure than the core isopleths, while garnet-biotite thermometry is well within the error of the core isopleths. At temperatures of ~700°C, it is likely that garnet-biotite maintained partitioning equilibrium. The discrepancy between GASP barometry and rim isopleths is likely spurious, as the composition of plagioclase in this sample (XAn = 0.05) falls outside of the composition range of the GASP calibration. As a result, P–T history interpretations for 21-AC-04 will not utilize the results of GASP barometry. Additionally, whereas phase equilibria modeling predicts some early melting within the kyanite stability field, any heating in that field would produce garnet rather than consume it (online Supplementary Figure S6). In order to produce the resorbed texture preserved in garnet, nearly isothermal decompression into the sillimanite + melt field is required.

While we are able to quantitatively determine the P–T conditions of anatexis, petrographic inference must be used to constrain pre-anatectic conditions. In all three samples, partial or complete diffusional resetting of garnet core and growth compositions has occurred. As a result, the determination of garnet nucleation and growth conditions from the method of intersecting isopleths [47, 71, 72] is impossible. Since we cannot quantitatively determine subsolidus P–T histories via methods that employ mineral chemistry, interpretations of subsolidus P–T histories are limited to the determination of “peak subsolidus assemblages” right before anatexis via MAD modeling and petrographic analysis. This peak subsolidus assemblage is similar to the “paleo assemblage” presented in Castro and Spear [73] but represents the mineral assemblage right before melting, rather than right before garnet nucleation. It is also important to note that MADs were calculated with the minimal water content needed to maintain fluid saturation. While the presence of prograde muscovite and biotite suggests this a reasonable assumption for subsolidus conditions, the calculated melt-in isograd is the fluid-present reaction. This represents the minimum possible melting temperature. Thus, the anatectic conditions referred to in this study represent minimum estimates. Nevertheless, the presence of retrograde muscovite and recrystallized leucosomes at the hand sample-scale suggests the chosen bulk compositions are sufficient.

8.3. P–T History and Comparison to Previous Studies

The peak anatectic conditions presented in this study are similar to that of Abu-Moustafa and Skehan [46] but 2–4 kbar deeper than what is presented in Stroud et al. [7]. Additionally, we present a history of early intermediate P/T metamorphism not recorded in previous studies. All three samples experienced clockwise P–T histories. JWB13-8 records exhumation and associated melting from peak subsolidus conditions of ~650°C–700°C and 7–10 kbar to ~715°C and 5.7 kbar. 21-AC-06 records coupled exhumation and anatexis from ~550°C to 650 °C and 6–9 Kbar to ~700°C and 5 kbar. 21-AC-04 records a history of exhumation and associated anatexis from ~650°C to 700°C and 8–12 kbar to ~700°C and 8.5 kbar. These paths imply a conservative minimum of ~1.3 kbar (~3 km) of exhumation between subsolidus and peak anatectic conditions. While the record of the earliest portion of these P–T paths is lost, we argue that regional metamorphism up to intermediate P–T conditions is the simplest scenario to produce the peak subsolidus assemblage of each sample. This contrasts with the period of early static-thermal metamorphism previously presented [5, 7]. Nashoba Formation sediments are thought to have been deposited within the Tetagouche-Exploits arc/back-arc basin west of the trailing margin of Ganderia before deformation during the Acadian orogeny [3]. The early, intermediate P/T metamorphism described here cannot be driven solely by contact metamorphism due to arc magmatism, as this would require an exceptionally thick sedimentary sequence of ~34 km for our deepest sample, assuming a density of 2.7 g×cm3, and emplacement of arc magmas at the peak P–T conditions determined here. Instead, clockwise intermediate P–T paths are typical of collisional tectonics [19, 20]. The presence of partially overprinted lower grade (e.g., andalusite grade) assemblages elsewhere in the Nashoba terrane [7, 74, 75] may reflect the exhumation of shallower, colder portions of the orogen that did not experience the same intermediate P–T paths explored here. Taken as a whole, our results suggest a three-stage clockwise P–T path characterized by: (1) thickening and burial up to kyanite grade conditions and incipient anatexis, (2) initially isothermal exhumation to sillimanite grade conditions, and (3) continued exhumation to present exposure (Figure 8).

8.4. Tectonic Implications

The results from this study allow us to bring finer detail to our understanding of the tectono-metamorphic history of the Nashoba Formation. Previous studies have suggested a low-P/T (less than ~5 kbar) history for Nashoba metamorphism due to the westward subduction under the trailing edge of Ganderia and associated arc magmatism, followed by anatexis due to the collision of Avalonia during the Acadian orogeny [5, 7]. The samples studied here do not preserve evidence for early, low P/T contact metamorphism due to arc magmatism but instead present compelling evidence for intermediate P/T metamorphism prior to anatexis with peak pressures ~2–4 kbar higher than previously published estimates (Figure 8; [7]). There are two scenarios that can explain our results: (1) early intermediate to high-P/T metamorphism is the result of regional scale crustal thickening and subsequent exhumation due to back-arc closure during the ~430 Ma Salinic orogeny (Figure 9(a)) or (2) the entire P–T history preserved in these rocks is the result of thickening and exhumation during the later Acadian and Neoacadian orogenies (Figure 9(b)).

Existing geochronological constraints make scenario (1) unlikely, as only Hepburn et al. [3] report metamorphic zircon U-Pb LA-ICPMS ages prior to the Acadian orogeny (older than 425 Ma), but weighted averages of metamorphic zircon ages are no older than ~433 Ma, so this metamorphism could be related to the onset of the Acadian orogeny. Intermediate to high P/T conditions due to crustal thickening have been observed in back-arc basin closures. For example, peak P–T conditions of the Imjingang belt, a relict back-arc basin within the Korean Peninsula, were reported to reach ~710°C and 11 kbar [21, 22, 24]. Metamorphic zircon ages in the Imjingang belt were found to match the timing of back-arc basin closure. Similarly [23, 76, 76] document clockwise P–T paths and peak P–T conditions of ~625°C and 6-10 kbar associated with back-arc closure in the Cordillera Darwin, Chile. Although those peak conditions and P–T paths roughly match those reported in this study, the bulk of reported metamorphic zircon and monazite ages postdate back-arc basin closure in Ganderia [3, 7].

The second scenario requires that the P–T histories presented (e.g., Figure 8) here are entirely the result of the collision and docking of Avalon during the Acadian and Neoacadian orogenies. In this scenario, the Nashoba formation experienced intermediate/high-P metamorphism, up to kyanite-grade as a result of crustal thickening and burial driven by the onset of the Ganderia-Avalonia collision during the Acadian orogeny. Continued collision and heating drove anatexis, which likely facilitated nearly isothermal exhumation of the Nashoba terrane through sillimanite-grade conditions. Taken as a whole, the entire P–T history recorded in the Nashoba Formation may be consistent with Devonian collision of Avalon onto the Ganderian trailing edge, (i.e., the composite Laurentian margin).

Our results also suggest upper amphibolite facies conditions and anatexis were more pervasive in the Nashoba Terrane than previously predicted. Such conditions are consistent with a proposed channel flow model of exhumation, in which more ductile, viscous rock is able to flow between two blocks of more rigid rock [77]. The amphibolite-facies Nashoba terrane is juxtaposed between two greenschist-facies blocks, the Merrimack belt to the west and the Avalonia terrane to the east (Figure 1), which could represent the bounding blocks on either side of the flowing channel [14]. Additionally, the shear zones between these terranes have yielded opposing shear sense indicators (normal sense on the Clinton-Newbury Fault Zone, reverse sense on the Bloody Buff Fault Zone) that are consistent with the predictions of the channel flow model [14]. However, detailed petrochronology across the Nashoba terrane is necessary to definitively determine if the terrane represents a tectonically coherent unit as predicted by the channel flow model or a collage of crustal slices that were juxtaposed after peak metamorphism. Regardless, the Nashoba terrane likely exposes rocks from greater structural depth within the orogenic hinterland than previously believed.

The application of petrographic analysis, thermodynamic modeling, and geothermobarometry has allowed us to refine our understanding of the metamorphism recorded in the Nashoba Formation. The samples in this study are mapped as nonmigmatite schists and gneisses, but we present evidence that they experienced anatexis that was predominantly driven by biotite-out melting. We also present novel evidence for intermediate-P/T metamorphism in the Nashoba Formation. Our samples experienced a clockwise P–T path, beginning with early, andalusite-grade metamorphism from burial and arc-related heat influx, followed by up to kyanite-grade metamorphism from crustal thickening, and finally coupled anatexis and exhumation from ~550°C to 700°C and ~6–12 kbar to ~700°C–715°C and 5–8.5 kbar. We outline two possible scenarios for this metamorphic history: early intermediate to high-P metamorphism from relict regional scale crustal thickening predating the Acadian orogeny, or a P–T history as a result of thickening and exhumation entirely within the Acadian orogeny. Our results are also generally consistent with previously purposed models of the ductile extrusion for the Nashoba Terrane. Further detailed petrochronology is necessary to definitively discern which model is correct and to determine if the metamorphic fingerprints of channel flow are present. Nevertheless, our results suggest that the Nashoba terrane exposes a hotter and deeper portion of Ganderia than previously believed.

The data that support the findings of this study are available either within the article itself (e.g., mineral chemistry and sample location tables) or from the corresponding author upon reasonable request (e.g., SPaC dataset).

The authors declare that there is no conflict of interest regarding the publication of this paper.

The authors wish to thank Freya George, Frank Spear, and Jay Thomas for thoughtful discussion of all things aluminosilicate and anatexis. We thank Wes Buchanan for graciously loaning sample JWB13-8 to A. Castro. We wish to thank Jay Thomas, Celine Martin, and Nicholas Tailby for expert help in collecting electron microprobe analyses. We thank Dan Brabander and Claire Hayhow for assistance with XRF analyses. We also thank J. Christopher Hepburn and Yvette Kuiper for thorough and helpful reviews on early versions of this manuscript. We thank two anonymous reviewers for insightful comments that helped improve the presentation of this paper. We thank Bo Wang for helpful editorial handling. This work was funded by the Hawk endowed fund in Geosciences at Wellesley College and the Wellesley College Provost’s Office.

Supplementary data