JGR–Solid Earth

Large bursts of non-volcanic tremor (“major” tremor episodes) correlated with geodetic deformation recur regularly in the Cascadia subduction zone and are often called episodic tremor and slip (ETS). Minor episodes of tremor between ETS are ubiquitous but have been understudied. This paper assesses time-invariant characteristics of tremor episodes of all sizes within northern Cascadia. We derive a catalog of tremor episodes ranging in size from 10 to >13,000 tremor events using the results of 17 years of tremor monitoring. Minor episodes represent ∼96% of all 896 tremor episodes and their occurrence varies on 10-km scales. Using estimates for the depth of the forearc Moho and subducting slab, we observe an association between the location of the forearc mantle corner (FMC) and tremor occurrence that leads to along-dip modality. Bimodality, present in southern Washington and Vancouver Island, represents the segmentation of major and minor episodes up-dip and down-dip of the FMC, respectively. Unimodality, present in Puget Sound, results when the FMC is located near the down-dip edge of the ETS zone and no segmentation occurs. We also use our extensive tremor episode catalog alongside three-dimensional regional tomographic velocity models to reassess the relationship between tremor activity and low Vp/Vs signatures in the forearc. We do not find a correlation between tremor episode recurrence intervals and Vp/Vs, contrary to some previous work, suggesting that controls on silica precipitation in the forearc crust are not dominant controls of tremor episode recurrence, or that the association is not widely observable.

This paper addresses one of the critical questions of scientific inquiry: How do we know when a given data set is representative of the phenomenon being examined? For paleomagnetists, the question is often whether a particular data set sufficiently averaged paleosecular variation (PSV). To this aim, we updated an existing PSV data set that now comprises 2,441 site mean directions from 94 individual studies (PSV10-24). Minimal filtering for data quality resulted in 1,619 sites from 90 publications. Fitting PSV10-24 with two newly defined parameters as well as two existing ones form the basis of a Giant Gaussian Process field model (THG24) consistent with the data. Drawing directions from THG24 yields directional distributions predicted for a given latitude allowing a comparison between empirical distributions and the cumulative distribution function generated by the model. This tests whether the observed data adequately averaged out PSV according to THG24. We applied these tests to five data sets from Large Igneous Provinces from the last billion years and find that they are consistent with the THG24 model as well. Sedimentary data sets that may have experienced inclination shallowing can be corrected using an (un)flattening factor that yields directions satisfying THG24 in a newly-defined, four-parameter space. This approach builds on the Elongation-Inclination (E/I) method of Tauxe and Kent (2004), https://doi.org/10.1029/145gm08, so the approach introduced here is called SVEI. We show examples of the use of SVEI and explain how to use this newly developed Python code that is publicly available in the PmagPy GitHub repository.

To investigate continental dynamics underneath the south-central Tibetan plateau, which composes the Himalayan, Lhasa, and Qiangtang blocks, we have conducted comprehensive examinations of seismic azimuthal anisotropy in the crust using receiver functions (RFs) and crustal and mantle anisotropy using teleseismic shear wave splitting (SWS) analysis. In the Qiangtang block, the observed predominantly E-W fast orientations from RF and SWS analyses with similar magnitude are interpreted as resulting from eastward crustal flow with minor contributions from the mantle. In the Lhasa block, the crustal anisotropy is approximately N-S oriented, which is parallel to the strike of rift basins and southward crustal flow. Anisotropy revealed by SWS demonstrates a rotation from E-W in the north to NE-SW in the south, which can be interpreted as reflecting mantle flow field induced by the northward movement of the subducting Indian plate. The addition of PKS and SKKS measurements and extension of epicentral distance range to 171.8° for SWS analysis revealed dominantly strong E-W oriented anisotropy in most parts of the Himalayan block, where most previous studies reported pervasively null measurements. The absence of azimuthal anisotropy is observed in two regions in the Himalayan block which is attributable to mantle upwelling through a previously identified slab window. A two-layered anisotropy structure with different fast orientations for the upper and lower layers can be constrained in the southern Qiangtang and the vicinity of the Main Boundary Thrust.

The 1999 basaltic eruption of Shishaldin volcano (Alaska) displayed a transition between Subplinian and Strombolian activity. Strombolian bubbles indicate the presence of a periodically unstable foam at the top of magma reservoir. In contrast, a long foam, whose rupture led to the eruptive column, was also able to collect in the conduit. Laboratory experiments show that long foams can be produced in a conduit by the spreading of a stable foam accumulated at the top of the reservoir. The existence of a Taylor bubble at the onset of the Subplinian phase, also reproduced by my laboratory experiments, suggests that the foam in the reservoir was just at the transition between stable and unstable. This constrains the bubble diameter prior to the Subplinian phase to be 0.034–0.038 mm when using the foam dimensionless analysis and the underlying gas flux (0.52–0.80 m3/s). The increase in bubble diameter and potentially gas flux prior to the Strombolian activity, 0.81–1.4 m3/s, is sufficient to explain the foam in transition to be unstable. The radius of the magma reservoir is small, 200–210 m, as expected. The bubble diameter is the smallest of those estimated for classical basaltic eruptions (Etna, Kı̄lauea, Erta 'Ale), while the gas flux is among the largest. A dilute suspension of small and isolated bubbles cannot explain the large gas flux at Shishaldin. This implies numerous bubbles with a gas volume fraction ≥0.63−2%, a regime for which the bubbles form bubble clusters. The diameter of these bubble clusters, 3.0–5.4 mm, is sufficient to explain large gas fluxes.

Despite significant progress in studying subduction earthquake cycles, the vertical deformation is still not well understood. Here, we use a generic viscoelastic earthquake-cycle model that has recently been validated by horizontal observations to explore the dynamics of vertical earthquake-cycle deformation. Conditioned on two dimensionless parameters (i.e., the ratio of earthquake recurrence interval T to mantle Maxwell relaxation time t M (T/t M ) and the ratio of downdip seismogenic depth D to elastic upper-plate thickness H c (D/H c )), the modeled viscoelastic deformation exhibits significant spatiotemporal deviations from the simple time-independent elastic solution. Caution thus should be exercised in interpreting fault kinematics with vertical observations if ignoring Earth's viscoelasticity. By systematically exploring these two parameters, we further investigate three metrics that characterize the predicted vertical deformation: the coastal pivot line (CPL), the uplift zone (UZ) landward above the downdip seismogenic extent, and the secondary subsidence zone (SSZ) in the back-arc region. We find that these metrics can all be time-dependent, subject to D/H c and T/t M . The CPL location and the UZ width are mainly controlled by D/H c and T/t M , respectively. The presence of the SSZ is prevalent during the interseismic phase due to viscous mantle flow driven by ongoing plate convergence. Contemporary vertical deformation in Nankai and Cascadia is largely consistent with the model predictions and features differences mainly related to contrasting D/H c values in the two margins. These findings suggest that vertical crustal deformation bears fruitful information about subduction-zone dynamics and is potentially useful for inversions of key subduction-zone parameters, deserving properly designed monitoring.

Pressure and stress perturbations associated with volcanic activity and geothermal production can modify the porosity and permeability of volcanic rock, influencing hydrothermal convection, the distribution of pore fluids and pressures, and the ease of magma outgassing. However, porosity and permeability data for volcanic rock as a function of pressure and stress are rare. We focus here on three porous tuffs from Krafla volcano (Iceland). Triaxial deformation experiments showed that, despite their very similar porosities, the mechanical behavior of the three tuffs differs. Tuffs with a greater abundance of phyllosilicates and zeolites require lower stresses for inelastic behavior. Under hydrostatic conditions, porosity and permeability decrease as a function of increasing effective pressure, with larger decreases measured at pressures above that required for cataclastic pore collapse. During differential loading in the ductile regime, permeability evolution depends on initial microstructure, particularly the initial void space tortuosity. Cataclastic pore collapse can disrupt the low-tortuosity porosity structure of high-permeability tuffs, reducing permeability, but does not particularly influence the already tortuous porosity structure of low-permeability tuffs, for which permeability can even increase. Increases in permeability during compaction, not observed for other porous rocks, are interpreted as a result of a decrease in void space tortuosity as microcracks surrounding collapsed pores connect adjacent pores. Our data underscore the importance of initial microstructure on permeability evolution in volcanic rock. Our data can be used to better understand and model fluid flow at geothermal reservoirs and volcanoes, important to optimize geothermal exploitation and understand and mitigate volcanic hazards.

The induced polarization (IP) method holds a strong potential to better characterize the critical zone of our planet especially in areas characterized by multi-phase flow. Power-law relationships between the bulk, surface, and quadrature conductivities versus the pore water saturation are potentially useable to map the subsurface water content distribution. However, the saturation exponents n and p in these power-law relationships have been observed to vary with the texture of geomaterials and the wettabilities of pore fluids. Traditional experimental setups in the laboratory do not allow to independently visualize the pore fluid distribution. Therefore, the physical interpretations of the two saturation exponents have remained unclear. We developed a novel milli-fluidic micromodel using clay-coated glass beads that exhibit excellent visibility and high IP response. Through laboratory experiments, we simultaneously determined the micromodel complex conductivity and acquired the corresponding pore-scale fluid distributions generated by drainage and imbibition through such class of porous materials. Finite-element simulations of complex conductivity based on the upscaling of the complex surface conductance of grains were conducted to determine the saturation exponents under ideal pore fluid distributions. Results indicate that saturation exponents n and p vary depending on the ganglia size of the insulating fluids. The saturation exponents n and p exhibit power-law relationships with the change rate of pore water connectivity with saturation, which is calculated through the computation of the derivative of Euler characteristics. These findings provide a new physical explanation to the relationships between the saturation exponents and the microscopic fluid distributions within the geomaterials.

Although blind normal faults are common in subduction environments, their rheology, kinematics and interaction with the upper crust are poorly constrained. A month-long shallow normal faulting sequence in the Ibaraki-Fukushima prefectural border (IFPB), northeast Japan, which followed the M w 9.0 Tohoku-Oki earthquake (TOE) and culminated in the M w 6.7 Iwaki earthquake, provides a window into megathrust-to-normal fault interaction. Stress change calculations indicate that direct triggering by the TOE co- and post-seismic slip does not provide a plausible explanation for the IFPB earthquake sequence. In quest for an alternative triggering mechanism, we analyzed post-TOE GNSS data from eastern IFPB. A key step in this analysis is the removal of the large-scale post-TOE displacement field, after which a distinct highly-localized strain along the coastline becomes apparent. The accumulation of this strain was mostly aseismic, and migrated with time prior to the Iwaki earthquake in a manner that correlates well with post-TOE local seismicity. We attribute the pre-Iwaki earthquake strain accumulation to aseismic slip along low-angle seaward dipping blind normal fault, activated by the TOE. Stresses transferred by this slip episode accelerated the failure along the IFPB shallow normal faults. This indirect triggering of the Iwaki earthquake sequence by the TOE highlights the complexity of stress transfers in subduction environments.

The geometrical prediction of multiple hydraulic fractures fed from a single fluid source is a major challenge due to transient stress interference among fractures and nonlinear coupling between rock deformation and fluid flow. Here we find that the evolution of multiple hydraulic fractures under isotropic stresses is controlled by a dimensionless toughness, which measures the ratio of the energy dissipated in rock fracturing to that dissipated in viscous fluid flow. The existence of a relation between the dimensionless toughness and the dimensionless length of the arrested fractures is demonstrated by using a bifurcation analysis. The numerical results show that the fractures tend to grow simultaneously in the viscosity-dominated regime, and the scaling remains effective when the number of fractures varies. This study provides a quantitative and efficient tool for predicting the fracture pattern in engineered and natural fracture systems.

The energy budget and the interplay between stable friction evolution and dynamic stick-slips are tested here under continuous slip along rough faults. We conducted 34 direct-shear experiments coupled with precise roughness measurements on diabase and limestone fault samples. The faults broad roughness ranges from highly rough and interlocked fractured interfaces to smooth polished surfaces. The analysis focuses on two slip phases: (a) the evolution of the shear strength of rough sample under stable, cumulative displacement; and (b) the dynamic of unstable stick-slip sliding. We found that the breakdown work during frictional strength evolution increases with roughness increase across multiple scales. The diabase samples are more sensitive to roughness increase than limestone samples in terms of the breakdown work implied by frictional evolution. We attribute this increased diabase sensitivity with fault roughness to its higher bulk elasticity and not to the fault shear stiffness. The diabase faults displayed multiple periodic system-size stick-slips, and the measured stored energy during the preparatory stage were surprisingly independent of the fault roughness. This finding suggests that during the preparatory stage a balance between the intracycle fault stiffness and stress drop govern the stored energy magnitude. Further, this energy balance suggests that some interface conditioning occurs before the spontaneous slip overcomes a sticking barrier.

Groundwater level tidal analysis is a powerful technique to monitor aquifer's permeability and hence its change over time. Earthquakes are known to affect aquifer's properties, in their vicinity through static stress changes but also further away through dynamic stresses. Most often changes are in the form of permeability increases, but sometimes decreases; the changes can be either permanent or transient. These observations are relatively well documented but the physical processes behind these changes are not well understood. By combining solid-earth and barometric tidal groundwater level responses in a borehole in a coherent poroelastic theoretical framework, and a bi-layer hydrogeological model, we recover a 15 years-long time series evolution of aquifer transmissivity and shear modulus. This study showcases the full potential of the tidal analysis method, coupling pore pressure diffusion and rock deformation, at the frontier of hydrogeology and rock physics. This unprecedented measurement of permeability and shear modulus evolution by tidal analysis reveals, during interseismic period, high sensitivity of this shallow aquifer to effective stress, and thus to pore pressure. Thanks to additional finite element simulation of seismic wave propagation, we explore the different mechanisms affecting permeability and shear modulus in the studied fractured andesite aquifer. This study confirms the predominant role of seismic dynamic stresses, and more precisely of dynamic shear stresses, in the change of permeability following an earthquake.

The Southwestern United States experiences active deformation, seismicity, and magmatism, remarkable in an intraplate setting. The Basin and Range and Colorado Plateau (CP) are inferred to differ in lithospheric thickness, but modeling geophysical properties of the lithosphere, in particular the depth of the Lithosphere-Asthenosphere Boundary (LAB), across the entirety of the region, has proved challenging. Here, we introduce a new model of 1-D depth profiles in shear wavespeed, determined through a probabilistic joint inversion of information from Sp receiver functions and Rayleigh wave phase velocity. From these profiles we quantify the locations and Vs contrast of wavespeed gradients that represent boundaries such as the Moho, the LAB, and intralithospheric discontinuities. We infer a lithosphere that is thinner and lower in Vs in the Basin and Range. In the CP and farther north, the LAB is more gradual, deeper, and intermittently observed. We also observe Mid-Lithospheric Discontinuities (MLDs) near the boundaries between the CP, Wyoming Craton, and Northern Basin and Range, as well as within the Craton. When both an MLD and LAB are observed, the Vs gradient associated with the LAB is narrower than expected. Finally, we image Positive Velocity Gradients beneath areas of thinner lithosphere, which are consistent with recent global observations that have been attributed to the base of a partially molten zone below the lithosphere. Overall, the picture of the lithosphere-asthenosphere system that emerges is one of considerable structural complexity with a strong dependency on tectonic regime and geological history.

Some rock and soil samples exhibit significant loss of magnetic susceptibility (χ) with increasing applied field amplitude even at relatively low (10–100s of A/m) fields, a behavior which remains unexplained. Exceptionally strong negative field-dependence of susceptibility (χ HD) is present in sandstones and altered intermediate-felsic igneous rocks in several cores from the northeastern Oklahoma subsurface. These same rocks also show elevated frequency-dependence of susceptibility (χ FD), with reasonable correlation of χ HD to χ FD, and frequency-dependent χ HD. Results from multiple characterization methods indicate that strongly negative χ HD in these rocks is linked to a yet-unidentified phase which begins the approach to magnetic saturation in low fields (<1 mT/800 A/m), shows elevated χ FD to low temperatures, is unstable at high temperatures, possesses significant anisotropy of magnetic susceptibility, and becomes paramagnetic above ∼83°C. Clear associations with fluid alteration features indicate that this material may be highly relevant to rock alteration, diagenetic, and environmental studies.

Determining the mechanisms responsible for intraplate volcanism - such as slab devolatilization melting versus active mantle plumes - remains a challenge. The greater South China Sea (SCS) region has experienced extensive intraplate Cenozoic volcanism across areas including Hainan, Southeast Indochina, northern Borneo, the northern SCS, and the post-spreading SCS basin. The prevalence of volcanism distributed widely across this region prompts fundamental questions about the key geodynamic processes driving such diverse magmatic activities. In this study, we elucidate the mantle transition zone (MTZ) discontinuities in this region using SS precursors, which helps to overcome the sparse seismic coverage due to its predominantly oceanic setting. We collected over 16,000 high-quality seismograms that sample the upper mantle and MTZ beneath this region from global earthquakes and stations. After correcting for the effects of shallow crustal variations and upper mantle heterogeneity on traveltimes of SS phases and their precursors, we unveil lateral variations in the MTZ boundaries (d410 and d660) and intricate features of the mid-MTZ reflectors (S520S). Significant MTZ thinning and normal S520S waveforms beneath Hainan provide compelling evidence for mantle upwelling through the MTZ. Conversely, the evident splitting of S520S beneath the northern SCS, Southeast Indochina, and northern Borneo, all characterized by stagnant subducted slabs, indicates that the volcanism in these regions likely originated from a mechanism distinct from the active upwelling beneath Hainan. Dehydration melting attributed to devolatilizing stagnant slabs in the MTZ is a potential cause for Cenozoic volcanism in these regions.

The East African Rift System (EARS) provides an opportunity to constrain the relationship between magmatism and plate thinning. During continental rifting, magmatism is often considered a derivative of strain accommodation—as the continental plate thins, decompression melting of the upper mantle occurs. The Turkana Depression preserves among the most extensive Cenozoic magmatic record in the rift. This magmatic record, which comprises distinct basaltic pulses followed by periods of relative magmatic quiescence, is perplexing given the lack of evidence for temporal heterogeneity in the thermo-chemical state of the upper mantle, the nonexistence of lithospheric delamination related fast-wave speed anomalies in the upper mantle, and the absence of evidence for sudden, accelerated divergence of Nubia and Somalia. We focus on the Pliocene Gombe Stratoid Series and show how lithospheric thinning may result in pulsed magma generation from a plume-influenced mantle. By solving the 1D advection-diffusion equation using rates of plate thinning broadly equivalent to those measured geodetically today we show that despite elevated mantle potential temperature, melt generation may not occur and thereby result in extended intervals of quiescence. By contrast, an increase in the rate of plate thinning can generate magma volumes that are on the order of that estimated for the parental magma of the Gombe Stratoid Series. The coincidence of large-volume stratiform basalt events within the East African Rift shortly before the development of axial zones of tectonic-magmatic activity suggests that the plate thinning needed to form these stratiform basalts may herald the onset of the localization of strain.

Shear strain localization refers to the phenomenon of accumulation of material deformation in narrow slip zones. Many materials exhibit strain localization under different spatial and temporal scales, particularly rocks, metals, soils, and concrete. In the Earth's crust, irreversible deformation can occur in brittle as well as in ductile regimes. Modeling of shear zones is essential in the geodynamic framework. Numerical modeling of strain localization remains challenging due to the non-linearity and multi-scale nature of the problem. We develop a numerical approach based on graphical processing units (GPU) to resolve the strain localization in two and three dimensions of a (visco)-hypoelastic-perfectly plastic medium. Our approach allows modeling both the compressible and incompressible visco-elasto-plastic flows. In contrast to symmetric shear bands frequently observed in the literature, we demonstrate that using sufficiently small strain or strain rate increments, a non-symmetric strain localization pattern is resolved in two- and three-dimensions, highlighting the importance of high spatial and temporal resolution. We show that elasto-plastic and visco-plastic models yield similar strain localization patterns for material properties relevant to applications in geodynamics. We achieve fast computations using three-dimensional high-resolution models involving more than 1.3 billion degrees of freedom. We propose a new physics-based approach explaining spontaneous stress drops in a deforming medium.

In order to identify relations between mechanical behavior, deformation mechanisms, microstructural properties, and H2O distribution, Tana-quartzite samples with added H2O ranging from 0 to 0.5 wt.% were deformed by axial shortening at constant displacement rates, at 900°C and 1 GPa, reaching up to ∼30% bulk strain. Samples with lower quantities of added H2O (0.1 and 0.2 wt.%) were in average ∼30 MPa weaker than the as-is samples with no added H2O. In contrast, samples with more than 0.2 wt.% added H2O revealed more variable mechanical behavior, showing either weaker or stronger trend. The weaker samples showed strain localization in their central parts in the vicinity of the thermocouple, that is, the hottest parts of the samples, whereas the stronger samples showed localization in their upper, slightly colder parts. Bulk deformation is accommodated by crystal plasticity and dissolution-precipitation processes. Distribution of H2O in our samples revealed systematic decrease of H2O content in the interiors of original grains, caused by increasing strain and H2O draining into grain boundary regions. With increasing content of added H2O, the quartz recrystallization gradually changes from subgrain-rotation-dominated to crack-induced nucleation, along with increasing quantity of melt/fluid pockets. The unexpected strain localization in the upper parts of stronger samples most likely results from mode-1-cracking that led to drainage of grain boundaries (GB) due to the crack dilatancy effect, and inhibited dissolution-precipitation in the hottest part of the samples next to the thermocouple. The locus of deformation is then shifted to colder regions where more H2O is available along GB.

We exploit nonlinear elastodynamic properties of fractured rock to probe the micro-scale mechanics of fractures and understand the relation between fluid transport and fracture aperture under dynamic stressing. Experiments were conducted on rough, tensile-fractured Westerly granite subject to triaxial stresses. We measure fracture permeability for steady-state fluid flow with deionized water. Pore pressure oscillations are applied at amplitudes ranging from 0.2 to 1 MPa at 1 Hz frequency. During dynamic stressing we transmit ultrasonic signals through the fracture using an array of piezoelectric transducers (PZTs) to monitor evolution of interface properties. We examine the influence of fracture aperture and contact area by conducting measurements at effective normal stresses of 10–20 MPa. Additionally, the evolution of contact area with stress is characterized using pressure sensitive film. These experiments are conducted separately with the same fracture and map contact area at stresses from 9 to 21 MPa. The measurements are a proxy for “true” contact area for the fracture surface and we relate them to elastic properties using the calculated PZT sensor footprints via numerical modeling of Fresnel zones. We compare the elastodynamic response of the fracture using the stress-induced changes in ultrasonic wave velocities for transmitter-receiver pairs to image spatial variations in contact properties. We show that nonlinear elasticity and permeability enhancement decrease with increasing normal stress. Additionally, post-oscillation wave velocity and permeability exhibit quick recoveries toward pre-oscillation values. Estimates of fracture contact area (global and local) demonstrate that the elastodynamic and permeability responses are dominated by fracture topology.

The array-based frequency-Bessel transform method has been demonstrated to effectively extract dispersion curves of higher-mode surface waves from the empirical Green's functions (EGFs) of displacement fields reconstructed by ambient noise interferometry. Distributed acoustic sensing (DAS), a novel dense array observation technique, has been widely implemented in surface wave imaging to estimate subsurface velocity structure in practice. However, there is still no clear understanding in theory about how to accurately extract surface-wave dispersion curves directly from DAS strain (or strain rate) data. To address this, we extend the frequency-Bessel transform method by deriving Green’s functions (GFs) for horizontal strain fields, making it applicable to DAS data. First, we test its performance using synthetic GFs and verify the correctness of extracted dispersion spectrograms with theoretical results. Then, we apply it to three field DAS ambient-noise data sets, two recorded on land and one in the seabed. The reliability and advantages of the method are confirmed by comparing results with the widely used phase shift method. The results demonstrate that our extended frequency-Bessel transform method is reliable and can provide more abundant and higher-quality dispersion information of surface waves. Moreover, our method is also adaptable for active-source DAS data with simple modifications to the derived transform formulas. We also find that the gauge length in the DAS system significantly impacts the polarity and value of extracted dispersion energy. Overall, our study provides a theoretical framework and practical tool for multimodal surface wave dispersion measurement using DAS data.

Ultrahigh-temperature (UHT; >900°C) metamorphism drives crustal differentiation and is widely recognized in the rock record, but its geodynamic causes are debated. Previous work on granulite-facies metapelite xenoliths from San Luis Potosí, Mexico suggests the lower crust experienced a protracted UHT metamorphic event that coincided with the onset of regional extension. To determine the duration, conditions, and heat sources of UHT metamorphism recorded by these xenoliths, this study characterizes the major-element, trace-element, and U-Pb isotopic systematics of quartz, rutile, feldspar, garnet, and zircon by in situ electron microprobe (EPMA) and laser-ablation inductively coupled-plasma mass spectrometry (LA-ICP-MS), and augments these data with detailed petrography, thermobarometry, phase equilibria modeling, and diffusion modeling. Thermobarometry and phase equilibria modeling suggest peak metamorphic conditions exceeded 0.7 GPa and 900°C. Zircon petrochronology confirms >15 Myr of UHT conditions since its onset at ∼30 Ma. A small population of zircon record elevated temperatures following transition from backarc compression to extension during the waning stages of orogenesis (60–37 Ma). Garnet preserves trace-element zoning and mineral inclusions consistent with suprasolidus garnet growth and subsequent compositional modification by intracrystalline rare-earth element diffusion during protracted heating, with diffusion chronometry timescales in agreement with zircon data, followed by fluid-driven remobilization of trace elements along now-healed fractures within ∼1 Myr of eruption. In sum, these data are most compatible with lithospheric mantle attenuation or removal as the dominant heat transport mechanism driving synextensional UHT metamorphism and crustal melting, which has bearing on models for crustal differentiation and formation of modern and ancient granulite terranes globally.