Lower crustal earthquake associated with highly pressurized frictional melts

Earthquakes at lower crustal depths are common during continental collision. However, the coseismic weakening mechanisms required to propagate an earthquake at high pressures are poorly understood. Transient high-pressure fluids or melts have been proposed as a viable mechanism, but verifying this requires direct in situ measurement of fluid or melt overpressure along fault planes that have hosted dynamic ruptures. Here, we report direct measurement of highly overpressurized frictional melts along a seismic fault surface. Using Raman spectroscopy, we identified high-pressure quartz inclusions sealed in dendritic garnets that grew from frictional melts formed by lower crustal earthquakes in the Bergen Arcs, Western Norway. Melt pressure was estimated to be 1.8–2.3 GPa on the basis of an elastic model for the quartz-in-garnet system. This is ~0.5 GPa higher than the pressure recorded by the surrounding pseudotachylyte matrix and wall rocks. The recorded melt pressure could not arise solely from the volume expansion of melting, and we propose that it was generated when melt pressure approached the maximum principal stress in a system subject to high differential stress. The associated palaeostress field demonstrates that a strong lower crust accommodated up to 1 GPa differential stress during the compressive stage of the Caledonian orogeny. Earthquakes in the lower crust may be facilitated by overpressure of frictional melts, according to pressure estimates from an analysis of quartz inclusions in garnets from pseudotachylytes in the Bergen Arcs.

R ecent studies demonstrate that the structural and metamorphic transformations taking place in the lower continental crust during plate convergence and orogeny are strongly coupled to earthquakes [1][2][3][4][5][6] . Field and microstructural evidence suggest that lower crustal pseudotachylytes (fault rocks interpreted as quenched frictional melts) are generated by earthquake slip in dry rocks, possibly reflecting gigapascal-level differential stress [7][8][9] . Seismic slip requires a coseismic weakening mechanism 10 . Transient high fluid or melt pressures generated due to localized shearing has been proposed as a likely cause of weakening 11 , but so far it has not been possible to constrain the pressure of the frictional melt during slip directly from observations of pseudotachylytes or other fault zone rocks. The application of phase equilibria is hindered because chemical equilibrium is often not reached during the extremely fast cooling of the melt 12 . Studies of subperpendicular injection veins into the wall rocks to constrain the conditions required for melt injection have also failed to determine the pressure of the frictional melt [13][14][15] .
In this article, we report pressure estimates of quenched frictional melts from two pseudotachylyte veins in lower crustal granulites from the Bergen Arcs, Western Norway. Abundant quartz inclusions were entrapped in dendritic garnets that grew during the early stages of melt cooling 16 . Using Raman spectroscopy, we calculated residual pressures in quartz inclusions that crystallized immediately after the earthquake. Following exhumation, the residual pressure was developed due to the different volumetric responses of quartz and garnet when subject to temperature and pressure changes 17,18 . An elastic model for the quartz-in-garnet system allows us to recover the pressure under which quartz and garnet crystallized. This provides a unique opportunity to reconstruct the palaeostress state on a fault plane immediately after a lower crustal earthquake.

Pseudotachylytes
We sampled pseudotachylytes from the Holsnøy Island, located in the Lindås Nappe of the Bergen Arcs in Western Norway. The Lindås Nappe is composed of granulite facies remnants of Proterozoic lower crust belonging to the former Jotun-Lindås microcontinent 19 . During Caledonian plate convergence between Baltica and Laurentia in the middle Silurian (approximately 430 Ma), ophiolites and Laurentian island-arc complexes were emplaced onto the Jotun-Lindås microcontinent. These compressive stages of the Scandian orogeny triggered earthquakes and subsequent fluid-induced metamorphism that locally transformed the anhydrous granulite facies rocks to eclogites and amphibolites in faults and shear zones 20 located near the leading edge of this microcontinent [21][22][23] . Maximum metamorphic conditions were estimated to 670-690 °C and 2.1-2.2 GPa 24 . However, on the basis of new geochronological and petrological data along with a simple mechanical model, it was argued that such pressures only arose locally during fluid-induced weakening of highly stressed granulites and do not represent the entire Lindås Nappe 4,8 .
At Holsnøy, eclogite-and amphibolite-facies pseudotachylytes cut pristine granulites 1 and show single rupture displacements of up to 2 m (ref. 25 ). The pseudotachylyte samples described here were collected at Ådnefjellet on the northeastern part of Holsnøy ( Supplementary Fig. 1). Two pseudotachylyte veins are studied: (1) a vein in intact wall rocks (A17) with no visible injection veins (Fig. 1a) and (2) a vein (~100 m away) in damaged wall rocks with associated injection veins (H19, Fig. 1b). Both intact and damaged wall rocks are composed of granulite facies anorthosite comprising plagioclase, pyrope-rich garnet, diopside and scapolite (meionite). Three cores were drilled from sample H19: H19-A located within the pseudotachylyte fault vein, H19-B located within the injection vein and H19-C located at the junction between the fault Lower crustal earthquake associated with highly pressurized frictional melts Xin Zhong 1,2 ✉ , Arianne J. Petley-Ragan 1 , Sarah H. M. Incel 1,3 , Marcin Dabrowski 1,4 , Niels H. Andersen 5 and Bjørn Jamtveit 1 Earthquakes at lower crustal depths are common during continental collision. However, the coseismic weakening mechanisms required to propagate an earthquake at high pressures are poorly understood. Transient high-pressure fluids or melts have been proposed as a viable mechanism, but verifying this requires direct in situ measurement of fluid or melt overpressure along fault planes that have hosted dynamic ruptures. Here, we report direct measurement of highly overpressurized frictional melts along a seismic fault surface. Using Raman spectroscopy, we identified high-pressure quartz inclusions sealed in dendritic garnets that grew from frictional melts formed by lower crustal earthquakes in the Bergen Arcs, Western Norway. Melt pressure was estimated to be 1.8-2.3 GPa on the basis of an elastic model for the quartz-in-garnet system. This is ~0.5 GPa higher than the pressure recorded by the surrounding pseudotachylyte matrix and wall rocks. The recorded melt pressure could not arise solely from the volume expansion of melting, and we propose that it was generated when melt pressure approached the maximum principal stress in a system subject to high differential stress. The associated palaeostress field demonstrates that a strong lower crust accommodated up to 1 GPa differential stress during the compressive stage of the Caledonian orogeny.
vein and the injection vein (Fig. 1b). The pseudotachylyte veins are subparallel to the wall-rock foliation ( Fig. 1b and Supplementary  Fig. 2). Wall-rock damage is asymmetrically developed and most pronounced on one side of the pseudotachylyte. This is believed to reflect the asymmetric stress distribution around the tip of the dynamic rupture preceding seismic slip 25 .
The pseudotachylyte vein in intact wall rock (A17) is ~2 cm thick, whereas the pseudotachylyte fault vein in the damaged wall rock (H19-A) is ~1 cm thick. The injection veins of H19-B are also ~1 cm thick. The pseudotachylyte matrix in A17 is made up of fine-grained omphacite (Jd 38 ), amphibole, sodic plagioclase (Ab 83 ), quartz and kyanite with minor alkali feldspar, rutile and sulfides 7 . The pseudotachylyte matrix in the H19 samples contains less omphacite, which was partly replaced by amphibole. The local equilibrium of omphacite, sodic plagioclase and quartz in the pseudotachylyte matrix constrains the ambient pressure-temperature (P-T) conditions to 1.5-1.7 GPa and 650-750 °C 26,27 (Fig. 2a, details in Supplementary Materials). This assumes that thermodynamic equilibrium is achieved upon quenching of the frictional melt. Additional evidence supporting the stability of plagioclase during seismic slip is reported by Petley-Ragan et al. 7,28 . Plagioclase aggregates filling coseismic microfractures in the wall rocks immediately adjacent to the pseudotachylyte veins show 'foam' textures, and individual grains are essentially dislocation-free in contrast to the dislocation-rich plagioclase in the wall rock 7,28 . Dislocation-free grains show a strong shape-preferred orientation (Supplementary Fig. 3) with the long axes preferentially aligned parallel to the pseudotachylyte vein regardless of the orientation of the fracture, suggesting plagioclase growth and stability during or immediately after the earthquake. Additional phases occurring within these microfractures include kyanite, clinozoisite, quartz and K-feldspar, which represent an equilibrium assemblage under the same pressure condition as the mineral assemblage found in the pseudotachylyte matrix (~1.5-1.7 GPa, Fig. 2a). This pressure is consistent with recent P-T estimates for 'recrystallized pseudotachylyte' at 670-680 °C and 1.5-1.6 GPa 24 . The preceding observations such as shape-preferred orientation and foam texture deserve special attention as they may provide evidence on the chemical equilibrium of the pseudotachylyte during or immediately after the earthquake.
Dendritic garnet is observed in all samples (H19 and A17), with a decrease in size and increase in abundance towards the wall rocks ( Fig. 3b,c). This was interpreted to reflect increasing cooling rates towards the wall rock 16,29 . Clasts of wall-rock garnets are observed in both the pseudotachylyte fault vein (Fig. 3d for H19-A) and a lateral injection vein ( Fig. 3e for H19-B). Inclusion-free garnet fragments are 'floating' in the pseudotachylyte and contain growth rims of dendritic garnet (Fig. 3e). Within the injection vein (Fig. 3e), the flow direction of frictional melt can be inferred from the redistribution of garnet fragments in the pseudotachylyte. Within the pseudotachylyte fault vein (Fig. 3a,d), flow structure is less prominent and wall-rock garnet clasts are less common, suggesting that the dendritic garnets may have crystallized from a stagnant melt after slip and flow.
In sample A17, the dendritic garnets are chemically homogeneous with an average composition near Alm 47 Pyr 32 Grs 21 . Inclusions of quartz, rutile, kyanite, omphacite, apatite and sulfides are found in the dendritic garnets (Fig. 3f). No plagioclase inclusions were observed, suggesting that they crystallized after the dendritic garnet. In sample H19, the garnet core composition is similar to sample A17, while the garnet rim is more Ca rich (Alm 44 Pyr 14 Grs 42 ). The garnet core contains only omphacite inclusions while the garnet rim contains exclusively quartz and kyanite inclusions (Fig. 3g). Mass balance calculation comparing the pseudotachylyte matrix composition with the bulk garnet, including its inclusions, suggests that initial garnet growth increased the SiO 2 content in the melt until quartz saturation, leading to subsequent entrapment of quartz inclusions in the garnet rim (Fig. 3h). Na 2 O was enriched in the melt near the garnet during its growth (Fig. 3i). By contrast, FeO was depleted around the garnet (Fig. 3j) as observed previously 16 .
The garnets are interpreted to have crystallized directly from the frictional melt on the basis of the following observations: (1) most of the garnets are dendritic and skeletal, similar to yttrium aluminium garnet formed from undercooled melt 30 ; (2) the grain size and abundance indicate increasing cooling rates towards the wall rock 16 . On the basis of mass balance calculation (Fig. 3h), chemical diffusion (particularly of the slowly diffusing species, for example, Al, Ca and Si 31 ) is not likely to be a limiting factor in controlling the rapid crystallization of dendritic garnets from the frictional melt.

High-pressure quartz inclusions
Raman spectroscopy was used to obtain the wavenumber shifts of the quartz inclusions in garnet. We measured nearly spherical, fully entrapped, isolated quartz inclusions (2-4 μm in diameter) situated at least 10 μm away from each other and the thin-section surface ( Supplementary Fig. 4). Residual quartz inclusion pressures were calculated on the basis of the wavenumber difference between pressurized inclusions and an unstrained gem-quality quartz crystal 32 ( Fig. 2b and Supplementary Fig. 5). Although the experimental calibration applied hydrostatic stress on quartz 32 , the effect of elastic anisotropy for the quartz-in-garnet system is insignificant for the application of elastic barometry 33 (see 'Method of Raman elastic barometry' in the Supplementary Information). The average quartz inclusion pressure is ~0.72 GPa for A17 and 0.64 GPa for H19 ( Supplementary Fig. 6). The maximal inclusion pressure is ~0.83 GPa for A17 and 0.74 GPa for H19. For H19, there is no clear difference in residual pressures for quartz inclusions in garnets from the pseudotachylyte fault vein (H19-A) and the injection vein (H19-B). There is no correlation between the residual pressure in a quartz inclusion and the distance of its host garnet to the wall rock. It is unlikely that plastic/viscous deformation has substantially modified the quartz pressure 34 because (1) no evidence of microfracturing around quartz inclusions has been observed, (2) the systematic residual pressure difference between A17 and H19 (located 100 m apart) cannot be explained if they had both been reset at peak P condition and (3) the measured quartz inclusion pressures are very consistent within each sample. The pressure difference between the samples is interpreted to be due to the formation of the injection vein in H19. A one-dimensional (1D) isotropic elastic model is used to recover the entrapment pressure of the quartz inclusions in the dendritic garnets 35 (Methods). To constrain the pressure, the crystallization temperature of garnet is required. Geothermometry based on the Fe-Mg exchange between dendritic garnet cores and omphacite inclusions (Supplementary Data and Fig. 3g) yields temperatures around 800 °C (ref. 36 ). However, due to the rapid quenching of frictional melt, the precise temperature of omphacite and dendritic garnet growth is uncertain; thus, additional constraints are needed. The crystallization kinetics of minerals that grow from a melt during rapid cooling is poorly understood. To our knowledge, no experimental data have been reported on the crystallization kinetics of Fe/Mg/Ca/Mn garnet growth from strongly undercooled melts. The 4D microtomography data on clinopyroxene crystallization in undercooled basaltic melt suggests that the clinopyroxene nucleation starts after a few minutes, and the crystal growth may last tens of minutes 37 . On the basis of a 1D thermal diffusion model (Methods), a lower estimate of the cooling rate of the frictional melt is calculated at the margin of a 1 cm melt vein initially at 1,500 °C (above the almandine melting point at ~1,300 °C 38 ) confined by 700 °C wall rocks (top x-axis in Fig. 2c). After 5 min, the melt temperature has dropped to ~800 °C (similar to the thermometric constraints based on omphacite inclusions in garnet).
Assuming that the dendritic garnets were fully crystallized and had gained their plastic strength after tens of seconds to a few minutes, the melt pressure estimate based on the residual quartz inclusions pressure is constrained to ~2.0-2.3 GPa for A17 and ~1.8-2.0 GPa for H19 (Fig. 2c). This pressure represents the most conservative estimate as the melt temperature cannot be lower than the ambient wall-rock temperature. At any entrapment temperature, the difference between the pressure obtained for the frictional melt and the pressure at plagioclase breakdown is always higher than 0.5 GPa (Fig. 2). If the entrapment of quartz in garnet was completed at, for example, 900 °C, the melt pressure would have been ~2.5 GPa, and the pseudotachylyte matrix would have equilibrated at a maximum pressure of ~1.8-2.0 GPa.

overpressure in frictional melt
Volume increase due to melting and heating of the wall rocks may generate melt overpressure 11 . We used Eshelby's solution to model the melt overpressure due to volume expansion (Methods) 39,40 . Using realistic isotropic elastic parameters (Fig. 4a) and assuming 10% volume expansion on melting, it is shown that low melt overpressure (<0.2 GPa) can be generated for a penny-shaped melt layer with length/thickness ratio (L/d) ~40. In the field, the pseudotachylyte veins can be traced up to 10 m, which leads to extremely high L/d (>500). For L/d = 100, the overpressure due to volume expansion is insignificant (<0.05 GPa). Therefore, it is not likely that the volume expansion during melting is the reason for the recorded overpressure. High pressure in a confined volume of melt or weak body may arise in a system where the wall rocks are subject to high differential  stress 8,41,42 . Considering that the melt has negligible shear modulus and low viscosity, the stress state in the melt is close to hydrostatic. When the confining wall rock is subject to differential stress (σ y -σ x , where σ y is the normal stress along the short axis of a penny-shaped melt layer, Fig. 4), the melt pressure must match this normal stress at the melt-wall-rock interface to keep force balance. On the basis of Eshelby's solution, the melt pressure is predicted to be higher than the confining pressure of the wall rock if the maximal principal stress is subperpendicular to the melt vein (Fig. 4b). For a differential stress of 1 GPa (σ y = 2.2 GPa, σ x = σ z = 1.2 GPa), the melt pressure exceeds 2 GPa as predicted by Raman barometry. In this case, the mean stress is ~1.5 GPa, consistent with the peak metamorphic pressure recorded in the wall rock and pseudotachylyte matrix assuming thermodynamic equilibrium. Here, we consider the mean stress to be the main controlling factor on thermodynamic equilibria under non-hydrostatic stress 43 . This mechanical set-up is similar to the model proposed by Jamtveit et al. 8 to explain eclogite formation in locally weakened rock volumes within a highly stressed granulite matrix. When injection veins form, the melt pressure P melt must be higher than the far-field stress normal to the injection veins (in this case, σ x if the injection is orthogonal to the fault vein) to inflate it 13 . In the case of sample H19, the injection vein is subperpendicular to the pseudotachylyte fault vein (Fig. 1b). For the injection vein, an elastic model can be applied to estimate σ x (ref. 13 ). To keep the injection vein open, the difference between P melt and σ x is related to the aperture/length ratio (w/L) of the injection vein 44 : (1-v) is the elastic stiffness of the wall rock calculated with shear modulus G and Poisson ratio ν. The aperture of the injection vein (H19-B) is ~1.0 ± 0.1 cm, and the length of the injection vein is ~38 cm. Assuming a wall-rock shear modulus of 40 GPa (ref. 45 ) and Poisson ratio of 0.2 (ref. 46 ), the difference between melt pressure and σ x is estimated to 0.66 ± 0.06 GPa. This represents an independent maximum value for the differential stress because no plastic yield of the host rock is taken into account 13 . Taking the previous melt pressure estimate for H19 (Fig. 2c), the far-field stress normal to the injection vein (σ x ) is in the range 1.2-1.4 GPa. This provides an additional independent constraint on the difference between the melt pressure and the far-field stress and is consistent with the equilibrium pressure recorded by the mineral assemblage that formed the pseudotachylyte matrix.

Dynamic versus tectonic stress
The average dip angle of the pseudotachylyte veins in the studied area is ~60° (Supplementary Fig. 2), consistent with previous measurements from other localities on the Holsnøy island 47,48 . Assuming that the precursory faults arise as shear fractures that follow a Mohr-Coulomb criterion, this indicates a near-vertical maximal principal stress axis before rupturing. Extrapolation of Byerlee's law would suggest a differential stress of ~3 GPa to produce ~2 GPa normal stress on the fault plane. This stress level could not have been supported by lower crustal lithologies over geological time scales 2 .  Possible triggering mechanisms, such as dehydration embrittlement or self-localization caused by thermal runaway, are judged unlikely on the basis of field and petrographic observations, such as the lack of hydrous minerals or carbonates in the host granulites and no evidence of shear predating dynamic rupture. Hence, we speculate that the observed seismic faulting may have been triggered by coseismic dynamic stress, possibly following in the wake of larger earthquakes in the shallower 'seismogenic' regime, that is, representing aftershocks 3,9 . The fault orientations are, however, to a large extent controlled by the foliation and layering of the granulitic host rocks ( Fig. 1 and Supplementary Fig. 2). A possible fault evolution scenario is schematically shown in Fig. 4c,d.
The heat released during the sliding of the wall rocks caused the formation of a thin film of frictional melt that facilitated further sliding 49,50 . Wall-rock minerals fragmented by the initial dynamic rupture or subsequent sliding of wall rocks are distributed in the frictional melt (Fig. 4d).
Coseismic dynamic stress has been invoked to explain the recorded aperture/length ratio of injection veins 13 . However, the absence of pronounced flow structure in the pseudotachylyte fault veins during garnet growth and the systematic grain-size variation (Fig. 3a) imply that the garnets crystallized in a stagnant melt. Hence, our reconstructed stress state is not likely to reflect a transient dynamic stress pulse, but rather to represent the far-field tectonic stress. This is schematically illustrated in Fig. 4e,f. The melt pressure recorded by the quartz inclusions (~2.0-2.3 GPa, Fig. 4e) is equal to the maximal principal stress acting on the wall rock. Melt injection into damaged wall rock (Fig. 4f) led to a local reduction in the melt overpressure. We hypothesize that this may be caused by a deflation effect due to the escape of the melt during the rapid quenching process. This is manifested by the systematic difference of residual quartz inclusion pressure measured in dendritic garnets from sample A17 (intact wall rock) and H19 (damaged wall rock) (Fig. 2c).
Our observations demonstrate that injection veins may form due to highly pressurized frictional melts rather than dynamic tension parallel to the fault (compare Rowe et al. 13 ). This may potentially provide a mechanism for coseismic fault weakening in systems subject to high differential stresses and is particularly relevant to explain unstable slip at lower crustal conditions. Our results also require that the lower continental crust experience differential stress levels in the range 0.5-1.0 GPa during continental collision. Melt overpressure is predicted to depend on the orientation of the slip surface relative to the relevant stress field. Hence, future work may focus on how the melt pressures depend on pseudotachylyte vein orientations.

online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41561-021-00760-x.

Methods
Raman spectroscopy. A Horiba Jobin-Yvon (T64000) Raman spectrograph located at the Department of Chemistry, the University of Oslo, was used. The following parameters were set: entrance slit and confocal pinhole of 100 microns, grating of 900 lines per mm, focal length of 64 cm and the detector was a 1,024 × 256 pixels CCD camera with 2.56 × 2.56 μm pixel size. The laser for the Raman instrument is a Spectra-Physics diode pumped Nd: YVO4 Millennia Pro SJ12 model with 532.1 nm wavelength. An Olympus LMPlanFl 50x/0.50 NA objective was used for backscattering Raman microscopy with 10 mW laser power at the sample. The acquired wavenumber range was 115 to 1,470 cm -1 . The wavenumber pitch between two nearby pixels is ~1.3 cm -1 . The spectral resolution (minimal wavenumber distance to distinguish two individual peaks) is ~5.1 cm -1 . The Gaussian-Lorentzian equation was used to fit the Raman spectra with the MATLAB least-square fitting function lsqcurvefit. The peak intensity equation is where A to F are unknowns to be fitted, B is the Raman shift and F is the background level. A Si standard has also been used to calibrate the Raman spectrograph at the beginning of the day (peak at ~520 cm -1 ).
For quartz inclusions, a window between ~455 and 490 cm -1 is taken to fit the 464 cm -1 peak. For the 128 cm -1 peak, the window is chosen between 120 and 140 cm -1 . No deconvolution is done to separate the signals between quartz and garnet because for the chosen two peaks (128 and 464 cm −1 ), the nearby garnet signal is very weak (less than 3% of quartz) and does not pose any difference if ignored. The same procedures can be found in Zhong et al. 51 .
The uncertainties of the fitted Raman band position are obtained together with the least-square fitting. By subtracting the Raman band position of inclusion from the position of a gem-quality quartz standard (for both the 464 and 128 cm −1 bands), we obtain the wavenumber shift that can be used to calculate the residual pressure. The uncertainty of the Raman band position is given in the supplementary data, where the uncertainty of Raman band position and Raman wavenumber shift compared with a quartz standard are almost identical because the uncertainty of quartz standard is typically very small. The uncertainties of the fitted Raman shifts varied depending on the widths and intensities of the Raman bands. The estimated uncertainty is ~0.2 cm -1 for the 464 and 128 cm -1 Raman bands (1σ). The gem-quality single quartz crystal was calibrated every hour. The standard deviation (1σ) of the quartz standard shift was less than 0.2 cm -1 during the measurement session, which lasts several weeks. The accuracy of the data is sufficient for dealing with the wavenumber shift of quartz inclusions at the level of ~6-7 cm −1 .

Microprobe analysis.
An electron microprobe analyser Cameca SX100 at the University of Oslo's Department of Geosciences was used for acquiring and quantifying mineral and pseudotachylyte compositions. For micrometre-scale mineral analysis, a beam current of 10-15 nA with an accelerating voltage of 15 kV was used for a 1 µm beam diameter. Bulk pseudotachylyte compositions were obtained by running 1-2 mm profiles across each vein using a 20 µm beam diameter for each point. The bulk pseudotachylyte compositions were averaged on the basis of the profile.
Elastic model. The 1D isotropic elastic model is based on Guiraud and Powell 35 : where P inc is the residual pressure of quartz inclusion, G is the shear modulus of garnet host, V 25 • C,Pinc inc is the specific volume of inclusion at room temperature and residual pressure, V Tetp,Petp inc is the inclusion volume at entrapment P-T and V 25 • C,1 bar host is the garnet host specific volume at room P-T. The almandine and pyrope PVT relationships are based on Milani et al. 52 , and grossular is from Milani et al. 53 . The composite volume of the garnet mixture is averaged on the basis of the molar percentage of each endmember. The quartz PVT relationship is based on Angel et al. 54 . The equation is nonlinear; thus, P inc is iteratively solved until convergence is reached with a residual less than 1 Pa.
As quartz is not an isotropic mineral, its residual stress after exhumation is not perfectly hydrostatic (demonstrated by the slight mismatch between the pressures determined on the basis of the 128 and 464 cm −1 shifts in Fig. 2a). However, Murri et al. 55 and Zhong et al. 33 showed that for the quartz-in-garnet system, the amount of deviatoric stress preserved in quartz inclusion at ~2.0 GPa and 700 °C is very low (<0.08 GPa). Therefore, the elastic anisotropy of quartz inclusion is not likely to substantially modify the reported frictional melt pressure. This is further supported by garnet synthesis experiment showing that an isotropic elastic model leads to errors <~0.1-0.2 GPa 56 . This is substantially lower than the melt overpressure >0.5 GPa reported here.
Another potential bias is the shape of the quartz inclusion. Although most of the quartz inclusions in dendritic garnet are non-spherical, we analysed only near-spherical inclusions (for example, in Fig. 2) for Raman analysis. This reduces the potential modification of the residual pressure due to shape dependence. In addition, finite element simulations were performed on an irregular inclusion ( Supplementary Fig. 7), which show that the deviation of residual pressure due to pronounced shape deviation from spherical is less than 0.01 GPa. This is similar to the uncertainty associated with the Raman spectroscopy and is much lower than the estimated overpressure of the frictional melt (>0.5 GPa).
Some of the quartz inclusions are close to each other and may interact elastically and affect the residual pressure. To test the effect of elastic interaction among inclusions, we performed finite element simulation to obtain the pressure distribution of two nearby quartz inclusions. This confirmed that the elastic interaction results in only ~1-2% variation in the residual pressure at the centre of each inclusion (Supplementary Fig. 8). The stress variation is concentrated mainly near the inclusion rim where the two inclusions are close to each other (but is still less than 5%). To avoid potential influence of this proximity, we always focused the laser spot (~1 μm in diameter) at the centre of the inclusions.
Eshelby's solution. The Eshelby tensor for a penny-shaped inclusion (a 1 = a2 ≫ a3, infinite ratio of length/thickness) is given by, for example, Mura 39 : where ν is the Poisson ratio of the host and a3 a1 is the ratio between the short and long axes of the penny-shaped inclusion. The equivalent eigenstrain e* follows the formula: where C and C are the host's and inclusion's stiffness tensor, respectively, s is the Eshelby's tensor, e is the true eigenstrain, equal to 10% in our simulations (thermal strain and phase transformation strain) and e ∞ is the applied far-field (anisotropic) strain. Finally, the stress in the inclusion can be calculated as follows: Details can be found in Mura 39  , where T m is the melt temperature (1,500 °C), T w is the wall-rock temperature (700 °C), a is the half-thickness of the initial temperature perturbation and κ = k cp ρ is the thermal diffusivity. The model is applicable for a flat T plateau in the middle of an infinite space (x is zero in the middle of the T plateau). The density ρ of the anorthosite rock is 2,600 kg m -3 . The heat conductivity k is 1.9 W m -1 K -1 (ref. 57 ) and the heat capacity c p is 800 J kg -1 K -1 (ref. 58 ). The thermal diffusivity κ is about 9.1 × 10 -7 m 2 s -1 . We assume that the melt has the same thermal properties as the wall rock.

Data availability
The authors declare that all the necessary data supporting the findings of this study are available in the article and its supplementary material files. Any further data are available from the corresponding author upon request.

code availability
The MATLAB code of Eshelby's solution is available from the corresponding author upon request.