Earthquake-Derived Seismic Velocity Changes During the 2018 Caldera Collapse of Kīlauea Volcano

The 2018 Kīlauea caldera collapse produced extraordinary sequences of seismicity and deformation, with 62 episodic collapse events which significantly altered the landscape of the summit region. Despite decades of focused scientific studies at Kīlauea, detailed information about the internal structure of the volcano is limited. Recently developed techniques in seismic interferometry can be used to monitor the internal structure of an active volcano more directly by detecting subtle spatiotemporal changes in seismic wave velocity, but their utility relies on accurate interpretations of the underlying phenomena causing those velocity changes. Here, we retrospectively apply repeating-earthquake-based seismic interferometry to the 2018 Kīlauea eruption sequence. We find that seismic velocities changed over two distinct time scales: a sudden increase followed by a slower decrease in velocity in the hours following each collapse event, and a gradual, long-term decrease in velocity over several weeks that ceased approximately 1 month prior to the end of the eruption. Modeling suggests that short-term changes can be explained by magma reservoir pressurization which specifically closed vertical ring fractures. Long-term changes are related to subsidence of the caldera and likely include the influence of inelastic strain from the formation of new fractures. These observations provide new insights into the evolution of Kīlauea during its progressive collapse and will inform future interpretations for near-real-time monitoring at hazardous volcanoes around the world using similar techniques, especially where a dominant fracture orientation is present. shows when interpreting the relation between velocity change and deformation.

reservoir. By the time of the final collapse in early August, the caldera had subsided as much as 500 m and grown to ∼0.8 km 3 (Neal et al., 2019). Over 60,000 earthquakes were located by the USGS Hawaiian Volcano Observatory (HVO, 1956;Shiro, Zoeller, et al., 2021), and we analyzed the waveforms of these earthquakes with coda wave interferometry to study temporal variations in the structure of the shallow summit.
Coda wave interferometry has been developed in recent years as a tool to monitor changes in seismic velocity and has been successfully applied at volcanoes (Brenguier et al., 2008;Haney et al., 2014). The technique uses small differences in the waveforms of the seismic coda (i.e., late-arriving scattered waves following the direct body and surface waves) of repetitive sources to derive temporal changes in the sources or subsurface structure, especially seismic velocity. The repeating sources can be earthquakes, for example, multiple ruptures of the same fault patch separated in time (Battaglia et al., 2012;Hotovec-Ellis et al., 2014;Poupinet et al., 1984;Snieder et al., 2002); correlations of ambient noise (Brenguier et al., 2008;Duputel et al., 2009;Rivet et al., 2014;Sens-Schönfelder et al., 2014;Taira & Brenguier, 2016;Ueno et al., 2012); or active sources recorded on common receivers (Wegler & Lühr, 2001). At volcanoes, velocity changes may occur in response to the opening and closing of microcracks in the subsurface due to volcanic processes such as the intrusion of magma (Hotovec-Ellis et al., 2015;Ueno et al., 2012). However, seismic velocity can also change in response to a variety of non-volcanic processes such as hydrologic loading , earthquake shaking (Battaglia et al., 2012;Brenguier et al., 2014;Richter et al., 2014), stress changes (Taira & Brenguier, 2016), and temperature changes (Richter et al., 2014), which makes interpreting results a challenge. To reduce ambiguity, seismic velocity changes are often compared with geodetic (e.g., Interferometric Synthetic Aperture Radar (InSAR), Global Navigation Satellite System (GNSS; Rivet et al. [2014]; Sens-Schönfelder et al. [2014]) and weather data .
Several studies using ambient noise correlations or autocorrelations have been conducted at Kīlauea in the years prior to (Donaldson et al., 2017), months leading up to (Feng et al., 2020;Flinders, Caudron, et al., 2020;Olivier et al., 2019), and then during (Flinders, Caudron, et al., 2020;Wu et al., 2020) the 2018 eruption. Seismic velocity changes at the summit were found to often highly correlate with ground tilt on timescales of hours to weeks (Donaldson et al., 2017;Olivier et al., 2019;Wu et al., 2020), with loss of correlation interpreted to be due to deeper pressure sources (Donaldson et al., 2017), damage (i.e., crack) accumulation (Feng et al., 2020;Olivier et al., 2019), and co-seismic shaking and/or stress changes due to nearby earthquakes (Feng et al., 2020;Flinders, Caudron, et al., 2020;Wu et al., 2020). In all studies, pressure in the magmatic system and the subsurface response of existing cracks and/or the creation of new cracks are recurring interpretations.
Unlike the previous studies at Kīlauea that used correlations of ambient seismic noise, we instead utilized a large number of repeating earthquakes during the latter portion of the 2018 eruption to derive velocity change. In this paper, we compare our timeseries of seismic velocity changes with other relevant geophysical data timeseries (e.g., tilt, near-realtime GNSS) throughout the course of the eruption to better understand the processes to which velocity changes are sensitive.

Monitoring Network
All permanent seismic stations and geodetic instruments used in this study are maintained and operated by HVO (HVO, 1956;Shiro, Burgess, et al., 2021;Thelen, 2014). We utilized 13 3-component broadband sensors, two single-component short period sensors, and one 3-component short-period sensor. We also used data from 12 temporary nodal seismometers (Lin & Farrell, 2018;Wu et al., 2020). Figure 1 shows the location of instruments used in the study.

Velocity Change Measurements
We followed the method of Hotovec-Ellis et al. (2014) to derive a continuous function of seismic velocity from a catalog of repeating earthquakes. The repeating earthquake catalog used was from the near-realtime repeating earthquake detection algorithm REDPy (Hotovec-Ellis, 2021a; Hotovec-Ellis & Jeffries, 2016), which automatically sorts triggered earthquakes into groups (or families) of similar events. Individual events were triggered from continuous waveform data with ObsPy's short-term average/long-term average coincidence algorithm (Beyreuther et al., 2010; https://www.obspy.org/) with REDPy's default triggering and filtering settings. An earthquake was considered as part of a family of events when it cross-correlated with at least one other event within a family at or above a normalized cross-correlation coefficient (i.e., where autocorrelation is 1.0 on an individual station) of 0.7 on four of the 5 summit stations used (RIMD, UWE, OBL, HAT, and BYL; Figure S1 in Supporting Information S1). Cross-correlation was computed over a 10 s window starting 1 s prior to the triggering time, which includes some of the early coda. The REDPy catalog spans 27 April 2018, to 6 August 2018, and contains 73,348 earthquakes in 4,644 families ( Figure 1). This catalog was generated separately from the HVO Advanced National Seismic System Quake Monitoring System catalog and subsequently contains no inherent hypocenter information. Hypocentral locations (i.e., clusters shown in Figure 1) were assigned to clusters in the REDPy catalog by associating the entire catalog with the relocated catalog (Shelly & Thelen, 2019) by event time. The times of each event in both catalogs were compared, and the closest match was assigned so long as the difference in time was <1 s. We then took the median location of matched earthquakes within each family as the centroid. If a location was unable to be determined, we assigned it to the centroid of deformation (19.4073°N, −155.278°E). Map of instrument and earthquake cluster locations. Earthquake clusters are derived from high precision relative earthquake locations (Shelly & Thelen, 2019) as grouped by REDPy (Hotovec-Ellis, 2021a; Hotovec-Ellis & Jeffries, 2016), colored by depth relative to summit (average summit elevation ∼1.2 km a.s.l.), and sized relative to the number of earthquakes in that cluster. Solid purple lines represent pre-eruption topographic features; the dashed red line shows the extent of the 2018 collapse; black and gray lines show mapped fault and fracture features at the summit (Neal & Lockwood, 2003), dashed where approximate and dotted where obscured. Dark red dashed circle is the location and diameter of the median pre-collapse magma reservoir geometry estimated by Anderson et al. (2019) and used here in strain modeling. Stations noted in the main text are labeled (full station labels are available in Figure S1 in Supporting Information S1); colored stations correspond to data shown in Figure 2; stations outlined in thick black are used in the average in Figure 3. Note that UWE (seismometer), UWD (tiltmeter), and UWEV (GNSS) are co-located. Side and bottom panels show earthquake cluster depths in cross-section, with dashed line separating "shallow" and "deep" clusters considered in the text. Inset on bottom right shows location of map on the Island of Hawaiʻi.

of 17
We used these locations to determine the approximate P-and S-wave arrival times, assuming the REDPy trigger time was at the time of the P-wave arrival on the nearest station with a constant P-wave velocity of 5 km/s. We considered all pairs of events within a family that were separated by at least 30 min and cross-correlated them in a 5 s window beginning 0.25 s before the approximate P-wave arrival after filtering the waveforms between 0.5 and 15 Hz. Because the clustering process agglomerates many related sources, this step ensures the two events were similar enough in location and mechanism to consider prior to further processing. If the pair correlated above 0.7 in this window as recorded on the station and component chosen, we aligned the original earthquake-station waveforms and stretched their coda to determine the relative velocity change (dv/v, expressed as a percentage) between the time of the two earthquakes. We considered the coda to be a 14 s window beginning 0.25 s after the approximate S-wave arrival. We used a time-axis stretching with the objective function method of Kanu et al. (2013) to determine dv/v ( Figure S3 in Supporting Information S1). In a final check, we cross-correlated the reference and stretched codas and kept the measurement if the stretched coda increased the cross-correlation coefficient and it correlated above 0.65. We chose 0.7 as the initial cutoff between pairs of earthquakes as a compromise between the number of available measurements (and therefore the continuity of the resulting timeseries in the following step) and error in the determination of velocity changes due to differences in the earthquakes' locations (e.g., Niu et al., 2003). We assume that while differences in the locations of the earthquake pairs introduce uncertainty in the individual velocity measurements, there is no bias toward an increase or decrease of velocity, and therefore the agglomeration of thousands of measurements approaches the "true" value of velocity change. Figure S4 in Supporting Information S1 shows the effect of increasing the minimum initial cross-correlation value from 0.7 to 0.8 and 0.9, as well as the effect of choice of filtering passband. Note that the number of discarded measurements at high cross-correlation values strongly impacts the continuity of the final inversion resulting in an incomplete record.
Finally, we used a damped linear least-squares inversion to solve for the best continuous function of velocity change that matched our pairwise measurements with a temporal resolution of 30 min . The inversion uses first and second order Tikhonov regularization to find the smallest and smoothest solution, where we iterate to find the optimal tradeoff between misfit and smoothness to prevent overfitting the data (i.e., the corner of the L-curve). We rely on the sensitivity tests of Hotovec-Ellis et al. (2014) for confidence in which timescales are reliably resolved; Figure S5 in Supporting Information S1 shows the result of sampling a known function at the times of our measurements on station BYL and component HHE, adding Gaussian noise with standard deviation 1% dv/v, and inverting the noisy data to see how well the original function is recovered. The added noise is larger than the misfit of our data to inverted solutions, which has an average standard deviation of approximately 0.37% dv/v across all stations and components considered. Given the number of long-lived, large families of repeating earthquakes (e.g., Figure S2 in Supporting Information S1) and their significant temporal overlap, we have high confidence in our ability to resolve long-term features; however, in the short-term, the several hours following each collapse event has some reduced resolution due to lack of seismicity. This process was repeated for each station and component individually. Where available, we took the average velocity change across all components; the only 3-component station we were unable to use all components for was KKO, where the east component was malfunctioning for half of our study period. Figure S4 in Supporting Information S1 shows the results from individual components for the same station (BYL). Although there is likely additional information contained in the differences between components, we do not consider them here.

Results
The method we used relies on a high temporal density of repeating earthquakes, best enabling us to resolve temporal changes in seismic velocity (%dv/v) from late May to early August 2018 when discrete earthquake rates were highest (>500 located events per day). The velocity changes observed at stations (Figure 1) around the caldera all show a similar pattern: a linear decrease in velocity with the time that becomes approximately flat starting in late June superimposed on a sawtooth of rapid increases in velocity followed by smoother decreases. Figure 2 shows examples of the velocity change timeseries as observed on four stations around the caldera; results from all stations are available in Figures S6 and S7 in Supporting Information S1.
The changes we observe at an individual station are a spatial function of coda wave sampling, resulting in an "apparent" change that likely underestimates the true magnitude of spatially heterogeneous changes. For example, the magnitude of the velocity decrease from late May to early August exceeded 10% dv/v on several stations nearest the caldera, but at greater distances, both long-term and cyclic velocity changes were smaller (though they can be resolved to at least 12 km distance). We hypothesize that this is because significant changes were concentrated within the growing caldera, and coda waves for distant stations sampled relatively less changed volume than nearer stations. Furthermore, the large apparent velocity changes we observe seem as if they should invalidate the assumptions that go into coda wave interferometry itself (i.e., subtle changes in velocity structure or scattering). However, although we observe these changes over periods where the structure of the volcano was massively changing, the individual measurements contributing to these timeseries are much smaller and more closely spaced in time. We have no pairs that span the entire eruption, likely because the structure and stress conditions within the caldera changed too much. Even earthquake pairs that span collapses have velocity changes of less than 1% at the nearest stations. Our technique is unique in that we accumulate large changes by combining many small changes. Whether the accumulation of these changes linearly combines with the values we present is an issue we leave unresolved pending additional constraints.
We demonstrate the consistent features of these velocity changes across many stations around the caldera by showing average velocity change (Figure 3a, bottom timeseries) over six stations (UWE, HAT, BYL, KKO, RIMD, and WRM; outlined thick black and labeled in Figure 1). We then further compare this velocity change with complementary timeseries of seismicity (i.e., seismic moment release) and deformation (i.e., tilt radial to Halemaʻumaʻu crater, vertical GNSS position of a station on the crater floor, line length between two GNSS stations spanning the caldera), as well as the timing of several important transitions noted through the rest of the discussion. We divide our analysis into two distinct time scales: "long-term" changes occurring over days to weeks, and "short-term" changes occurring over many hours, including sawtooth cycles between and across individual collapse events.
The long-term velocity change observed to varying degrees at all stations (Figures S6 and S7 in Supporting Information S1) can be divided into two periods. During the first period (1 June to ∼26 June), velocity decreased at an approximately constant rate that was highest at stations near the caldera (dv/v between −0.2% and −0.4% Figure 2. Relative velocity change timeseries observed at four stations around the caldera (BYL, UWE, WRM, and AHUD; mean of all components shown; distance from station to centroid of deformation noted below name). Absolute vertical position is arbitrary. Triangle color corresponds to station color in Figure 1. Line color corresponds to number of measurements (i.e., half of a pair) within each 30-min time segment; line thickness is based on the number of pairs of measurements that span across the endpoints of each line segment such that thin or missing lines are less well-constrained by measurements before and after; gray area corresponds to ±1 standard deviation of measurements to the solution for each time segment. Blue dashed lines correspond to the time segments and naming conventions used for linear fitting in Figure 6. 6 of 17 per day), whereas during the second period (∼26 June to 5 August), the velocity did not change significantly. The exact timing of the transition between the two periods is obscured by the short-term changes, but it is clear that the approximate timing (within ±2 collapse events) is coincident with changes in several other observations. Changes in horizontal line length between GNSS sensors on opposite sides of the caldera (Figure 3a, blue timeseries) closely match the velocity change timeseries, with steady shortening that slowed around 21 June. The vertical position of CALS, located on the caldera floor inboard of the final collapse area, transitions from relatively little vertical subsidence to multiple meters of vertical subsidence each collapse; we highlight the collapse on 25 June where CALS dropped more than 1 m as an approximate transition. We note that although CALS and the crater floor subsided >100 m after this date, we observed almost no long-term seismic velocity change associated with this significant deformation. Aerial and satellite observations show that the peripheral fault scarp that had started to form to the north, east, and south of Halemaʻumaʻu crater completed its formation sometime between 18 June and 30 June ( Figure S8 in Supporting Information S1).
Several changes in the character of the seismicity also occurred in late June. For example, the 32nd collapse on 21 June had a distinctly different waveform with more high frequency energy compared to previous collapse events (Shiro, Burgess, et al., 2021), moment release showed a clear boundary between the 34th -36th collapses (24 and 26 June; Shiro, Burgess, et al., 2021), and the moment rate (i.e., a combination of earthquake rate and moment release for seismicity leading up to each collapse event) peaked prior to the 37th collapse on 27 June (Tepp et al., 2020). In the REDPy catalog ( Figure S2 in Supporting Information S1), many of the large repeating earthquake families (>200 members) formed around this time, representing the rapid creation of new earthquake sources. Afterward, relatively few new sources were created, and repeating seismicity largely reactivated sources that formed during this transition period. Indeed, when many seismic, geodetic, and geologic timeseries are compared in concert, the collapse on 24 June was a natural division between phases in the eruption (Flinders, Caudron, et al., 2020;Shiro, Burgess, et al., 2021;Tepp et al., 2020) and the onset of main-block subsidence following the complete formation of the peripheral fault. We elaborate on the relation between this fault and the velocity changes in more detail in the Discussion.
On shorter time scales, we observed that seismic velocity increased at the times of the collapse events and decreased between them. Collapse events are noted in Figure 3b as vertical bars, which correspond to the origin times of ∼M5 earthquakes, along with rapid subsidence of the GNSS sensor CALS on the caldera floor, extra-caldera inflationary deformation (as seen on the radial component of the UWD tiltmeter and GNSS line length), and abrupt reductions in seismicity (seismic moment). The observed velocity increase is somewhat surprising because large earthquakes are usually observed to cause an abrupt decrease in seismic velocity that slowly recovers (increases) to its previous value in the hours to days following the earthquake, sometimes even at great distance from the earthquake itself (Battaglia et al., 2012). These changes are often attributed to shallow, nonlinear responses to strong ground motions, perhaps involving fluids . While many of the largest earthquakes (M3.8-4.0) leading up to collapse events at Kīlauea are coincident with resolvable decreases in velocity (two examples noted in Figure 3b), the collapse events themselves have the opposite response even though they produced significantly stronger ground motions.
We cannot confidently resolve whether the increase in velocity occurred as rapidly as the geodetic offsets (i.e., over the course of <1 min) due to limited temporal resolution during the periods of low seismicity following collapses: the first 30 min following each collapse were almost completely aseismic, and few earthquakes occurred even 4-6 hr after the collapse. Furthermore, the inversions that transform pair-wise observations into continuous timeseries use a temporal smoothing parameter that prevents abrupt changes in velocity unless they are constrained by many measurements. However, a study with uniform temporal sampling (i.e., not reliant on earthquake rates) using ambient noise autocorrelations with a dense nodal array around the caldera  found a ∼4-hr delay between collapse and the maximum increase in velocity following each collapse event (Figure 3c). Wu et al. (2020) attributed the delay to the superposition of two competing effects: an instantaneous decrease in velocity due to shaking-induced structural damage followed by an exponential recovery due to healing and an instantaneous increase in velocity due to pressurization of the magmatic system followed by an exponential recovery as magma was withdrawn to supply the eruption.
We more directly compare our measurements to those of Wu et al. (2020) by estimating the average velocity change over a single collapse cycle. Instead of averaging our results after the inversion step, we performed a separate inversion with the time of each pair of individual velocity measurements referenced relative to the time of collapse to better mirror the method of Wu et al. (2020). For each pair, we set the time of the first earthquake relative to the time until the next collapse (i.e., negative time since collapse) and the second earthquake relative to time since the previous collapse. We also included copies of these measurements where both measurements were in time until the next collapse and time since the previous collapse to help more fully resolve the complete cycle. These measurements were then fed into the same inversion scheme as before, with a temporal spacing of 30 min. Despite the reduced resolution in the first few hours of the following collapse, our results are remarkably similar to Wu et al. (2020), but at a higher amplitude. Given the similarity between our observations, we assume the same delay and causative effects are present but less adequately resolved in our measurements.
In the next section, we discuss possible interpretations for the origin of these velocity changes and how they may relate to pressure changes in the subcaldera magma system.

Short-Term Velocity Changes
We first examine the mechanism by which seismic velocities abruptly increased and then decreased during and between caldera collapse events, respectively. Wu et al. (2020) attributed these changes to reservoir pressurization; here, we explore in greater detail the possible relations between changes in the volcanic system and velocity change. The opening and closing of microcracks caused by reservoir pressure changes is perhaps the most obvious mechanism for producing the relatively large and abrupt changes observed here, although numerous other mechanisms could also affect seismic velocities. To explain the temporal evolution of velocity changes during collapse cycles would require that cracks abruptly closed during (or immediately after) collapse events and then progressively opened prior to the next event. Here, we discuss potential sources of cyclic crack opening in detail.
First, we address the possibility that the velocity changes could be related to the stresses driving the occurrence of the prolific seismicity leading up to each collapse (the same earthquake sources used to derive the velocity timeseries themselves). A high correlation between the rate and amplitude of seismicity and seismic velocity change was noted during the start of the 2004 eruption of Mount St. Helens (Hotovec-Ellis et al., 2014), with the same sense of change as noted here at Kīlauea: velocity decreased as seismicity rates increased, then increased when seismicity abated. In rock mechanics experiments, increases in acoustic emissions and sample volumes and reductions in seismic velocities preceding sample failure are interpreted as the creation and opening of dilatant microfractures perpendicular to the applied compressive stress (Brace et al., 1966). There are parallels between these experiments and our system: fractured rock is put under stress, which created increasing amounts microearthquakes (i.e., acoustic emissions, illustrated in our case as an increase in pre-collapse seismic moment release) leading up to the next failure of the specimen (the collapse events). Velocity would then recover during (or just after) collapse events; the collapse event itself would have to fully destroy the dilatancy and reset the system. Although there is some experimental evidence (Holcomb, 1981) that rock samples under cycles of stress exhibit memory of velocity to stress state (i.e., seismic velocity is a recoverable function of stress), it is unclear whether this is an appropriate parallel. Additionally, while we conceive that stress cycles reloaded faults and could control near-fault dilatancy, the same relation does not apply to the pre-eruptive period discussed in Donaldson et al. (2017).
We hypothesize rather that observed inter-collapse seismicity and velocity changes were separate consequences of deformation due to pressure changes in the shallow magmatic system. We found a strong correlation between seismic velocity and short-term (∼hourly) tilt, which has been previously identified at Kīlauea on a longer (∼weekly) time scale (Donaldson et al., 2017), as well as with GNSS displacements. Ground tilt at Kīlauea's summit has long been interpreted in terms of pressure changes in the shallow magmatic system (e.g., Cervelli & Miklius, 2003;Fiske & Kinoshita, 1969), a relation that was clearly demonstrated in the correlation between ground tilt and lava lake surface height during the Puʻuʻōʻō eruption (Patrick et al., 2015). Caldera collapse events were manifested in tilt and high-rate GNSS data recorded outside of the collapse crater as outwards-and-upwards ("inflationary") offsets. These geodetic offsets are consistent to first order with pressurization of the Halemaʻumaʻu reservoir caused by the episodic, piston-like collapse of roof rock into the reservoir (e.g., Anderson et al., 2019;Kumagai et al., 2001;Segall et al., 2020), leading to down-rift surges in eruption rate . Collapse-induced ground deformations outside of the caldera were largely recovered between events.
The appearance of collapse offsets in extra-caldera geodetic data demonstrate that strain was imparted to extra-caldera rock. This strain is believed to have been largely imparted by collapse-induced pressurization of the magma reservoir, although ring fault slip likely played a role (Segall et al., 2020). Donaldson et al. (2017) suggested that changes in volumetric strain could explain changes in velocity, where for a homogeneous mix of crack orientations or spherical pore spaces the cracks will on average open when subjected to volumetric extension, thus decreasing the velocity (and vice versa). We modeled the strain due to a pressure increase in the subcaldera magma reservoir using a reservoir geometry previously inferred from geodetic (tilt, GNSS, and InSAR) and lava lake observations during the initial draining of the magmatic system before the onset of caldera collapse in 2018 by Anderson et al. (2019). In that study, the magma reservoir was modeled as a slightly prolate spheroidal body located east of Halemaʻumaʻu crater (Figure 1), centered at 19.4073°N, −155.278°E with centroid depth 1,940 m (−0.84 km a.s.l.), horizontal radius 914 m, and vertical radius 1,120 m. This model is the median of the suite of models which fit the deformation data. We used COMSOL Multiphysics software to calculate all components of strain due to a pressure increase of 3 MPa (the collapse-induced pressure change estimated by Segall et al. [2020]) in a spheroidal cavity with these dimensions embedded in a uniform elastic medium with a shear modulus of 3.0 Gpa, Poisson's ratio of 0.25, and assuming radial symmetry Segall et al., 2020). Figure 4 shows modeled strains caused by collapse-induced pressure increase in the magma reservoir (between collapse events, the sign of strain change would be opposite and occur over hours to days instead of seconds). We find that collapses induce relative volumetric extension in most of the regions around the reservoir, implying a decrease in seismic velocity-opposite of our observations. Volumetric strain transitions to compression at distance, but at lower amplitudes and too far away to explain our near-source observations. In contrast, Donaldson et al. (2017) predicted that a large portion of the subsurface directly to the sides of a pressurized reservoir at depth should be in relative compression when that reservoir is sufficiently shallow. Our methodology uses finite element modeling with an updated reservoir geometry, which more accurately reproduces strains at depth than the where blue corresponds to compressive strain (i.e., an expected increase in velocity due to crack closure) and red to extension. Dimension of the reservoir is in black; earthquake cluster locations relative to the centroid of the reservoir are yellow dots scaled in size to the number of members in that cluster. Triangles at top represent distances of nearby seismic stations used to measure velocity change. Earthquake clusters and station locations are projected onto this arbitrary radial plane; while a specific pair of an earthquake cluster and station may not have the distance and volume between them accurately portrayed in this view, the sense of change in the immediate vicinity of any individual cluster or station should be representative. Dashed lines correspond to depth division used in Figure 6 (6b)-(6d) Individual components of strain, modeled using the same reservoir as (a): (b) radial, (c) tangential, and (d) vertical strains with color scheme matching (a). analytic approximations used in Donaldson et al. (2017). Therefore, our short-term velocity change observations require a different explanation than changes in volumetric strain for a spheroidal reservoir.
While complexities in the reservoir geometry could conceivably change the sign of strain near the reservoir, we instead observe that while volumetric strain for a spheroidal reservoir does not predict our observations, the modeled radial strain (Figure 4b) does. In our paradigm where opening and closing cracks cause the velocity to decrease or increase, we should consider the orientation of those cracks. If there is a strong preference for crack orientation, crack opening or closing would be sensitive to strain parallel to their opening mode. In our radial cylindric coordinate system, radial strain opens vertical ring fractures concentric to the reservoir centroid. In Figure 1, we show that the majority of the mapped fault and fracture features at Kīlauea (Neal & Lockwood, 2003) are indeed concentric ring fractures. Although much of the macroscale fracture and fault features are in alignment with this conceptual model, it is unclear if the microscale structure follows the same trend, although studies have suggested that fracture network geometry may be to some extent scale-invariant (Bonnet et al., 2001).
If the small-scale structures important to changing the seismic velocity are of the same orientation as the larger ring faults, the pattern of radial strain predicts the sign of change for both our study and several others (Donaldson et al., 2017;Wu et al., 2020) without the need to invoke any further assumptions. Furthermore, our model predicts that the opposite relation between pressurization and velocity change could be observed at volcanoes without a dominant fracture pattern or with primarily radial fractures, even if they have similar reservoir geometries. For example, the fracture geometry at Piton de la Fournaise is more dominated by radial dike features (Carter et al., 2007) than Kīlauea, which may contribute to the opposite response of seismic velocity to reservoir pressurization (e.g., Rivet et al., 2014) in addition to the difference in reservoir depths between the two volcanoes.

Long-Term Velocity Changes
Next, we address whether this conceptual model applies to the long-term velocity changes. In Figure 3, we see that the timeseries of line length change between CRIM and UWEV is similar to velocity change on both the short-and long-term timescales. Relative line length change between two stations is a measure of average linear strain along this line. The line length between stations CRIM and UWEV, spanning the caldera, increases during collapses. On the other hand, the line between extra-caldera stations CRIM and AHUP is oriented approximately radial to the reservoir ( Figure S9 in Supporting Information S1) and decreases in length during collapses, consistent with the modeled radial compression outside the caldera. While the sense and timing of these changes qualitatively correspond well to our velocity change timeseries on both the hourly and monthly timescales, suggesting that the relation between radial strain and velocity change likely holds for both to some extent, the underlying phenomenon causing the strain need not be the same (i.e., pressure changes in the reservoir).
We evaluate whether the two timescales are different by examining their apparent locations. Although we do not formally solve for the spatial extent of the velocity changes (e.g., as in Obermann et al. [2013]), we can make some simple relative comparisons between the two. First, we analyzed how the timeseries of velocity change differ when using subsets of source earthquakes based on their hypocentral locations. If the long-and short-term changes were localized in the same place, we would expect to see a similar change in amplitude between two subsets of earthquake clusters regardless of which subset we chose. The most robust and easiest to interpret difference we found was between using earthquake clusters shallower and deeper than a given depth, which we arbitrarily chose at 1.5 km below the surface ( Figure 5); using only the shallower clusters in the inversion results in more long-term change than using only deeper clusters, while the short-term amplitudes remained approximately the same (see inset, Figure 5). We suggest that the short-term changes are more distributed in depth and the long-term changes are more concentrated near the surface. However, since the earthquake clusters are not uniformly distributed in either epicentral location or depth (e.g., the greatest concentration of shallowest clusters being on the southeast side of the caldera), we cannot uniquely determine the spatial extent of the long-term change using subsets of earthquake clusters alone. Additionally, there is significant uncertainty in the provenance of the coda and its spatial sensitivity, especially in depth. Got and Coutant (1997) suggest that coda waves at Kīlauea are multiply scattered in low velocity layers or dominated by converted surface waves, and would therefore have dominantly shallow sensitivity. Whether the depth of earthquakes is sensitive to the depth of changes or not, the difference in sensitivity to the long-term change is indicative of a difference in location.
Next, we examine how the amplitude of apparent velocity changes from station to station (here using all earthquake clusters). Generally, the stations closest to Halema'uma'u crater have larger amplitudes on both time scales than those that are farther away ( Figure S6 and S7 in Supporting Information S1), although the relation between amplitude and distance is more complicated in detail. Consider, for example, the difference in the timeseries for BYL and UWE (Figure 2) both at approximately 2 km distance from the centroid of deformation. BYL has more pronounced long-term changes in June compared to UWE, though they have similar short-term amplitudes. The two stations differ in elevation (i.e., UWE is on the caldera rim while BYL is on the caldera floor) and azimuth. The difference in azimuth suggests that the location of maximum velocity change during June may be centered further east than the centroid of deformation. In Figures 6a-6c, we show how the amplitude of change on both timescales decays as a function of distance from the centroid of deformation we used for our strain model to test whether this difference holds when considering all stations. For the long-term, we divided the timeseries into two time periods based on when the approximate break in slope for velocity change across many stations occurred (10 June-27 June and 27 June-1 August, noted as vertical dashed lines in Figures 2, S6, and S7 in Supporting Information S1), and plot the slopes from a linear regression of velocity change from those time periods ('first' and 'second'). Short-term amplitude was derived by removing the slope from the second time period, sorting the absolute values of the remaining timeseries, and plotting double the value corresponding to the 90th percentile as an approximation of the total change per average collapse event over that second time period.
Despite some scatter, we see that there is roughly an inverse relationship between amplitude and distance for at least the first long-term segment and short-term amplitudes. We resolve some similarities to Wu et al. (2020), who found that there was slightly more short-term amplitude apparent on the north side of the crater compared to stations on the south side of the same distance, though this relation is less clear for the long-term. The second long-term segment is mostly scattered around zero change, with some small but positive and negative slopes, and has two notable outliers at a near distance. For all three, we found the best fitting function of amplitude to Figure 5. Mean velocity change inversion results when only using earthquake clusters located above (red) or below (blue) 1.5 km depth. Line color corresponds to number of measurements (i.e., half of a pair) within each 30-min time segment, and line thickness corresponds to number of measurement pairs spanning before the start and after the end of each time segment. The two timeseries are overlaid to emphasize the difference in the early long-term rate of change and the similarity of short-term amplitudes. Inset shows similarity of short-term results in more detail.
inverse distance (i.e., a/R, where R is the distance from the station to the model centroid of deformation, and a is the best-fitting amplitude of scale), shown as a dashed line in Figures 6a-6c along with the r 2 coefficient of determination for each regression. This functional form was chosen because it is the simplest function that fits the first order decay of the amplitudes of velocity change with distance, even though it implies that velocity change should be infinite at zero distance from the source location as well as imposing symmetry is not reflected in the data. Despite this clear oversimplification, the first long-term slope and short-term amplitudes are relatively well-fit, with this functional form accounting for approximately 51% and 77% of the variability in the data respectively. As expected, the r 2 for the second long-term slope is quite low because there is not a clear relation between velocity change and distance. Although some additional scatter is likely due to the true location(s) of velocity change having finite volume instead of an unrealistic point source, the distance to each station is also a function of where the reference location is, which we have imposed explicitly.
Instead of assuming that the centroid of deformation is the correct point of reference, we did a grid search for all possible reference locations in the area to find the point with the highest possible r 2 for each of the time spans. Figures 6d-6f show the fits for the reference location with the highest r 2 , and Figures 6g-6i show how that value changes as a function of location. For the short-term, moving a few hundred meters north of the original reference location improves the explained variance to 91%, though many of the locations within the caldera have similarly high values. For both long-term segments, r 2 is improved by moving the reference location significantly to the north. Of course, this is largely controlled by the high, apparently outlier values at two stations: OBL and UWB, both on the caldera floor on the north side of Halema'uma'u crater. While the first long-term slope is not as well-resolved on those stations as many of the others, the second slope is well-resolved, so the unusually large second long-term amplitudes at those two stations are likely real. Whether the long-term changes are concentrated on the north side of the caldera due to a localized but pre-existing difference in the underlying system or that the process causing the changes itself is concentrated there remains unresolved.
We previously mentioned that decrease in the rate of long-term velocity change occurred at about the time that the peripheral fault formed, and a larger portion of the crater floor began to subside. First, let us consider that the first longterm slope is due to strain from loss of pressure in the shallow Halema'uma'u magma reservoir that, during the time of the second long-term slope and following the formation of the peripheral fault at the surface, was accommodated instead by sliding on that fault. If this were the case, we would expect the short-term amplitude and first long-term slope to be related, yet we have demonstrated that the short-and long-term changes appear to be concentrated in different locations. While it makes sense that the long-term changes are related to the subsidence of the piston given the timing of those changes, elastic strain alone may not be sufficient to explain our observations. We propose that part of this change could be due to inelastic strain from new fracture generation during the formation of the peripheral fault. These new fractures would be concentrated near the fault interface rather than distributed through the subsurface and could open more at shallower depths where confining stress is lower, potentially explaining the observed differences in location between the two timescales. Once the fault fully formed and was able to slip to the surface, fracture formation would have stopped or greatly slowed, and the velocity stabilized.
The mix of small but resolvable long-term increases and decreases in velocity observed on separate stations following the formation of the fault could be due to imperfections in the geometry of the piston as it subsided, allowing cracks to continue to open on one side of the piston (i.e., north side, with highest decreases in velocity) and close on the other (south and east sides). Another possible explanation for the long-term velocity decrease is a volumetric expansion of the collapsing overburden above the depleting reservoir, known as "bulking" (e.g., Poppe et al., 2015). In this paradigm, a decrease in velocity change would be concentrated in the volume above the reservoir and caused by an increase in void space due to inelastic processes during collapses. Physical experiments indicate that bulking occurs during progressive collapse (Poppe et al., 2015); however, our observations would suggest that bulking stopped when the peripheral fault forms. Figure 6. Analysis of velocity change amplitudes over different timescales with distance. First long-term amplitude corresponds to the slope of velocity change from 10 June to 27 June (a),(d); second long-term amplitude corresponds to the slope of velocity change from 27 June to 1 August (b),(e); and short-term amplitude is an estimate of the average total change across a single collapse event (c),(f). For each time scale, we show the best-fitting function of a/R (a, amplitude scaling factor; R, distance) as a dashed line, with r 2 goodness of that fit displayed in the corner. (a)-(c) have distance referenced to the centroid of the magma reservoir used for modeling, and (d)-(f) are referenced to the centroid with highest r 2 . (g)-(i) Maps of how r 2 varies as a function of location; centroid for (a)-(c) denoted with red star, and best fit centroid with a red cross. (j) Key for station colors used in (a)-(f) as well as throughout the rest of the paper. Stations UWB and OBL are denoted here as notable outliers.

Conclusions
The prolific, repeating seismicity of the 2018 Kīlauea eruption provides a unique opportunity to measure changes in seismic velocity with coda-wave interferometry. Our analysis yields a velocity change timeseries with high amplitudes (>8% total dv/v; >1% over several hours) and remarkable temporal detail during the last 2 months of partial summit collapses. We found that velocities correspond to deformation over two different time scales. Short-term changes (∼hours) were coincident with collapse events, where velocity increased following the collapse and decreased leading up to the next collapse. Long-term changes (∼weeks) can be roughly divided into two segments, where the first segment had a steady decrease in velocity, and the second segment had almost no velocity change. For both long and short time scales, we interpret velocity change as a consequence of the opening, closing, or creation of cracks in strain from different processes (Figure 7).  Figure 3C), annotated with the start and end times of each frame in red. Each cross-sectional frame shows the interpreted location and sign of change across each time frame as well as the major fractures involved (black), magma reservoir at depth (orange), relative internal pressure in that reservoir (red arrows), and relative motion (purple arrows). Note that although we presume the opposite sense of change above the reservoir here, we do not explicitly observe it. (b) Same as (a) but focusing on changes in long-term. Actual timeseries spanning June and July 2018 in gray at top, general trend overlaid in black, and times of start and end of frames in red. Focus of interpretation is on the formation of the eastern peripheral fault outlining the outer collapse margin. In both panels, presumed locations of increased velocities are denoted as blue, and reduced velocities as red.
Specifically, we attribute the increase in velocity following each collapse to vertical ring fracture closure due to the radial strain from pressurization of the magma reservoir (Figure 7a). We concur with Wu et al. (2020) that the ∼5-hr delay in peak velocity change is due to the superposition of velocity change caused by ground motion from the ∼M5 earthquake accompanying collapse and with pressurization of the magmatic system caused by the collapse itself. We propose that the presence of a dominant fracture fabric strongly affects the relation between pressure changes in the magma system and seismic velocity change. Further work is necessary for the modeling, experimental, and observational realms to test whether this interpretation holds. For example, we have not considered how the inclusion of a more complex geometric configuration (e.g., the inclusion of a piston, non-spheroidal reservoir geometries, or spatially variable material properties) would affect the strain, nor have we considered the possible role of fluids. However, our first-order modeling is consistent with our observations and several other previous studies at Kīlauea with few assumptions. Despite uncertainty in the strain field, the consideration of crack geometry may be important when interpreting velocity changes at volcanoes and other locations where a dominant fracture fabric may be present.
We attribute the long-term decrease in velocity (over days to weeks) to the growth of the peripheral caldera ring fault as it approached and ultimately reached the surface (Figure 7b). Unlike the recoverable, primarily elastic changes occurring over the span of a single collapse cycle, the longer-term changes likely incorporate some inelastic damage. Whether this is due to fracture formation from the growth of the peripheral fault or bulking of the overburden in the piston, the long-term changes are due, at least in part, to a different process than elastic strain from loss of pressure in the reservoir.
Coda wave interferometry and observations of velocity changes continue to hold promise for volcano monitoring. The technique requires few stations, can resolve signals several kilometers away from a volcano of interest, and are useful for identifying changes in the absence of other monitoring timeseries (e.g., direct observations of deformation). Although our work suggests velocity change may be useful as a proxy for strain, the geometry of the system must be carefully considered. This study highlights the continued challenge of unraveling a complex genesis of velocity changes in detail even when excellent complementary data are available.