Garnet sand reveals rock recycling processes in the youngest exhumed high- and ultrahigh-pressure terrane on Earth

The placer deposit contains (sub-)angular detrital garnet grains derived from a proximal source. Electron microprobe analyses of garnet from outcrop samples (compiled from the literature and this study) and detrital garnet (this study) allow a comparison of garnet compositions to assess potential sediment sources (Fig. 3 and SI Appendix). Ternary plots of garnet compositions (Fe + Mn, Ca, and Mg end-members) for crystalline rocks (n = 881) compared with detrital garnet compositions (n = 716) indicate that the detrital garnet grains cover the entire compositional range of potential crystalline source rocks.

Sixty-two percent of the detrital garnet grains investigated contained mineral inclusions. These include 83 garnets from the 63- to 125-µm-sized fraction, 117 garnets from the 125- to 200-µm-sized fraction, and 154 garnets from the >200-µm-sized fraction, with the largest grain having a long axis of ∼600 µm (SI Appendix). Detrital garnet mineral inclusions (≥2 µm) identified using Raman spectrometry include omphacite, kyanite, epidote, oxides (SiO₂ polymorphs, rutile), feldspars (alkali feldspar, plagioclase feldspar, kokchetavite), phyllosilicates (paragonite–muscovite series, phlogopite–biotite, chamosite–clinochlore, and pyrophyllite), and accessory phases (apatite, zircon, titanite) (Fig. 4 and SI Appendix). Rutile is the most abundant mineral inclusion, occurring in 73% of the inclusion-bearing detrital garnet grains studied. At least 2% of inclusion-bearing detrital garnets are inferred to be sourced from mafic eclogite as indicated by omphacite inclusions in high-Mg garnets (Fig. 3).

Three SiO₂ polymorphs (coesite, cristobalite, and quartz) were identified as inclusions in detrital garnets (Fig. 4). Quartz is the most abundant SiO₂ inclusion, occurring in 30% of the inclusion-bearing detrital garnets. Two garnets in the >200-µm fraction were sourced by UHP protoliths as indicated by a 4.0 × 3.0-µm polyphase SiO₂ inclusion, consisting of quartz and coesite (Fig. 4A), and a 6.5 × 5.5-µm monomineralic coesite inclusion (Fig. 4B). The garnet containing a polyphase coesite/quartz inclusion also contains carbonate, epidote, and CO₂ inclusions (Fig. 4A), whereas the monomineralic coesite inclusion coexists with numerous graphite inclusions in the same garnet grain (Fig. 4B).

Eleven detrital garnets contain cristobalite, the metastable highest-temperature polymorph of SiO₂. Raman spectra (Fig. 4C) indicate the inclusions are likely α-cristobalite (20). Nine (of 13) cristobalite crystals occur as polyphase inclusions, often together with mica (muscovite–paragonite, pyrophyllite, or chamosite–clinochlore) ± carbonate ± kokchetavite (Fig. 4C and SI Appendix). Kokchetavite, a metastable polymorph of K-feldspar (21), occurs in 4 of the 11 cristobalite-bearing garnets either in the same polyphase inclusion or as an adjacent inclusion. Five additional garnets were found to contain kokchetavite as polyphase inclusions (± quartz ± carbonate ± muscovite–paragonite ± rutile) (SI Appendix). The compositions of cristobalite- and kokchetavite-bearing garnet grains span the entire compositional range of the measured detrital garnets (Fig. 3).

Carbon-bearing inclusions consist of carbonate minerals (26% of inclusion-bearing garnets), graphite (22%), and CO₂ (19%) (Fig. 4 and SI Appendix). The compositions of carbon-bearing garnet grains also span the entire compositional range of the measured detrital garnets. However, graphite inclusions do not occur in high-Mg garnets. Therefore, we infer that carbon-bearing detrital garnets were derived from felsic sources (Fig. 3).

Six detrital garnets contained monomineralic elastically isolated quartz and zircon inclusions (i.e., away from grain boundaries, cracks, or other inclusions). The remnant strains of elastically isolated monomineralic quartz and zircon inclusions were used to constrain the P-T conditions of elastic equilibration within garnet hosts (Fig. 5). Application of elastic geothermobarometry requires that the inclusions remained isolated and fully contained within the host mineral after entrapment in order to preserve elastic strains from inclusion entrapment (22). However, inclusion strains can be modified through elastic, brittle, or ductile deformation (23–25), leading to modification of the calculated P-T conditions of elastic equilibration. Small (<10-μm-diameter) monomineralic inclusions contained within six detrital garnets from the >200-μm-size fraction were analyzed using Raman spectroscopy, ensuring that there is sufficient garnet volume to elastically isolate the inclusions (25–27) (SI Appendix). Each inclusion was verified to be more than three times the inclusion radii beneath the surface of the garnet and away from any visible inhomogeneities in the garnet host (i.e., other inclusions, fractures, or grain boundaries). The 1,008-cm⁻¹ peak width, in terms of full width at half-maximum intensity, of less than 5 cm⁻¹ was used as a criteria in the selection of nonmetamict zircon inclusions for Raman thermobarometry (27) (SI Appendix). Remnant strains were determined for quartz and zircon inclusions from measured Raman shifts using their respective Grüneisen tensor (28). No evidence of microfractures around inclusions was observed from visual inspection of each measured inclusion. Therefore, the remnant strains of the isolated inclusions within the garnet hosts are interpreted as the P-T conditions of elastic equilibration. Strains were used in conjunction with the elastic stiffness tensors of quartz and zircon to calculate the average remnant stress within the inclusions for use in an isotropic elastic model to constrain the range of possible elastic equilibration conditions [i.e., isomekes (29)]. P-T conditions of elastic equilibration between the host and inclusions were determined from the intersection points of isomekes calculated for coexisting quartz and zircon inclusions in six detrital garnets that met the criteria as stated above. Results indicate P-T conditions of elastic equilibration ranging from 0.5 to ∼1.3 GPa and from 550 °C to 750 °C. Despite the large uncertainties, the data suggest elastic equilibration under primarily granulite facies conditions and subordinate amphibolite facies conditions (Fig. 5).

In the eastern Papuan (U)HP terrane outcrop samples of felsic gneisses that preserve UHP indicator minerals have yet to be found. Only one sample of mafic eclogite enclosed within retrogressed felsic gneiss has been found to preserve coesite as a partially transformed inclusion in omphacite (16). The coesite crystallized at 763 °C (+72 °C to −61 °C) and 3.3 GPa (+0.77 to −0.60 GPa) (30) during UHP metamorphism at ∼8 Ma (4, 17, 31) (Fig. 5). Garnet growth in coesite eclogite occurred in two stages (32): 1) core garnet growth at ∼650 °C and >2.7 GPa within the coesite stability field and 2) garnet rim growth following partial exhumation from UHP conditions at <7.1 Ma during a transient (0.3-My) heating event at 1.5 GPa (Fig. 5).

Although it is not possible to extract the detailed P-T paths followed by individual detrital garnet grains, their compositions and inclusion suites can be compared with those from outcrop samples to assess the P-T paths followed by protoliths with a range of bulk compositions, during transit through the rock cycle (Figs. 3 and 5 and SI Appendix). The discovery of coesite inclusions in detrital garnets provides direct mineral evidence that UHP metamorphism also affected felsic lithologies in the eastern Papuan (U)HP terrane (Figs. 2–4). Compared with garnet from the coesite eclogite outcrop, the two coesite-bearing detrital garnets have lower Mg contents and contain carbon-bearing inclusions (i.e., graphite, carbonate minerals, CO₂). These detrital garnets are inferred to have formed from metamorphism of organic-rich sedimentary rocks (Fig. 3 and SI Appendix), providing petrologic evidence for subduction of protoliths originating at the Earth’s surface. No microdiamonds were found as inclusions in detrital garnet. Therefore, protolith conditions are inferred to have reached the coesite stability field but not the stability field for diamond (Figs. 2 and 5). Our findings expand the range of bulk compositions, notably low-density bulk compositions (felsic gneiss and metasedimentary rocks), that preserve evidence for UHP conditions in eastern Papua New Guinea.

Many mineral inclusions in the placer deposit garnets are commonly found as matrix minerals in basement rocks of the D’Entrecasteaux Islands (e.g., eclogite, felsic host gneiss, and amphibolite). However, cristobalite and kokchetavite are discoveries that have not yet been documented in basement rocks. Preservation of hydrous phases in the gneisses and shear zones on Goodenough Island (i.e., peak phengite, retrograde muscovite, biotite and chlorite, serpentinite) and mica and fluid inclusions preserved in detrital garnet span the entire compositional range (Fig. 3 and SI Appendix). We infer that fluids affected all lithologies during transit through the rock cycle. In order to be preserved, metastable mineral inclusions (i.e., coesite, cristobalite, kokchetavite) must have remained isolated systems, trapped within mechanically strong garnet hosts during rapid exhumation. These metastable inclusions also remained insulated from fluids during erosion and transport of sediment from limited catchment areas and deposition to form the placer deposit.

Detrital garnet mineral inclusions sourced by subducted carbonates and organic materials include graphite, carbonate minerals, and CO₂. We infer that subducted carbonaceous material converted to graphite via rapid recrystallization and devolatilization (33). However, subducted carbonaceous material apparently did not reach sufficient depths to form diamond. The abundance of graphite and CO₂ inclusions (22 and 19% of inclusion-bearing garnets, respectively) provides evidence for deep carbon recycling in the subduction channel (34). At present, shallow seeps offshore Normanby and Dobu Islands are releasing CO₂ (>98%) and CH₄ (87 to 4,360 ppm) (35) (Fig. 1B), providing evidence for the return of recycled carbon to the oceans and atmosphere.

Evidence for high-temperature conditions occurs in detrital garnets of all compositions. Polyphase inclusions of cristobalite + kokchetavite + mica (Fig. 4C) are interpreted as trapped partial melt crystallized during rapid decompression (i.e., “nanogranitoids”) (36, 37) (SI Appendix). Elastic thermobarometry of individual detrital garnet grains also suggests elastic equilibration of quartz and zircon inclusions under granulite facies conditions and for one grain in the amphibolite facies (Fig. 5). Such high thermobaric ratios (T/P) are not compatible with equilibration along a prograde path associated with low-T/P subduction gradients but instead, provide additional evidence for heating during exhumation.

Although we are unable to ascertain the potential effects of chemical variations in quartz and zircon inclusions in the detrital garnets used for elastic thermobarometry, we found no evidence for elastic equilibration of quartz and zircon inclusions in detrital garnet corresponding to eclogite facies conditions. Instead, consistent P-T conditions in the high-temperature, lower-pressure metamorphic facies are obtained and cannot be attributed to processes that would variably affect individual inclusions, such as changes in the degree of elastic isolation or brittle deformation of inclusions. Therefore, we infer that the elastically isolated quartz and zircon inclusions in detrital garnets reequilibrated by viscous creep in their garnet hosts at temperatures >600 °C (25). Amphibolite facies conditions determined for one detrital garnet (DG-35) provide additional evidence for retrograde garnet growth, in agreement with quartz-in-garnet and Ti-in-quartz thermobarometry for a quartzofeldspathic gneiss in the catchment area (38) (Fig. 5).

The formation and preservation of both coesite and cristobalite as inclusions in garnet require extreme metamorphic conditions indicative of transient geothermal gradients (i.e., from low T/P to higher T/P). The rock record in eastern Papua New Guinea preserves evidence for Late Miocene low-T/P metamorphism reaching UHP conditions and higher-T/P metamorphism during exhumation, as well as highest-T/P gradients associated with decompression melting of asthenospheric mantle to form basalts in the Woodlark Basin (Figs. 1 and 2). The UHP locality occurs along strike of a sinistral strike–slip transfer fault system from the placer deposit (39); the faults displace Pliocene–Quaternary basaltic lavas (Fig. 1B). Such petrologic variation preserved in the geologic record over short spatial and temporal scales indicates that steady-state processes are not operative within this part of the obliquely convergent AUS–PAC plate boundary zone (11).

The detrital garnet data support models for the transport and transformation pathways of heterogeneous lithologies, including organic and inorganic carbon, in the forearcs of subduction zones (34). Results reinforce the observation that subducted rocks do not have to remain dry during transport through the rock cycle in order to preserve evidence for UHP metamorphism, as long as mineral inclusions (such as coesite) in strong host minerals, like garnet, remain isolated from metamorphic fluids (5). In addition, detrital garnet data provide observational evidence for H₂O and CO₂ storage in the forearc and volatile and melt migration within the forearc during rifting. We infer that the detrital garnet SiO₂ inclusions (coesite, quartz, and cristobalite) capture processes including subduction to upper mantle depths along low-T/P gradients and exhumation associated with an increase in geothermal gradients (i.e., to higher T/P) as the Woodlark Rift formed ahead of the westward-propagating Woodlark Basin seafloor-spreading system (Fig. 1A).

The Papuan (U)HP terrane preserves evidence for geochemical cycling between the surface and the upper mantle, including UHP rock formation and exhumation, that can be directly linked to the rapid spatial and temporal evolution of the obliquely convergent PAC–AUS plate boundary zone. Microplate rotation of the Woodlark plate relative to the AUS plate (40), upper plate divergence (41), and buoyancy forces (42) contributed to rapid (U)HP exhumation during the transition from subduction to rifting in eastern Papua New Guinea (Fig. 1). Following the onset of upper plate divergent motion, the oceanic upper plate and subducted continental lower plate will decouple, providing accommodation space (41) for subducted metasediments, felsic gneisses, and their partial melts to rapidly exhume within the subduction channel. Numerical models that use constraints provided by the geologic record in eastern Papua New Guinea (43, 44) predict that asthenospheric mantle will flow upward toward the subduction channel to fill the accommodation space left by the buoyantly exhuming felsic (U)HP rocks, resulting in isotherm advection, an increase in geothermal gradients, and heating of the (U)HP terrane during exhumation.

The investigation of erosional products at the catchment scale expands the potential for finding additional evidence of UHP metamorphism in the eastern Papuan (U)HP terrane where pervasively overprinted felsic rocks and their partial melts have yet to yield additional evidence in outcrop since coesite was discovered in a mafic eclogite. Results demonstrate that inclusions in detrital garnets from modern sediments preserve evidence for UHP metamorphism. Such geologically short timescales (i.e., since the Late Miocene) for rock recycling within this active plate boundary zone imply <10-My petrologic lag times for metamorphic rocks exhumed from approximately ≤120-km depths, eroded, and deposited as sediment. High-spatial resolution spectroscopic imaging of heavy mineral fractions of sediment derived from the erosion of the Papuan (U)HP terrane provides a means to efficiently search for evidence of UHP metamorphism and to further constrain the thermal evolution and geodynamics of continental subduction. A better understanding of subduction and exhumation processes in active plate boundaries is needed before modern processes can be accurately extrapolated and applied to the geologic record. Raman spectroscopy of mineral inclusions in detrital garnet in both modern sediments and stratigraphic sections, especially when integrated with single-grain Sm–Nd detrital garnet ages (45), has the potential to provide insight into the tempo of rock recycling processes throughout Earth’s history, including on early Earth when plate velocities and transport through the rock cycle were at least as fast as present day.

The beach placer sample (PNG08002c) was dry sieved into sized fractions (63 to 125, 125 to 200, and >200 µm), split by coning and quartering, mounted in epoxy, ground with SiC, and polished with Al₂O₃ to expose grain surfaces. The compositions of 572 detrital garnet grains were obtained by electron microprobe analysis at the Geosciences Center of the University of Göttingen using a JEOL JXA 8900 RL equipped with five wavelength dispersive spectrometers. Measurements were performed with an accelerating voltage 15 kV and a beam current of 20 nA. Counting times were 15 s for Si, Mg, Ca, Fe, and Al and 30 s for Ti, Cr, and Mn. Measurement spots were preferentially set to the garnet centers but shifted toward the rim when inclusions or fractures are located in the center. For the coesite-, cristobalite-, and kokchetavite-bearing garnet grains, nine spots per garnet were set: one at the center, four at the mantle, and four at the rim.

From the investigated 572 garnet grains, 354 contain inclusions ≥2 µm, which were identified using a Horiba Jobin Yvon XploRA Plus spectrometer equipped with an Olympus BX41 microscope at the Geosciences Center of the University of Göttingen (10). Analytical conditions include a 532-nm excitation laser, a 1,800-L mm⁻¹ grating, a 100× long working distance objective with a numerical aperture of 0.8, a confocal hole diameter of 100 µm, and a slit of 100 µm. Two-dimensional Raman images of monomineralic coesite, bimineralic coesite/quartz, and polyphase inclusions containing cristobalite and kokchetavite were collected with a WITec alpha300R ultrahigh-throughput Raman spectrometer at the Geosciences Center of the University of Göttingen. The spectral images were acquired with a 532-nm excitation laser, an automatically controlled laser power of 20 mW, a 1,200-L mm⁻¹ grating, and a 100× long working distance objective with a numerical aperture of 0.75. Spectra were collected at a step size of between 50 and 200 nm using an acquisition time of between 1 and 2 s per spectrum, depending on the specific inclusion. Automated cosmic ray correction, background subtraction, spectral averaging/smoothing, and supervised component analysis were performed using the WITec Project software.

A Renishaw inVia Raman spectrometer with a confocal optical microscope at Syracuse University was used to measure the Raman shifts of quartz and zircon inclusions in detrital garnet for elastic thermobarometric calculations following procedures for quartz-in-garnet inclusions (38, 46) and zircon-in-garnet inclusions (27). All measured inclusions were elastically isolated and located more than three inclusion radii from other inclusions, grain boundaries, and visible defects in the host garnet. Because metamictization of zircon modifies the zircon bulk elastic properties (27), only nonmetamict zircon inclusions were selected for analysis as determined using the 1,008-cm⁻¹ peak width, in terms of full width at half-maximum (FWHM) intensity, of less than 5 cm⁻¹ (SI Appendix). A 532-nm laser was focused onto specimens with a 100× microscope objective with a numerical aperture of 0.9. The Raman scattered light was statically dispersed with a 1,800-L mm⁻¹ grating onto a charged couple device, resulting in a spectral resolution of 0.5 cm⁻¹. The spectrometer was calibrated against neon lines and a silicon standard. Spectral accuracy and linearity were checked throughout each analytical session by measuring the Rayleigh scattered light from the 532-nm laser, the 520.5-cm⁻¹ Raman band of a silicon standard, the Raman bands of a synthetic quartz reference material from the Westinghouse Corporation, and the Raman bands of a synthetic zircon reference material. All Raman spectra were acquired for 20 s and measured at room conditions of 23 °C and 1 bar. Spectra were not processed or corrected prior to peak fitting using a pseudo-Voigt function in the Renishaw software. Errors on fitted band positions are ∼0.2 to 0.3 cm⁻¹.

For strain calculations, Raman band frequencies of quartz were measured at ∼127.5, 205.9, and 464.8 cm⁻¹, and the Raman band frequencies of zircon were measured at 213.4, 224.0, 356.0, 438.8, 974.8, and 1,008.7 cm⁻¹. Synthetic quartz and zircon reference materials were measured between sample measurements to derive calibration curves in the case of any spectral drift. Measurements were performed near the centers of the inclusions, where the effects of inclusion geometry and elastic anisotropy have the least effect (26). Furthermore, it has recently been demonstrated that when elastically anisotropic inclusions, such as quartz and zircon, are contained within nearly isotropic host minerals such as garnet, the effect of elastic anisotropy on calculated entrapment conditions is <2% and therefore, negligible as compared with the errors on Raman measurements (47). The remnant strains (ε₁ + ε₂, ε₃) in inclusions were calculated using the fitted Raman band positions using the stRAinMAN software and the respective Grüneisen tensors of quartz and zircon (28). Errors on fitted band positions of 0.3 cm⁻¹ were used in stRAinMAN to calculate the estimated standard deviation (esd) and covariance matrix needed for error calculation of the remnant inclusion pressure (48). Remnant strains of quartz and zircon inclusions were used in conjunction with the respective elastic tensor to calculate the pressure within the inclusion. Elastic modeling of inclusion pressures and entrapment isomekes was performed for quartz and zircon inclusions in an almandine host using an isotropic elastic model implemented in the EntraPT program (49). Errors on remnant strains were propagated to entrapment isomekes using the covariance matrix and esd on inclusion strains according to the method described in ref. 49. Errors on calculated entrapment isomekes, especially for zircon inclusions in garnet, varied over a large temperature range. Inclusions with the largest errors were those that had the least number of fitted Raman bands and were subsequently excluded from P-T calculations and interpretation.

S.L.B. acknowledges support from NSF Division of Earth Sciences Grant 0709054, NSF Major Research Instrumentation Program Grant 1625835, and the Thonis Endowment at Syracuse University. J.S. and H.v.E. acknowledge support from German Research Foundation Grant EY 23/27-1. J.P.G. acknowledges support from a Syracuse University Science Technology Engineering Math fellowship and a European Research Council Postdoctoral Fellowship (Grant 714936 to M. Alvaro). We thank R. M. Davies for garnet electron microprobe analyses and P. G. Fitzgerald for providing comments on an earlier version of the manuscript. We acknowledge the constructive reviews of Editor Peter Kelemen and two anonymous reviewers, which helped us improve our manuscript. Collaborative discussions on detrital analysis with S. Andò, M .G. Malusà, E. Garzanti, and P. G. Fitzgerald are greatly appreciated.