Long-period microseismicity reveals cryptic earthquake-triggered fluid activity can facilitate caldera eruptions

Domain: neuroscience
Published: 2026-01-26
Editorial outcome: minor_revision

Abstract
Spatiotemporal correlations between moderate-to-large earthquakes and volcanic eruptions indicate that earthquake-induced stress changes can trigger eruptions. Triggering is rarely instantaneous, with time-dependent processes inferred to mediate magma ascent to the surface, such as fracture propagation or fluid pressurization. While various processes have been proposed, observational constraints on the specific triggering mechanisms remain limited. Here we integrate an earthquake catalogue, generated by machine-learning and template-matching techniques to enhance microseismicity detections, with geodetic data to show that seismic triggering of the 2018 eruption of Sierra Negra Volcano, Galápagos Islands, was facilitated by a cryptic phase of fluid activity. Following 13 years of near-continuous magmatic inflation, totalling >6.5 m, stress changes induced by a moment magnitude 5.4 intra-caldera earthquake on 26 June 2018 did not immediately trigger magma release from the pressurized magma reservoir. Instead, long-period earthquake swarms revealed post-seismic fluid activity along reservoir-bounding faults in the northwestern caldera, locally promoting edifice failure and facilitating magma intrusion that initiated 8 hours later. These observations demonstrate that even when a pressurized reservoir experiences significant stress perturbation, cascading processes could potentially be essential to initiate volcanic eruptions.

Introduction
Intermediate-to-large magnitude earthquakes are spatiotemporally correlated with volcanic eruptions 1–4. This correlation can be explained by earthquake-induced changes in crustal static and/or dynamic stresses, coupled with magmatic processes that increase magma overpressure within the system 4,5. Static stress changes are commonly invoked, as volcano-magmatic systems are often aligned with the local to regional stress field, parallel to the principal horizontal crustal stress (\sigma_1) (e.g., refs. 6–8). Under such stress conditions, a reduction in normal stress, or ‘unclamping’ of magmatic plumbing systems, allows for magma migration through crustal pathways, facilitating effusive to explosive eruption episodes 6,8,9,10. Dynamic stress changes, induced by the passage of seismic waves through magmatic systems, are also hypothesized to induce a host of magmatic processes that increase magma overpressure, leading to system failure and eruption 5. However, geophysical observations that shed light on the specific mechanisms underlying static and/or dynamic stress triggering of eruptions are often unavailable, perpetuating ambiguity in the processes connecting earthquake triggers and volcano eruptions.

Sierra Negra Volcano in the Galápagos Islands, Ecuador, is one of the few examples of a basaltic caldera undergoing resurgence 11–13 and offers an exceptional site to investigate the interaction between a magmatic system, volcano-tectonic structures, and the volcanic edifice. The caldera hosts an intra-caldera fault system, the Trapdoor Fault system (TDF). This fault system is likely the reactivation of structures formed during caldera formation and is hinged in the northeast (Fig. 1). Reverse motion along the TDF has uplifted the C-shaped Sinuous Ridge by ~150 m above the caldera floor in the southwestern region of the caldera. Spatially correlated with the southwestern limb of the TDF is the Minas del Azufre hydrothermal area, a solfatara with a maximum temperature of ~200 ˚C and that crosses at least two segments of the fault (Fig. 1b). Hydrothermal activity here has resulted in alteration of the lava pile and hydrothermal mineralization 14,15. Displacement on the TDF occurs by discrete earthquakes following meters-scale magmatic uplift centered in the caldera during pre-eruptive periods 13,16,17.

The three most recent eruptions in 1979, 2005, and 2018 were preceded several hours by earthquakes near and along the TDF. Specifically, the 22 October 2005 eruption began 2 hour and 48 minutes following a moment magnitude (M_w) 5.5 earthquake with surface rupture on the TDF 17,18, while the 13 November 1979 eruption initiated 1 hour and 15 minutes after a magnitude (M) 4.3 earthquake 19. Although precise constraints on the earthquake source location are lacking for the 1979 event, the TDF remain the only volcano-tectonic structure that could have hosted the earthquake, as the caldera ring fault shows no evidence of recent displacement 17. Therefore, compared to most reported cases of volcanic eruptions triggered by tectonic earthquakes, which typically act over timescales of days to weeks or longer (e.g., refs. 1, 6), the approximately one to three hour-long intervals between earthquakes on the TDF and eruption suggest near-instantaneous opening of magmatic pathways by static stress changes from TDF fault slip 12,20.

Unlike the 1979 and 2005 eruptions, the 2018 eruption of Sierra Negra was well-documented by densely distributed seismic and geodetic monitoring networks (Fig. 1), providing important high-resolution insights into the eruption initiation mechanisms. After 13 years of near-continuous magmatic inflation totalling >6.5 m, a M_w 5.4 earthquake occurred along the southern TDF at 09:15 (all times UTC) on 26 June 2018, resulting in at least 1.83 m of co-seismic vertical displacement of the Sinuous Ridge 12 (Figs. 1, 2a). Coulomb Failure Stress and normal stress changes promoted slip and/or opening on structures in the northwestern and northern segments of the TDF following the M_w 5.4 earthquake 16,21. An intense seismic swarm and contemporaneous surface deformation initiated 8 hours (~17:00) after the M_w 5.4 earthquake, indicating failure of the edifice surrounding the shallow magma reservoir and multi-directional magma propagation within the northwestern and northern caldera 12,22. The emergence of intense seismic tremor at 19:35 marked the eruption onset 12,23. In contrast to the onsets of the 1979 and 2005 eruptions, in which the eruptions occurred less than 3 hours after medium-to-large magnitude earthquakes (M>4) on the TDF, the 2018 event exhibited an 8-hour delay between the TDF earthquake and initiation of magma migration, with eruption beginning 10 hours after the TDF earthquake. This delay raises critical questions regarding mechanisms governing the activation of magmatic activity after co-seismic stress perturbations.

The initial manually-picked earthquake hypocenter catalogue 12 documented seismicity preceding and during the 2018 Sierra Negra eruption (Fig. 1b). This catalogue predominantly comprised volcano-tectonic earthquakes (VTs) characterized by relatively high-frequency energy (>5 Hz), indicating brittle failure within the volcanic edifice along the TDF 24. However, despite the intense pre- and syn-eruptive seismicity, this catalogue lacks detail during the initiation phase, as manually-picked catalogues likely miss earthquakes with local magnitude (M_L)<1.0 (i.e., microseismicity 25,26). Such microseismicity provides critical insights into magma-volcanic and/or volcano-tectonic processes that lead to eruption. This omission particularly affects detections of events with dominant low-frequency energy (1–5 Hz) and emergent phases, known as long-period earthquakes (LPs) that are closely related to magmatic and fluid activities 25,27–29.

In this study, we integrate machine-learning and template-matching techniques to generate a new high-resolution earthquake catalogue covering the period during the initiation of the Sierra Negra eruption on 26 June 2018. The resulting catalogue allows us to identify and track swarms of LP microseismicity in the northwestern corner of the caldera 2 hours after the M_w 5.4 earthquake, which culminate in failure of the edifice surrounding the magma reservoir within the same region 6 hours later. The statistics and waveform properties of these LPs strongly suggest seismicity driven by cryptic fluid activity within the northwestern caldera ring fault and TDF region. These observations reveal cascading processes that led to eruption: co-seismic stress changes initiated localized fluid-driven stress perturbations, which eventually facilitated failure of the edifice and drainage of the sill-like magma reservoir.

Results
Seismicity and deformation associated with eruption initiation
We define four sequential phases of volcanic activity on 26 June 2018 at Sierra Negra caldera based on our analysis of seismic and geodetic data (Fig. 2). These phases are: continued pre-eruptive sill inflation (uplift) and seismicity culminating in the M_w 5.4 TDF earthquake (Phase 1); post-earthquake aftershocks and pre-intrusion unrest (Phase 2); co-intrusive seismicity and deformation (Phase 3); and syn-eruptive subsidence and seismicity (Phase 4). Frequency index (FI) analysis 29–33 (see Methods) allows systematic classification of seismicity during these 4 phases into 1224 VTs and 711 LPs (Fig. 3a, Supplementary Fig. S1). Phase 1 exhibits sustained moderate rates of VTs in the vicinity of the TDF (Fig. 4a), concurrent with surface uplift at ~4 mm/day in the center of the caldera. At 09:15, a M_w 5.4 earthquake induced a minimum 1.83 m uplift of the Sinuous Ridge along the southern segment of the TDF (Fig. 2a, c). This earthquake was immediately followed by a sequence of VT aftershocks located on the southern segment of the TDF, and repeating VTs north of Minas del Azufre, in the northwestern caldera (Figs. 3, 4b; Phase 2a). Two hours later, a new population of seismicity appeared in the same region (Phase 2b), predominantly composed of LP microseismicity (Fig. 3a, d), which has not been previously reported in historical records from Sierra Negra or other Galápagos volcanoes (Fig. 1b, c). It is notable that LPs were not evenly distributed across the TDF, but spatially clustered in the northwestern caldera, north of Minas del Azufre (Figs. 1b, 3a, 4b). Throughout Phase 2, vertical displacement data from continuous Global Navigation Satellite System (cGNSS) stations across the caldera showed no significant changes above noise levels or signals that could be associated with magmatic or volcano-tectonic processes (Fig. 4e). Seismicity rates and magnitudes increased at 17:00 (Fig. 2b, c), marking failure of the edifice in the northwestern caldera. Magma evacuated from the sill and propagated bidirectionally into the northern and northwestern sectors of the edifice (Phase 3, Fig. 4c). Displacements in high-rate (30 s) cGNSS time series tracked magma migration in the northern caldera and progressive sill deflation starting at ~17:30 12 (Fig. 2a). The eruption commenced at approximately 19:40 (Phase 4), accompanied by relatively moderate rates of VTs across the caldera during the first few hours and rapid deflation of the caldera (Figs. 2, 4d).

Families of earthquake multiplets
Earthquake waveform clustering applied to the 26 June enhanced catalogue identifies 29 distinct earthquake families, each involving at least 10 multiplet earthquakes (see Methods), comprising a total of 714 earthquakes (Fig. 3b). Earthquake families are classified into LP and VT types by the average FI values of clustered events. Two short-lived (duration <20 minutes) VT families (average FI ≥ −2) were detected in the early hours of 26 June along the southern segment of the TDF, where the fault ruptured 9 hours later, resulting in the M_w 5.4 earthquake (Fig. 3). Three distinct VT families emerged concurrently and immediately after the M_w 5.4 earthquake in the northwestern corner of the caldera (Phase 2a). Approximately 2 hours into Phase 2, this VT activity evolved into 5 LP families, located near the Minas del Azufre hydrothermal system and northwestern TDF (Figs. 1b, 3c) and characterized by average M_L <1 (Phase 2b; Fig. 3d). Multiplets within families exhibited repetitive occurrences with time intervals between successive events, defined as recurrence time (T_r), spanning 3–11 minutes on average. These events indicate localized stationary sources that rapidly reactivate over time 34,35. The 5 LP families persisted for approximately 6 hours and ceased immediately with failure of the edifice surrounding the sill-like magma reservoir, marking the end of Phase 2 at 17:00 12. Three short-lived VT families subsequently appeared on the northwestern corner of the TDF at the beginning of Phase 3 (Fig. 3). Their event magnitudes increased from 0.1 to 2.3, reflecting energy release potentially from dynamic magma intrusion 36–38 (Fig. 3d). This process was followed by 4 persistent LP and 1 VT families in the same region, which are characterized by an average M_L of 2.8 (Fig. 3a, d). Throughout Phase 3, only one VT family originated on the southern TDF, persisting into Phase 4. No other earthquake families were identified on either northern or southern TDF before the eruption (Fig. 3c). The eruption commenced at ~19:40 with new VT and LP multiplet families appearing across all TDF segments, persisting at least until the end of the day, concurrent with rapid deflation of the central caldera (Figs. 2a, 3).

Repeating VTs in volcanic regions typically result from repeated brittle failure, slip, and healing of asperities 24, while the nature of volcanic LP multiplets requires non-destructive source mechanisms dominated by low-frequency energy, such as fluid activity 25,27,28 or slow rupture processes 39. Frequency content perturbations from path or site effects (i.e., intrinsic seismic attenuation/scattering in the shallow edifice 40,41) could play a role at Sierra Negra. However, we identify VT sources within 0.4 km of LPs (Supplementary Fig. S2), consistent with the typical location uncertainties for earthquakes involved in the catalogue (see Supplementary Methods). This observation suggests that the high-frequency energy deficit in LPs here is primarily related to source processes instead of path or site effects.

Cryptic fluid activity as part of the eruption initiation process
Further insight into the nature of the LP source processes can be gained from the observation that 11 LP earthquake pairs occur as anti-correlated multiplets, also known as anti-repeaters (Fig. 5c and Supplementary Fig. S3a), based on waveform similarity (normalized cross-correlations (NCC) < −0.7) and manual inspection (see Methods). These anticorrelated multiplets demonstrate systematic waveform polarity reversals, requiring rapid local stress reversals, which are commonly associated with pronounced local stress heterogeneity or fluid activities 42,43. Interestingly, these anti-correlated multiplets are only identified during Phase 2, in LP families located within the northwest caldera. Anti-correlated LPs are not identified in other phases, and anti-correlated VTs are not detected at all throughout Phases 1–4. Therefore, the unique occurrence of anti-correlated LPs in the northwestern caldera prior to magma migration provides crucial evidence for rapid local stress variations 42, most likely associated with fluid-driven processes.

Robust positive correlations between seismic moment release (M_0) and T_r (R^2=0.7, Fig. 5a) in VT families during Phases 3–4 indicate repeated brittle failure and stress reloading coupled with stick-slip instabilities 44 and fault creep 45, similar to observations at Mount Agung during the pre-eruptive phase 46. However, LP families in Phase 2b lack statistically significant M_0–T_r scaling (R^2=0.07; Fig. 5b), reflecting stress release processes distinct from classical quasi-static slip models, where fault strength recovery is proportional to the duration of stationary contact 35,44. These findings eliminate slow rupture as a viable process for generating the LP signals 39 and provide further evidence that a fluid-related process is the source of LP multiplets, such as volatile-filled crack pressurization 47,48.

Periodicity of multiplet families designates sequences with values >1 as exhibiting relatively lower variance in T_r than expected for a Poissonian (temporally random) occurrence process (see Methods). Overall, the 5 Phase 2 LP families exhibit a maximum periodicity of 1, indicating that LP multiplets occurred randomly or slightly temporally clustered in time. In contrast, 12 of 13 Phases 3–4 VT families exhibit periodicity above 1, indicating the involvement of a temporal renewal process (Supplementary Fig. S4a). The lack of significant correlations (r=0.03) between the average T_r of multiplet families and periodicity excludes seismicity rate variations as confounding factors (Supplementary Fig. S4b). Collectively, these observations demonstrate that, unlike the potential periodic stick-slip processes identified by VT families, the 5 LP families in Phase 2b are most likely the result of a cryptic phase of fluid activity preceding the edifice failure surrounding the sill-like magma reservoir.

Methods
cGNSS data analysis
Continuous GNSS data are processed following the methods in Bell et al. 12. Specifically for this study, high-rate (1 epoch/30 s) kinematic positions are estimated using the precise point positioning (PPP) method implemented in the GIPSY-OASIS II version 6.3 software 55. Kinematic analysis followed Larson and Miyazaki 56. We apply ocean loading corrections using FES2004 57, and model wet and dry tropospheric zenith delays with VMF1 mapping functions 58. We use the 30 s data to produce kinematic time series to investigate co-seismic displacements for pre- and syn-eruptive earthquakes and deformation related to magma migration. Co-seismic displacements are calculated by differencing the average of positions one hour before and after the earthquake. All data are available through the EarthScope archive (see Data Availability).

Seismic data
To produce an enhanced earthquake catalogue, we use data from 15 three-component broadband seismometers located across the Sierra Negra caldera during 26 June 2018 (Fig. 1). Specifically, station VCH1 at Volcán Chico has been operating since 2012 as part of the 6-station Galápagos permanent seismic monitoring network run by the Instituto Geofísico, Escuela Politécnica Nacional (IGEPN) 59. The 14-station IGUANA network, deployed on Sierra Negra in April 2018, recorded continuous seismic data from 9 operational stations during the 2018 eruption onset 12. In addition, the seismic station PAYG on the neighbouring island of Santa Cruz, part of the Global Seismograph Network (GSN) since 1998, documented the seismicity associated with the onset of the 2005 eruption.

Earthquake catalogue
We utilize Volpick, an EQTransformer-based phase picker for volcano seismicity developed by Zhong and Tan 33, to pick P- and S-arrival times for VTs and LPs across seismic stations on 26 June 2018. A total of 6261 P- and 5159 S-arrivals are subsequently associated into individual earthquakes through PyOcto 60 (see Supplementary Methods). With HypoInverse software 61 and HypoDD double-difference method 62, 694 earthquake hypocenters are successfully located with average location uncertainties determined to be 0.3 km and 0.4 km in horizontal and vertical dimensions, respectively (see Supplementary Methods). M_L calculations for individual earthquakes and the determination of magnitude of completeness (M_c) for the earthquake catalogue are detailed in the Supplementary Methods.

To further enhance our catalogue, we apply a template matching algorithm through EQcorrscan 63, a Python package for cross-correlation-based earthquake detection. We first select 669 templates from the machine-learning based earthquake catalogue and cross correlate them with continuous seismic data recorded on 26 June 2018 (see Supplementary Methods). This process identifies 1241 new earthquakes, which are assumed to be co-located with their template events. Following Shelly et al. 64, we then estimate M_L for template matching detections (M_Detection) based on relative magnitudes with respect to their associated templates (M_Template), using a scaling factor of 0.7 (see Supplementary Methods, Supplementary Fig. S6). With template-matching detections, the M_c decreases by 0.25 for both Phases 1–2 and 3–4 earthquakes (Supplementary Fig. S7).

Earthquake classification and clustering
Earthquakes are classified using pre-filtered vertical-channel seismograms (2–20 Hz), analyzed within 5 s time windows, starting 1 s before P arrivals at each station. We first label 694 well-located earthquakes as LP and VT by visually analyzing their velocity waveforms and spectral characteristics (see Supplementary Methods; Supplementary Fig. S1a, b). Considering potential subjectivity of manual labels 31, we reclassify earthquakes by further calculating FI of each earthquake following:
FI = log_2 (A_higher / A_lower) (1)
where A_higher and A_lower represent mean spectral amplitudes of seismic signals in high and low frequency bands at individual stations 30. A scale of log_2 provides an equal indication of the relative amount of high- and low-frequency energy within seismic signals by positive and negative FI values 32. We use frequency bands of 9–15 and 2–5 Hz as they can differentiate 5 calibration LPs and VTs within the resulting earthquake catalogue reasonably well (Supplementary Fig. S1b). To mitigate station effects and systematic path attenuation during seismic wave propagation, we take the average FI of seismic signals with signal-to-noise ratio (SNR) ≥ 2 as final estimates for individual events 31. We find that manually-labeled events exhibit bimodal FI distributions, indicating that manually-labeled LPs show systematically lower FI compared to manually-labeled VTs (Supplementary Fig. S1c). This result is consistent with previous studies at Alaskan, Hawai’ian, and Japanese volcanoes 29,31,33. We then applied the local minimum FI of −2 as a threshold to systematically classify events with FI below this value as LPs, which also signifies that low-frequency energy in LP signals is 4 times the high-frequency energy. Since only 47% of template-matching detection signals are well recorded (SNR ≥ 2) which limits FI measurement, we thereby assume that all template-matching detections FI values are equivalent to their associated templates (Fig. 3a). This assumption is reliable as we do not observe systematic artefacts when comparing the assumed FI values to 1124 detections containing measured FI (Supplementary Fig. S8). Eventually, 1214 VTs and 711 LPs are identified within the enhanced earthquake catalogue.

Earthquakes exhibiting highly similar waveform characteristics are further clustered into families. Specifically, seismic signals with SNR ≥ 2 in the 2–20 Hz frequency band are extracted from vertical channels for P and horizontal channels for S arrivals, with 5-s time windows starting 1 s before arrivals. We calculate the NCC for all event pairs at common stations, defining events as multiplets if their network-averaged NCC values exceed 0.7 34,65 (Supplementary Fig. S3). Multiplet events are further clustered into sequences if different multiplet event pairs share common events. Template-matching detections are assigned to the same sequence as their templates.

Statistical methods
To investigate the underlying mechanisms for earthquake families, we then compare the characteristics of families between LP and VT multiplets. We first quantify the periodicity of each family as the ratio of the mean to the standard deviation of T_r following Bell et al. 66. Sequences with periodicity >1 are defined as periodic, indicating a lower variance in T_r than expected for a Poissonian occurrence process.

Multiplet families originating from stress release and reloading cycles during creep on rupture patches should theoretically exhibit positive correlations between M_0 and T_r 67. To quantify this relation, we estimate M_0 of each multiplet event within families using the empirical conversion,
M_0 = 10^(1.5*M_w + 9.105) (2)
Where M_w is defined as M_w = M_L for M_L > 3 68 and M_w = 2/3*M_L + 1 for M_L ≤ 3 69. We perform linear regression on M_0–T_r relations across the 5 pre-intrusion LP families (Phase 2) and post-intrusion VT families (Phases 3–4). To ensure statistical robustness in model fitting, we perform Jackknife resampling by systematically leaving out one sample at a time and repeating the analysis across all subsets (Supplementary Fig. S9). Our analysis suggests that VT families demonstrate significant positive M_0–T_r correlations (R^2=0.7; Fig. 5a) consistent with slow-slip fault systems 70. In contrast, LP families do not exhibit statistically significant scaling (R^2=0.07; Fig. 5b).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This is a very well written and well conceived piece that is perfectly suited for Nature Communications. The observation of unusual earthquake activity in an unpredicted location during an interval between faulting and eruption is very important and yields insight into how eruptions might be nucleated from a pressurized magma chamber. This is an outstanding data set that is directly related to that problem.\n\nI am not an experienced seismologist, but I found the explanations in the manuscript to be very clear, especially the changes in spectral signal. The authors do a good job of laying out all of the evidence that the unusual timing and frequency contents are likely driven by stresses in a porous water-saturated medium.\n\nI have two suggestions that might improve the article:\n\n1. A brief review of the hydrothermal system, especially the spatial distribution in relation to the seismic foci is in order. The best information comes from:\nGoff, F., McMurtry, G.M., Counce, D., Simac, J.A., Roldán-Manzo, A.R. and Hilton, D.R., 2000. Contrasting hydrothermal activity at Sierra Negra and Alcedo volcanoes, Galapagos Archipelago, Ecuador. Bulletin of Volcanology, 62(1), pp.34-52.\nAiuppa, A., Allard, P., Bernard, B., Lo Forte, F.M., Moretti, R. and Hidalgo, S., 2022. Gas leakage from shallow ponding magma and trapdoor faulting at Sierra Negra volcano (Isabela Island, Galápagos). Geochemistry, Geophysics, Geosystems, 23(2), p.e2021GC010288.\n\n2. I find the conclusions to be a bit vague, and to really get at the point the authors were making, I only figured it out by careful look at Figure 6. The text would be improved with a rhetorical summary. Something like:\n1. Pressure from the subcaldera sill caused inelastic failure along the TDF.\n2. Transient stress changes due to this slip caused pore pressure changes 2 km to the NW of the epicenter.\n3. The pore pressures weakened the roof of the sill (or increased deviatoric stresses by reducing normal component), triggering the propagation of a dike from the sill.\n4. How all this relates to the seismic measurements (what processes cause what signal)\n\nI made a number of less-important suggestions as comments on the manuscript pdf.", "reviewer_2": "This is an interesting paper from an incredible set of observations. The technical work looks to be high quality and the paper is generally well organized and well written. While the interpretation that local fluid processes could help promote failure and magma intrusion is very reasonable, I don’t think it is well supported. There is no attempt to understand the source mechanism for the earthquakes and the separation into distinct VT and LP categories oversimplifies the complexity here. Given that, I would argue that there is not enough support for the idea of a cryptic fluid being a part of the lead-up to the eruption. In short, I think this paper has great potential and with additional analysis, it could become a much more impactful contribution.\n\nTo begin with, the FI data for the full dataset (Fig 3a) does not appear bimodal , which might make it easier to interpret VT and LP separately, but is instead continuous across the range. It would be interesting to see a histogram like that in Figure S1b for the full catalog to see. And even the data presented in S1b are weakly bimodal. With this dataset, it is very hard to determine which are VT and LP and so interpreting them as due to brittle failure or fluid involvement is very challenging.\n\nFirst-motion mechanisms would go a long way toward helping explain the cause of the seismicity across the FI range. I realize this is a lot of additional work, but I see it as essential. It may be that many of the events with FI < -2 still have mixed first motions that could be explained by a double-couple mechanism. Furthermore, understanding the sense of slip on these faults would allow for comparison with the FEM modeling from Gregg et al.\n\nWhile the paper is generally well written, there are sections of the paper that are not clear to someone who isn’t familiar with all the literature cited about Sierra Negra volcano. It would be relatively easy to correct this.\n\nFor example, could you be more clear in the introduction about the TDF mechanism. Which way does the ‘door’ open? The trap door faults are not plotted on a map.\n\nThe term ‘sill failure’ is used multiple times. In descriptions of the sequence of events, I think this may not be appropriate since it is an interpretation. The idea of a large sill (e.g., Gregg et al., Bell et al,) is a model and is somewhat less definite that an earthquake hypocenter (also a model), so to refer to the sill failure as an observation, rather than an interpretation, seems incorrect. It is fine in the discussion, but what does ‘sill failure’ mean precisely?\n\nThe waveform shown in Figure 5 is interesting. It certainly looks to be dominated by relatively low frequencies. Based on the P wave, it seems like a nodal arrival. The black waveform has the picked P wave on the Z channel later than the first break on the N channel. Are these picks from EQTransformer or some other automated process? Also, it would make sense to remove the numbers on the vertical axis for part c since they are meaningless.\n\nIn the discussion of phase 2 beginning on line 269, the authors describe an absence of deformation and intense VTs. This does not seem well supported although it may be true. In Figure 4, that period is the only time with VTs in the caldera – though I agree there aren’t many – and the GNSS data show deformation across all 4 stations plotted. The problem is that the displacement data are normalized so one can’t tell what the scale of deformation is. I assume the reason for normalizing is to fit them all on one plot and show the consistency. In this case it would be helpful to have 4 scale bars showing the range of deformation represented by each of the stations.\n\nIn Figure 2, could you mention why there is missing data for GV03?\n\nRegarding the caption for Figure 4, Phase 3 is mentioned as an ‘eruption episode’, but the eruption did not begin until 19:40 which marks the start of Phase 4\n\nIn Figure 5, the caption for part a should start with Seismic moment, not Moment magnitude\n\nLine 212: I would argue that a R-squared of 0.07 is not a weak correlation, but is not a correlation at all\n\nLine 245: I suggest a better term than ‘near instantaneous’ for something close to 3 hours", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "We have enhanced the text according to the comments raised (below), including clarifying terminologies and conclusions, to make the details of this study more readable, comprehensible, and robust.\n\n1. A brief review of the hydrothermal system, especially the spatial distribution in relation to the seismic foci is in order. The best information comes from:\nGoff, F., McMurtry, G.M., Counce, D., Simac, J.A., Roldán-Manzo, A.R. and Hilton, D.R., 2000. Contrasting hydrothermal activity at Sierra Negra and Alcedo volcanoes, Galapagos Archipelago, Ecuador. Bulletin of Volcanology, 62(1), pp.34-52.\nAiuppa, A., Allard, P., Bernard, B., Lo Forte, F.M., Moretti, R. and Hidalgo, S., 2022. Gas leakage from shallow ponding magma and trapdoor faulting at Sierra Negra volcano (Isabela Island, Galápagos). Geochemistry, Geophysics, Geosystems, 23(2), p.e2021GC010288.\n\nResponse: We have expanded the introduction of the Minas del Azufre hydrothermal system at Sierra Negra and its spatial correlation with the Trapdoor Falt system (TDF) on lines 73-77: “Spatially correlated with the southwestern limb of the TDF is the Minas del Azufre hydrothermal area, a solfatara with a maximum temperature of ~200 ˚C and that crosses at least two segments of the fault (Fig. 1b). Hydrothermal activity here has resulted in alteration of the lava pile and hydrothermal mineralization (Aiuppa et al., 2021; Goff et al., 2000).”\n\n2. I find the conclusions to be a bit vague, and to really get at the point the authors were making, I only figured it out by careful look at Figure 6. The text would be improved with a rhetorical summary. Something like:\n- 1 -\n1. Pressure from the subcaldera sill caused inelastic failure along the TDF.\n2. Transient stress changes due to this slip caused pore pressure changes 2 km to the NW of the epicenter.\n3. The pore pressures weakened the roof of the sill (or increased deviatoric stresses by reducing normal component), triggering the propagation of a dike from the sill.\n4. How all this relates to the seismic measurements (what processes cause what signal)\n\nResponse: Thank you for your constructive suggestion to clarify the key conclusions. We have revised the Discussion section (lines 300-307) to concisely outline the logical links between fluid activity and associated processes: “The source-dominated low-frequency energy, unique presence of anti-correlated multiplets, non-periodic interevent times, and absence of M0-Tr scaling within these Phase 2 LP families, collectively indicate localized non-destructive processes most readily explained by fluid activity. High-temperature volatiles or fluids could change the friction coefficient along faults and fractures surrounding the reservoir by increasing pore fluid pressure (Caricchi et al., 2021) and promoting edifice failure and magma intrusion in the same area. We therefore propose that this cryptic phase of fluid activity plays an important role in weakening the volcanic system, facilitating edifice instability and eventually triggering the 2018 Sierra Negra eruption”. The revised text emphasizes fluid activity as the primary driver of the observed seismic signals, while highlighting its cascading effects that contributed to triggering the 2018 eruption.\n\nI made a number of less-important suggestions as comments on the manuscript pdf.\n\nLine 47: I believe the pioneering paper is Nakamura (1977)\nResponse: We have included this citation in the manuscript (line 48).\n\nLine 73: Perhaps clarify that the uplift is not centered on the TDF.\nResponse: We have clarified this point on lines 77-79: “Displacement on the TDF occurs by discrete earthquakes following meters-scale magmatic uplift centered in the caldera during pre-eruptive periods (Bell et al., 2021a; Chadwick et al., 2006; LaFemina et al., 2025)”.\n\nLine 79: Kind of vague and I’m skeptical: what regional structures and how do they point to the TDF?\n- 2 -\nResponse: As reflected in the revised text (lines 84-86): “Although precise constraints on the earthquake source location are lacking for the 1979 event, the TDF remains the only volcano-tectonic structure that could have hosted the earthquake, as the caldera ring fault shows no evidence of recent displacement (Chadwick et al., 2006)”, we have further specified the volcano-tectonic structures. The only fault structures at Sierra Negra capable of producing an earthquake of ~M5 are associated with the caldera. There is no indication that the outer ring fault has been active in the recent geologic history of the volcano.\n\nLine 80: “the hours-long intervals between earthquakes on the TDF and eruption”\nResponse: We have revised the text to further specify the time intervals on lines 89-90: “…the approximately one to three hour-long intervals between earthquakes on the TDF and eruption…”.\n\nLine 91: “failure along the margins of…” (a sill can’t really fail).\nLine 177: “… failure of the sill” – see previous comments.\nLine 273: Again, I don’t like this terminology, but perhaps it is ingrained. The sill never really fails. The wall rock does.\nLine 275: change “failure” to “propagation”.\nLine 262: In other words, the eruption was initiated by a sequence of mechanical disruption of the wall rocks instead of within the magma body itself.\nResponse: Since all these comments pertain to the terminology of \"sill failure\" used in our manuscript, we have consolidated them here for a unified response. We agree that “sill failure” can be potentially misleading. In reality, the sill itself doesn’t fail, instead, failure occurs in the edifice surrounding the shallow magma reservoir, which enables magma migration. Therefore, we have replaced the term “sill failure” with more precise phrases, such as “failure of the edifice surrounding the shallow magma reservoir” (line 101) and “failure of the edifice in the northwestern caldera” (line 150). Accordingly, we deleted \"failure\" from the original line 275 and changed it to \"magma migration initiation\" (line 299) to emphasize the spatial correlations between pre-intrusion LP multiplet families and magma initiation processes. We have reviewed the relevant terminology throughout the entire manuscript to ensure no ambiguity regarding \"sill failure\" remains.\n- 3 -\nLine 94-96: Perhaps restate the delay intervals for 79 and 05 for comparison\nResponse: Text on lines 103-106 has been revised: “In contrast to the onsets of the 1979 and 2005 eruptions, in which the eruptions occurred less than three hours after medium-to-large magnitude earthquakes (M > 4) on the TDF, the 2018 event exhibited an 8-hour delay between the TDF earthquake and …”.\n\nLine 245: I know it’s semantic but I wouldn’t call hours “near-instantaneous”\nResponse: The term ‘near-instantaneous’ here refers to the migration of magma from the sill in response to the triggering earthquake. We agree that “hours-interval” might not align perfectly with the strict meaning of “near-instantaneous” when evaluating the temporal link, at the first glance. However, as reflected in the text “…global spatiotemporal correlations between tectonic earthquakes and volcanic eruptions, with the time interval from the earthquake to eruptions spanning from hours to years…” (lines 259-260), from the perspective of the complete triggering process (i.e., tectonic earthquake → stress perturbation → edifice failure/magma migration → eruption), most such triggering processes typically operate on timescales of days to weeks or longer (e.g., Chesley et al., 2012; Bebbington & Mazocchi, 2011). Hence, by comparison, the approximately 3-hour windows between the M >4 earthquake and volcanic eruptions at Sierra Negra in 1979 and 2005 would qualify as “near-instantaneous”.\n\nAs with the eruption in 2018, we envisage that a period of a few hours was subsequently required for that magma then to ascend the 2 km to the surface. To address the semantic concern, we have carefully clarified this definition in lines 87-91 (“…compared to most reported cases of volcanic eruptions triggered by tectonic earthquakes, which typically act over timescales of days to weeks or longer (e.g., Bebbington & Mazocchi, 2011; Chesley et al., 2012), the approximately one to three hour-long intervals between earthquakes on the TDF and eruption suggest near-instantaneous opening of magmatic pathways by static stress changes from TDF fault slip (Bell et al., 2021b; Gregg et al., 2018)”) and lines 265-269 (“… eruptions of Sierra Negra Volcano were preceded less than 3 hours (i.e., 75 minutes and 168 minutes) by M 4.3 and Mw 5.5 earthquakes on the TDF, indicating near-instantaneous eruption triggering by static stress changes on the volcano-tectonic and magmatic systems (i.e., unclamping of faults and/or opening\n- 4 -\nof magmatic conduits, magma migration from the sill-like magma reservoir at 2 km depth and eruption initiation)”).\n\nBebbington, M. S., & Marzocchi, W. (2011). Stochastic models for earthquake triggering of volcanic eruptions. Journal of Geophysical Research: Solid Earth, 116(B5), doi: 10.1029/2010JB008114.\nChesley, C., LaFemina, P. C., Puskas, C., & Kobayashi, D. (2012). The 1707 Mw8.7 Hoei earthquake triggered the largest historical eruption of Mt. Fuji. Geophysical Research Letters, 39(24), doi: 10.1029/2012GL053868.\n\nLine 279: I don’t see how it is “stochastically activated”. It was activated by transient strain related to the quake\nResponse: We agree that the underlying process (i.e., fluid activity) was initially activated by the transient strain from the M 5.4 earthquake. However, the varying recurrence times of LP multiplets within identified families (with periodicity < 1) suggest that these activated processes subsequently drove LPs to recur stochastically. To address potential ambiguity, we have removed the expression and revised lines 301-303 as: “…and absence of M0-Tr scaling within these Phase 2 LP families, collectively indicate localized non-destructive processes most readily explained by fluid activity”.\n\nLine 286: Gregg et al.\nResponse: Corrected (line 309).\n\nLine 295: Clarify. There was increased pressurization over the course of years, but no sudden pressure change leading to eruption.\nResponse: We have revised the sentence on lines 318-320 to “…pre-eruptive processes devoid of increased pressurization in the magma reservoir immediately prior to the 2018 Sierra Negra eruption offer critical insights into volcano-tectonic interactions”. Thank you for clarification.\n\nLine 297: A bit awkward. Maybe “The findings discussed above…”\n- 5 -\nResponse: Text on line 321 has been revised as suggested: “…the findings discussed above collectively…” .\n\nReviewer #2 (Remarks to the Author):\nThis is an interesting paper from an incredible set of observations. The technical work looks to be high quality and the paper is generally well organized and well written. While the interpretation that local fluid processes could help promote failure and magma intrusion is very reasonable, I don’t think it is well supported. There is no attempt to understand the source mechanism for the earthquakes and the separation into distinct VT and LP categories oversimplifies the complexity here. Given that, I would argue that there is not enough support for the idea of a cryptic fluid being a part of the lead-up to the eruption. In short, I think this paper has great potential and with additional analysis, it could become a much more impactful contribution.\n\nResponse: We appreciate the positive evaluation of the dataset, the technical work, and the manuscript. However, we respectfully disagree with two of the comments: 1) that the interpretation is not ‘well supported’; and 2) that source mechanisms would help the interpretation. To address these concerns, we have undertaken additional analyses, and revised the text, specifically with reference to other studies that support our approach and interpretation. We have added detailed explanations of methodologies in this letter and further interpretations in the manuscript, ensuring all the relevant points have been carefully treated.\n\nTo begin with, the FI data for the full dataset (Fig 3a) does not appear bimodal , which might make it easier to interpret VT and LP separately, but is instead continuous across the range. It would be interesting to see a histogram like that in Figure S1b for the full catalog to see. And even the data presented in S1b are weakly bimodal. With this dataset, it is very hard to determine which are VT and LP and so interpreting them as due to brittle failure or fluid involvement is very challenging.\n\nResponse: We agree that FIs in Fig. 3a exhibit continuous distribution of values. A pronounced bimodal FI distribution is not a reasonable expectation for dataset containing ‘LPs’ and ‘VTs’, nor is it necessary for defining event categories or critical to our interpretations. The frequency content of volcanic seismicity commonly displays a wide spectrum of properties, which depend on varying combinations of source processes and path effects. Similarly continuous FI distributions have been described\n- 6 -\nat Augustine volcano (Buurman & West, 2010) and Mount Agung (Wellik et al., 2021). Meanwhile, previous studies from volcanic settings in Japan (Zhong and Tan, 2024), Hawai’i (Matoza et al., 2014), and Alaska (Song et al., 2025, Augustine volcano: Buurman & West, 2010; Mount Agung: Wellik et al., 2021) have also used this \"weakly-bimodal\" FI distributions to categorize LPs and VTs. Key to this approach is the consistency between FI-based categorization and manual labelling. As noted by Matoza et al. (2014), volcanic earthquakes can be misclassified if one relies solely on subjective visual inspections of waveform characteristics. Therefore, although FI distributions are rarely distinctively bimodal, they provide an established, systematic, and objective means for earthquake classifications.\n\nOne of the crucial observations underpinning our interpretation of local fluid activation is that the seismicity in the northwestern corner of the caldera during Phase 2b has anomalously low FI values. In response to the question about FI values for the full catalogue, we have extended this analysis to include one month of seismicity before the eruption (published in Bell et al., 2021b). These data show rare occurrences of events with FI < -2, but with no spatial or temporal pattern or clustering (Figure R1). We would not be surprised if fluid-driven seismicity occurs in other locations or times at Sierra Negra, but this does not reduce the significance of seismic clusters in Phase 2b during the eruption initiation process.\n\nIn the revised manuscript, we have added sample waveforms of typical LPs and VTs in Fig. S1a (attached below) to help convey the FI analysis more intuitively, and revised the method on lines 386-389 (“We find that manually-labeled events exhibit bimodal FI distributions, indicating that manually-labeled LPs show systematically lower FI compared to manually-labeled VTs (Fig. S1c). This result is consistent with previous studies at Alaskan, Hawai’ian, and Japanese volcanoes (Matoza et al., 2014; Song et al., 2025; Zhong & Tan, 2024)”). We hope this can help address your concerns on the FI analysis.\n- 7 -\nFigure R1. Spatio-temporal distributions of FI for earthquakes that occurred one month before the day of the 2018 Sierra Negra eruption. (a) Temporal distribution of earthquakes (circles) that occurred one month before the eruption. Circles are colored by their FI values. The horizontal dashed line marks the FI threshold (-2) applied in this study to identify LPs. (b) Spatial distribution of earthquakes (circles) that occurred one month before the eruption. Circles are colored by their FI values. Gray circles are the seismicity identified by Bell et al., (2021b) at Sierra Negra caldera on 26 June 2018.\n\nSupplementary Figure 1. (a) Sample waveforms of a calibration long-period earthquake (LP; red) and volcano-tectonic earthquake (VT; blue) recorded by the vertical channel of station SN05. Waveforms are filtered within the 2-20 Hz frequency band. To enhance microseismicity classification, we selected 5 largest-magnitude manually-labeled LPs and 5 largest-magnitude manually-labeled VTs (all magnitudes < 3) as calibration events. (b) Stacked spectra of 5 manually classified calibration LPs (red) and VTs (blue). Spectra for each event are calculated across stations with high-quality recordings (signal-to-noise ratio ≥ 2). The low and high frequency bands applied for frequency index (FI) calculation are highlighted by red and blue shaded regions, respectively. (c) Histograms of FI values for manually classified earthquakes, with red and blue colors representing LPs and VTs, respectively. The vertical dashed line indicates local minimum at FI = -2, defining the threshold for systematic earthquake classification.\n- 8 -\nBell, A. F., La Femina, P. C., Ruiz, M., Amelung, F., Bagnardi, M., Bean, C. J., Bernard, B., Ebinger, C., Gleeson, M., Grannell, J., Hernandez, S., Higgins, M., Liorzou, C., Lundgren, P., Meier, N. J., Möllhoff, M., Oliva, S.-J., Ruiz, A. G., & Stock, M. J. (2021). Caldera resurgence during the 2018 eruption of Sierra Negra volcano, Galápagos Islands. Nature Communications, 12(1), 1397, doi: 10.1038/s41467-021-21596-4.\nBuurman, H., & West, M. (2010). Seismic precursors to volcanic explosions during the 2006 eruption of Augustine Volcano. The 2006 Eruption of Augustine Volcano. 1769.\nMatoza, R. S., Shearer, P. M., & Okubo, P. G. (2014). High-precision relocation of long-period events beneath the summit region of Kı̄lauea Volcano, Hawai‘i, from 1986 to 2009. Geophysical Research Letters, 41(10), 3413-3421, doi: 10.1002/2014gl059819.\nSong, Z., & Tan, Y. J. (2025). Characteristics of deep long-period earthquakes at Alaska volcanoes from 2005 to 2017. Journal of Geophysical Research: Solid Earth, 130, doi: 10.1029/2024JB030444.\nWellik, J. J., Prejean, S. G., & Syahbana, D. K. (2021). Repeating Earthquakes During Multiple Phases of Unrest and Eruption at Mount Agung, Bali, Indonesia, 2017. Frontiers in Earth Science, 9, doi: 10.3389/feart.2021.653164.\nZhong, Y., & Tan, Y. J. (2024). Deep-Learning-Based Phase Picking for Volcano-Tectonic and Long-Period Earthquakes. Geophysical Research Letters, 51(12), doi: 10.1029/2024GL108438.\n\nFirst-motion mechanisms would go a long way toward helping explain the cause of the seismicity across the FI range. I realize this is a lot of additional work, but I see it as essential. It may be that many of the events with FI < -2 still have mixed first motions that could be explained by a double-couple mechanism. Furthermore, understanding the sense of slip on these faults would allow for comparison with the FEM modeling from Gregg et al.\n\nResponse: We understand why focal mechanism solutions (FMS) could be considered as an additional source of evidence to interpret the nature of the seismicity in Phase 2b. We did attempt to construct focal mechanisms for these earthquakes, but chose not to include them in the original submission. The low absolute magnitude of the events, low\n- 9 -\nsignal-to-noise ratio (SNR), sparse station coverage, and emergent P-arrivals mean that the solutions are poorly constrained. However, we also argue that FMS are not particularly diagnostic of the underlying source process.\n\nThe average magnitude of LPs identified by FI in Phase 2 is 0.54, which is substantially lower than that of LPs used for focal mechanism analysis in previous studies. For example, Wech et al. (2020) stacked waveforms of LPs with magnitudes ranging from 1.2 to 1.5 to derive reliable focal mechanisms. Meanwhile, Bell et al. (2021) selected 80 events around the caldera with well-recorded P-wave polarities from at least 7 stations over a period of several months, but in our study, only 31 (8%) of the LPs in Phase 2 meet this threshold. Of these 31 LPs, only 22 are associated with 3 LP multiplet families in Phase 2b.\n\nIn response to these comments, we renewed our efforts to derive the FMS for these 22 LPs using P-wave polarities determined by QuakeMigrate (Winder et al., 2025). We then applied MTfit, a Bayesian approach to seismic moment tensor inversion (Pugh & White, 2018). We first analysed the source assuming fully double-couple (DC) mechanisms, and the results (Figure R2) suggest that, due to the challenges related to station coverage, the uncertainties of these solutions are poorly constrained in terms of the fault plane determination. Specifically, the maximum probability of DC solutions for the 22 LPs averages 0.1% with an averaged Kullback-Leibler divergence of 5.3. While the divergence provides a moderate concentration of the probability density function relative to a non-informative reference distribution, uncertainties of the derived FMS persist. We then analysed the data in the full moment tensor space. The Hudson plot (Figure R3) further reveals substantial uncertainties associated with DC mechanism assumptions, and specific source mechanisms remain poorly constrained.\n- 10 -\nFigure R2. Focal mechanism solutions of LP multiplets that are well-recorded (SNR ≥ 2) by at least 7 stations within families No. 5, 6, and 9 during phase 2b. Stations with signals showing positive and negative polarities are indicated by white circles and black triangles, respectively.\n- 11 -\nFigure R3. Probability density functions of source types for 22 LP multiplets, which are plotted on the Hudson plot (Hudson et al. 1989). These LP multiplets were well-recorded (SNR ≥ 2) by at least 7 stations, belonging to families No. 5, No. 6, and No. 9 during phase 2b. Several source types are shown: compensated linear vector dipole (CLVD) sources, implosive and explosive sources, opening and closing tensile crack sources (TC and TC ), and double-couple source (DC) sources. Light yellow indicates\n+ −\nhigh probability, and dark purple indicates low (non-zero) probability.\n\nWe also tried to stack waveforms within each family to increase the SNR following Wech et al. (2020). However, we found that stacking waveforms did not consistently improve the SNR, particularly associated with the initial arrivals. In summary, the small amplitudes and emergent phase arrivals as well as poor focal sphere coverage render the FMS too poorly constrained to derive a reliable result.\n\nFurthermore, we consider that even if they were available, reliable FMS would not be definitively diagnostic of whether the LP sources are fluid related. LP sources are known to exhibit both DC and non-double couple mechanisms (e.g., Oikawa et al., 2019;\n- 12 -\nWech et al., 2020). Even for LPs with well constrained DC mechanisms, their potential association with the local stress field does not exclude the involvement of fluid activity in their source process (e.g., Oikawa et al. 2019).\n\nWech, A. G., Thelen, W. A., & Thomas, A. M. (2020). Deep long-period earthquakes generated by second boiling beneath Mauna Kea volcano. Science, 368(6492), 775-779.\nBell, A. F., La Femina, P. C., Ruiz, M., Amelung, F., Bagnardi, M., Bean, C. J., Bernard, B., Ebinger, C., Gleeson, M., Grannell, J., Hernandez, S., Higgins, M., Liorzou, C., Lundgren, P., Meier, N. J., Möllhoff, M., Oliva, S.-J., Ruiz, A. G., & Stock, M. J. (2021). Caldera resurgence during the 2018 eruption of Sierra Negra volcano, Galápagos Islands. Nature Communications, 12(1), 1397, doi: 10.1038/s41467-021-21596-4.\nPugh D. J., & White R. S. (2018). MTfit: A Bayesian Approach to Seismic Moment Tensor Inversion. Seismological Research Letters, 89(4), 1507–1513, doi: 10.1785/0220170273.\nHudson, J. A., R. G. Pearce, and R. M. Rogers (1989), Source type plot for inversion of the moment tensor, J. Geophys. Res., 94(B1), 765–774, doi:10.1029/JB094iB01p00765.\nOikawa, G., Aso, N., & Nakajima, J. (2019). Focal mechanisms of deep low-frequency earthquakes beneath Zao volcano, Northeast Japan, and relationship to the 2011 Tohoku earthquake. Geophysical Research Letters, 46, 7361–7370, doi: 10.1029/2019GL082577.\nWinder, T., Bacon, C., Smith, J., Hudson, T., Greenfield, T., & White, R. (2020). QuakeMigrate: a Modular, Open-Source Python Package for Automatic Earthquake Detection and Location, doi: 10.17863/CAM.83083.\n\nWhile the paper is generally well written, there are sections of the paper that are not clear to someone who isn’t familiar with all the literature cited about Sierra Negra volcano. It would be relatively easy to correct this.\n\nFor example, could you be more clear in the introduction about the TDF mechanism. Which way does the ‘door’ open? The trap door faults are not plotted on a map.\n- 13 -\nResponse: We have added text to the Introduction (lines 70-79) that includes further details on the description of Sierra Negra volcano: “The caldera hosts an intra-caldera fault system, the Trapdoor Fault system (TDF). This fault system is likely the reactivation of structures formed during caldera formation and is hinged in the northeast (Fig. 1). Reverse motion along the TDF has uplifted the C-shaped Sinuous Ridge by ~150 m above the caldera floor in the southwestern region of the caldera. Spatially correlated with the southwestern limb of the TDF is the Minas del Azufre hydrothermal area, a solfatara with a maximum temperature of ~200 ˚C and that crosses at least two segments of the fault (Fig. 1b). Hydrothermal activity here has resulted in alteration of the lava pile and hydrothermal mineralization (Aiuppa et al., 2021; Goff et al., 2000). Displacement on the TDF occurs by discrete earthquakes following meters-scale magmatic uplift centered in the caldera during pre-eruptive periods (Bell et al., 2021a; Chadwick et al., 2006; LaFemina et al., 2025)”.\n\nWe have also added the Trapdoor Fault in Figure 1b as requested (attached below).\n\nFigure 1. (a) Map view of seismic stations (inverted triangles) within and around the Sierra Negra caldera, Galápagos Islands, used in this study. (b) The distribution of seismic events (circles) recorded on 26 June 2018 (Bell et al., 2021b), the Trapdoor Fault (dashed black lines), the seismic stations (inverted triangles), and continuous\n- 14 -\nGlobal Navigation Satellite System (GNSS) stations (white squares) are overlain on topography (Tozer et al., 2019). The blue star marks the epicenter of the moment magnitude (M ) 5.4 earthquake that occurred at 9:15 UTC on 26 June 2018 (Bell et al., 2021b). The circles and the blue star are color-coded by event origin time in UTC. The first eruptive fissures were located north-northwest of Volcán Chico (Vasconez et al., 2018). (c) Vertical profile of earthquake hypocenters along a North-South cross-section. The sill-like magma reservoir is highlighted by the red-shaded zone (Bell et al., 2021b).\n\nThe term ‘sill failure’ is used multiple times. In descriptions of the sequence of events, I think this may not be appropriate since it is an interpretation. The idea of a large sill (e.g., Gregg et al., Bell et al,) is a model and is somewhat less definite that an earthquake hypocenter (also a model), so to refer to the sill failure as an observation, rather than an interpretation, seems incorrect. It is fine in the discussion, but what does ‘sill failure’ mean precisely?\n\nResponse: Based on comments from Reviewer 1, we have deleted the terminology ‘sill failure’ in the revised manuscript and have opted for more precise phrases, such as “failure of the edifice surrounding the shallow magma reservoir” (line 101) and “failure of the edifice in the northwestern caldera” (line 150), to avoid conceptual ambiguity and ensure the rigor of wording. We have reviewed the relevant terminology throughout the entire manuscript to ensure no ambiguity regarding \"sill failure\" remains.\n\nThe waveform shown in Figure 5 is interesting. It certainly looks to be dominated by relatively low frequencies. Based on the P wave, it seems like a nodal arrival. The black waveform has the picked P wave on the Z channel later than the first break on the N channel. Are these picks from EQTransformer or some other automated process? Also, it would make sense to remove the numbers on the vertical axis for part c since they are meaningless.\n\nResponse: As reflected in the text, “We utilize Volpick, an EQTransformer-based phase picker for volcano seismicity developed by Zhong and Tan (2024), to pick P- and S-arrival times…” (lines 354-355). The P- and S-wave picks annotated in Fig. 5c are derived directly from the phase picker. These waveforms are then aligned based on time shifts from cross-correlation analysis between waveform pairs.\n\nWe understand that the first break of the P-wave on the HHZ channel (black waveform) appears later than that on the HHE/N channels, which raised the question of a nodal\n- 15 -\narrival. This discrepancy arises from the time alignment method applied in the figure: to highlight the characteristics of different phases, the vertical channel (HHZ) was aligned relative to the P-wave arrival while the horizontal channels were aligned relative to the S-wave arrival. Thus, the \"0 s\" on these channels do not represent the same absolute time, leading to the apparent delay of the P-wave on HHZ. This alignment strategy was intended to clarify phase characteristics but may have caused confusion. We therefore have added a note in the figure caption on lines 245-247 to explicitly explain this: \"P- and S-wave arrivals, marked by white triangles, are derived from vertical and horizontal channels, respectively, and aligned by cross-correlation analysis between waveform pairs\".\n\nAlso, as suggested, we have removed the meaningless numbers on the vertical axis for Fig. 5c. The updated figure is attached for reference.\n\nFigure 5. (a) Seismic moment (M0) versus recurrence time (Tr) scaling for post-intrusion volcano-tectonic earthquake (VT) families. (b) Corresponding scaling for earthquakes in the 5 pre-intrusion long-period (LP) families. Linear regressions (black lines) and event distributions (circles color-coded by FI values) are shown for both panel a and b. Statistical descriptions are detailed in the text. (c) Example waveforms of anti-correlated LP multiplets recorded by station SN07. Red and black traces represent the pre-filtered (2-20 Hz) original and polarity-reversed waveforms of paired multiplets, respectively. P- and S-wave arrivals, marked by white triangles, are derived from vertical and horizontal channels, respectively, and aligned by cross-correlation analysis between waveform pairs.\n- 16 -\nIn the discussion of phase 2 beginning on line 269, the authors describe an absence of deformation and intense VTs. This does not seem well supported although it may be true. In Figure 4, that period is the only time with VTs in the caldera – though I agree there aren’t many – and the GNSS data show deformation across all 4 stations plotted. The problem is that the displacement data are normalized so one can’t tell what the scale of deformation is. I assume the reason for normalizing is to fit them all on one plot and show the consistency. In this case it would be helpful to have 4 scale bars showing the range of deformation represented by each of the stations.\n\nResponse: We have updated the figure caption to specify the normalization scale for deformation (lines 178-180): “Normalization of displacements is based on maximum displacement at each station, with values of ~14 cm for GV06 and ~10 cm for the other stations”. However, we hope to clarify that we do not share the view that the observation of a lack of deformation is insufficiently supported. The maximum amplitude of vertical displacement data from GNSS stations varies between 10~14 cm throughout Phase 2, and as stated in the text “… showed no significant changes above noise levels or signals that could be associated with magmatic or volcano-tectonic processes…” (lines 148-149). If this were a magmatic or volcano-tectonic signal, it would not have comparable amplitudes at stations within (GV06, 08, and 09) and certainly outside (GV10) of the caldera. To clarify our point that this is noise in the high-rate kinematic time series, we have added the following text to the caption for Figure 4: “The similarity in the amplitude of deformation across stations inside and outside of the caldera suggest noise in the high-rate solutions, and not magmatic or volcano-tectonic signals” (lines 180-182).\n\nFigure 4. Evolution of pre- and syn-eruptive seismicity on 26 June 2018. (a-d)\n- 17 -\nEpicentral distributions of seismicity during the 4 phases of activity at Sierra Negra caldera. Seismic events are color-coded by their FI. (e) Normalized vertical displacements of high-rate (30 s sampling) GNSS data recorded by continuous stations GV06, GV08, GV09, and GV10 during Phase 2 and the initial stage of Phase 3. Normalization of displacements is based on maximum displacement at each station, with values of ~14 cm for GV06 and ~10 cm for the other stations. The similarity in the amplitude of deformation across stations inside and outside of the caldera suggest noise in the high-rate solutions, and not magmatic or volcano-tectonic signals.\n\nIn Figure 2, could you mention why there is missing data for GV03?\n\nResponse: We have added additional text in the Figure 2 caption stating: “Missing data for GV03 are due to a data recording issue in the GNSS receiver” (lines 158-159).\n\nRegarding the caption for Figure 4, Phase 3 is mentioned as an ‘eruption episode’, but the eruption did not begin until 19:40 which marks the start of Phase 4\n\nResponse: We have revised the figure caption by deleting the words “eruption episodes” (line 178).\n\nIn Figure 5, the caption for part a should start with Seismic moment, not Moment magnitude\n\nResponse: We have corrected the typo in this figure caption (line 240).\n\nLine 212: I would argue that a R-squared of 0.07 is not a weak correlation, but is not a correlation at all\n\nResponse: We have revised the text on lines 232-233 (“LP families in Phase 2b lack statistically significant M -T scaling …”) and line 419 (“LP families do not exhibit\n0 r\nstatistically significant scaling …”) to ensure precise statistical clarification.\n\nLine 245: I suggest a better term than ‘near instantaneous’ for something close to 3 hours\n\nResponse: We refer the reviewer to our response to reviewer 1 regarding the terminology used.\n- 18 -\nReviewer #3 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\n\nResponse: We would like to extend our sincere thanks for your participation in the collaborative review process and your dedication to supporting the fair and comprehensive assessment of our work.\n- 19 -"}, {"round": 2, "reviewer_comments": {"reviewer_2": "As this is a re-review, I will keep my comments brief. I believe the paper is worthy of publication after the revisions. I commend the authors for the extensive work done to address the comments of the reviewers and feel satisfied with the responses to my comments."}, "author_response": ""}]

Decision letter