Complex magmatic-tectonic interactions during the 2020 Makushin Volcano, Alaska, earthquake swarm

are often associated with the onset of unrest that can lead to eruption. However, determining whether seismicity reﬂects magmatic rather than tectonic stresses is often challenging, although critical for hazard assessments and risk management strategies. To investigate the triggering mechanisms of the recent Makushin seismicity, we integrate information from space-time patterns of the earthquake hypocenters with their fault-plane solutions. We relocate the swarm events using double-difference relocation techniques and a 3D velocity model and ﬁnd that the earthquakes, although they seem to follow two predominant orientations (NW-SE and SW-NE), do not show clear clustering into preferred alignments. Similarly, we do not observe pronounced migration in time and space. Fault-plane solutions (FPS) for all but one M2.5+ earthquakes have P-axis orientations consistent with subhorizontal NW-SE oriented regional maximum compression, whereas many of the lower-magnitude earthquakes have P-axes perpendicular to regional maximum compression. This provides evidence for the presence of a local stress ﬁeld likely induced by magma intrusion. Results from Coulomb stress modeling are also consistent with dike inﬂation modulated by stresses induced by the M4+ earthquakes. The seismic swarm is thus likely linked to a superposition of driving stresses from both magmatic and tectonic processes on pre-existing faults. The case of the 2020 Makushin swarm, with its unusual characteristics, challenges traditional swarm classiﬁcation schemes and suggests that a reconsideration of the deﬁnition of seismic swarms as having the maximum magnitude event in the middle of the swarm is warranted. © 2022 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Volcano-tectonic (VT) earthquakes are brittle-failure earthquakes directly or indirectly associated with volcanic activity (Mc-Nutt, 2005).They owe their name to their similarity to tectonic earthquakes, although the stresses that induce them are usually derived from volcanic processes such as magmatic intrusions (Mc-Nutt, 2005;Roman and Cashman, 2006;Chouet and Matoza, 2013), (Benoit and McNutt, 1996;Passarelli and Brodsky, 2012).However, although most eruptions are preceded by seismic swarms, VT swarms do not necessarily indicate an impending eruption (Moran et al., 2011;Newhall et al., 2017).
Determining whether a VT swarm reflects magmatic rather than tectonic stresses is challenging.This is especially true in areas lacking dedicated monitoring and/or in the absence of observations complementary to seismology, in particular crustal deformation and/or gas emissions.Fault-plane solutions (FPS) and patterns hypocenters locations have helped to distinguish between different swarm-triggering processes, as well as to clarify understanding of magma transport pathways (e.g., Kīlauea volcano [Karpin and Thurber, 1987;Roman et al., 2021]; Mount Hood [Jones and Malone, 2005]; Jailolo Volcano [Passarelli et al., 2018a]).
In this study we investigate the processes driving the intense earthquake swarm in June-December 2020 on Unalaska Island, Alaska, ∼12 km southeast of the summit of Makushin Volcano.The earthquake swarm produced two M > 4 events that were widely felt in the town of Unalaska, ∼15 km west of the epicenters.This is the strongest seismic activity at Makushin Volcano since its minor eruption in 1995 and since instrumental monitoring began in 1996.As of this writing, no eruptive activity or other surface changes have been observed at the volcano in webcam images, GPS or In-SAR.
Here we use data recorded by the Alaska Volcano Observatory (AVO) seismic network to relocate the 2020 swarm earthquakes.The improved precision of earthquake locations allows for wellconstrained FPS and Coulomb stress change modeling, which shed new light into the processes driving the swarm, in particular the relationships among stress state, faulting and magma transport.

Makushin Volcano
Makushin Volcano is located on Unalaska Island in the central Aleutians (Fig. 1).Makushin is a large ice-capped stratovolcano whose oldest lava flows date to the early Pleistocene (McConnell et al., 1997).It had multiple large eruptions in the early Holocene, including two caldera-forming events between ca.9-10 ka (Mc-Connell et al., 1997).Makushin's reported eruptions are exclusively focused at the central summit area and include >18 small-tomoderate eruptions since the late 1700's (Miller et al., 1998).The most recent eruption occurred in January 1995, producing a ∼2.5 km high steam and ash plume (McGimsey and Neal, 1996).In addition to aviation hazards, it also poses a potential local hazard to the ∼2,500 permanent residents of Unalaska Island, including ongoing commercial fishing activities and development of a new geothermal power plant.As a consequence, Makushin was identified by Ewert et al. (2018) as one of 18 very high threat volcanoes in the United States.

Seismic network and background seismic activity
AVO operates eight seismic stations (network code AV) around the summit area of Makushin, in addition to stations at neighboring Akutan and Okmok volcanoes.Currently station MREP is the only 1-component short-period station on the island.The other seven stations were already equipped or have been recently upgraded with three-component broadband seismometers (Power et al., 2020).The Alaska Earthquake Center operates an additional broadband station, UNV, located south of Unalaska (Fig. 1).Since 1996, background seismicity elongates in a NW-SE pattern from the volcano summit to ∼15 km E of the caldera.Earthquake depths shallow from ∼15 km to 0-7 km closer to the volcano summit (Fig. 1a,c).Previous seismic tomography studies (Syracuse et al., 2015;Lanza et al., 2020) have identified the seismogenic area as a possible brittle-ductile transition zone, with an aseismic magma body overlain by a brittle area that can support rock fracture (Fig. 1c, d).This interpretation is also consistent with geodetic models based on InSAR data by Lu et al. (2002) and Lu and Dzurisin (2014), who located an inflating source at 7 km BSL 5 km east of the caldera.Deep (∼15-37 km) long-period events (DLPs) associated with magma movement have also been observed below the summit area (Fig. 1c, d) and could represent migration of deeper magma towards a shallower crustal reservoir (Power et al., 2004).

June-December 2020 seismic activity
In June 2020, a seismic swarm began ∼12 km ESE of Makushin's summit caldera and ∼15 km west of the city of Unalaska (Fig. 1).The swarm began with a M4.2 earthquake at 21:16 UTC on June 15 at a depth of ∼8.6 km (https://earthquake .usgs.gov/ earthquakes /eventpage /us6000acjk /executive) that was shortly followed by a M4.1 at 00:34 UTC on June 16 at a similar depth of ∼8.7 km.Seismic activity lasted for several months until December 2020; bursts of seismicity were also observed in October 2021 and January 2022.However, most of the seismic unrest was concentrated in the first weeks of the swarm, so we restrict our analysis to June-December 2020.The earthquake swarm was felt locally by the population of Unalaska, including eight M > 3 events in addition to the two events greater than magnitude 4 (Fig. 2a).AVO located 1,509 VT events from June 2020-December 2020 of which 1,354 are within 3.5 km from the swarm centroid; in addition to the VTs, a few long-period (LP) seismic events were also detected between July 1 and September 8, 2020.LPs were emergent and characterized by energy concentrated at around 2 Hz, with the highest amplitude on station MREP (∼6 km south of the swarm volume), followed by station MGOD (∼10 km SE of the swarm volume).Examples of waveforms for VT and LP events are shown in Fig. S1.Initial catalog depths placed the VT swarm earthquakes between 7 and 10 km BSL at the termination of the lineation of background seismicity that extends towards the SE from the summit (Fig. 1b).
The 2020 swarm is not easily characterized as a simple swarm or mainshock/aftershock sequence.On the one hand, the largest event (M4.2) occurred as the first event in the swarm, consistent with the definition of a mainshock-aftershock sequence.This behavior resembles the 2006 Mount Rainier swarm, which started with a M4.5 located 10-12 km east of Rainier's summit and gradually died off over 30 days (Hartog et al., 2008).Moreover, as for the Mount Rainier swarm, the number of aftershocks as a function of time follows well an Omori-type decay (Fig. 2b) further suggesting a typical mainshock-aftershock behavior.On the other hand, the absence of a correlation between the total duration of the swarm and its seismic moment is a first-order indication that the 2020 seismic activity at Makushin deviates from a normal mainshock-aftershock sequence where instead either a linear or cubic moment-duration scaling is expected (Passarelli et al., 2018b(Passarelli et al., , 2021)).In the following sections, we investigate this complex seismic behavior through analyses of temporal evolution of b-values, earthquake relocations, and FPS.

Earthquake magnitudes and b-value estimation
Assessing the statistical and quantitative properties of earthquake swarms/sequences can aid in distinguishing their driving processes.One of the best-known descriptions of earthquake statistics is the Gutenberg-Richter distribution, which describes the relative numbers of small and large earthquakes as log(N) = a −bM  (Gutenberg and Richter, 1942), where N is the cumulative number of earthquakes larger than a given magnitude M, a describes the rate of earthquakes at or above M = 0, and b describes the relative proportions of small and large earthquakes.This latter 'bvalue' for a given fault distribution is thought to be controlled by the effective stress on the faults (Wiemer and McNutt, 1997;Bachmann et al., 2012;El-Isa and Eaton, 2014).For example, if the effective stress on a fault is decreased due to a change in crustal stress conditions or the presence of fluids in a volcanic or hydrothermal system, then the b-value would be higher (reflecting a greater proportion of low-magnitude earthquakes) than for typical tectonically-driven seismicity.High b-values can also be due to greater heterogeneity in crustal structure (i.e., many small faults as opposed to a few large faults) and/or increased thermal gradients or release/exsolution of magmatic gases (e.g., Wyss et al., 1997); this condition is often found in crustal volumes located above and around a magma chamber, as a result of applied pressure from expansion (Bridges and Gao, 2006).b-values can be therefore very useful at providing an additional constraint on whether seismicity is of volcanic or tectonic origin.
The choice of magnitude scale has significant implications for the b-value estimated from populations of low-magnitude events (Shelly et al., 2016;Hudson et al., 2022).In particular, for a swarm composed of very-low magnitude events, the estimated b-value may be up to 50% higher when using Mw (Moment magnitude) compared to M L (Local magnitude).In fact, while Mw provides an estimate of the actual seismic moment or energy release of an earthquake, M L is expected to scale directly with the logarithm of the seismic moment; therefore, a b-value derived from M L is perhaps more similar to the β parameter relating frequency and moment (M 0 ) for certain magnitude ranges (Deichmann, 2017), log(N) = a − β log(M 0 ), where b = 1.5β (Wyss, 1973;Kagan, 1991).
In this study, we calculate Mw using the spectral method described in Stork et al. (2014) and the open-source code SeisSr-cMoment by Hudson (2020).For each event, we first correct for instrument response and convert the velocity seismograms to dis- placement.This step is followed by calculating the spectrum of the displacement signal using a multi-taper spectrum method (Prieto et al., 2009).By fitting a Brune source model to the resulting spectrum, we derive the long-period displacement spectral amplitude, which is then used, together with the earthquake-receiver hypocentral distance and a radiation pattern correction term, to calculate the seismic moment release, M 0 .Mw can then be calculated from Hanks and Kanamori (1979) as: For the b-value calculation we use the B-Value Stability (BVS) criterion method of Roberts et al. (2015).The BVS method uses the maximum likelihood method to estimate the b-values, then it evaluates the estimated b-value as a function of the cut-off magnitude.Finally, to estimate the temporal b-value variation and corresponding uncertainties, we follow the method described in Roberts et al. (2016), in which the chronological earthquake catalog is sampled at different windows of random lengths and the b-value and associated Probability Density Function (PDF) are calculated for each window and then stacked to produce the overall b-value and associated PDF through time.Because changes to location result in changes to hypocentral distance and hence computed magnitudes, we perform the statistical analysis described above with the relocated catalog (see Section 4).
Fig. 3 shows the full evolution of swarm magnitudes for 888 events in our relocated catalog for which we were able to obtain a reliable Mw.The magnitude of completeness (Mc) is equal to 1.56 and is one magnitude unit larger for Mw compared to M L , for which we obtained an Mc of 0.36.The Mw b-value of 1.56 is significantly greater than 1, and also in stark contrast to the M L b-value of 0.74 (Fig. 3).These differences are not unexpected for very small events as the Mw and M L moment scales deviate from ∼2/3 ratio for Mw < 3 (Deichmann, 2017;Hudson et al., 2022, see also Fig. S2).
Estimated b-values from 5000 random windows of sizes 50 to 500 (for events with Mw ≥ Mc) show a consistent behavior through time with values always above 1 for the duration of the swarm (Fig. 3b).The swarm is located in an area close to a region of low Qp (50-150) and low Vp/Vs ratio (Lanza et al., 2020;Fig. S3) and in proximity to the tomographic magma centroids (Fig. 4, Syracuse et al., 2015;Lanza et al., 2020).On the one hand, theoretical and empirical studies (O'Connell and Budiansky, 1977;Peacock et al., 1994) have shown that increased fracturing and crack density lead to an increase in attenuation (low Q).In addition, the presence of gas would decrease the bulk modulus of the rock, leading to a decrease in the P -wave velocity.However, this would not significantly affect the S-wave velocity, resulting in a low Vp/Vs ratio (Nakajima and Hasegawa, 2003).As a consequence, the region of the swarm occurrence could be interpreted as an ensemble of seismically-active magma/fluid-filled sills, dikes, and/or fractures, something that is also reflected in the complexity of the earthquake relocation patterns (Section 4).As an alternative to exsolution of gas coming from a magma body stored nearby or stresses created by a deeper magmatic fluid source on preexisting fracture zones, the high b-values could be linked to the presence of hydrothermal fluids circulating in or above a supercritical fluid reservoir (e.g., Kirishima Volcano [Chiba and Shimizu, 2018]; Uturuncu Volcano [Hudson et al., 2022]).The known geothermal field of Makushin is not far away, located on the east and southeast flanks of Makushin Volcano, at the head of Makushin Valley (Fig. 1), although this would require the hydrothermal reservoir to extend quite deep in the crust and it may not be physically feasible.

Relocation procedure
The initial AVO earthquake catalog used for relocation purposes contains 1,414 earthquakes between June 14, 2020 and December 7, 2020 and within 25 km of the summit of Makushin that was prepared using the methodology described by Power et al. (2019).For our relocation analysis, we use AVO's manually picked P -and S-wave arrival times.These phase readings amount to 8,980 P picks, and 7,373 S picks, which were weighted according to a standard scheme (from 0 to 4) in which the highest weight (0) was assigned to picks with uncertainties of <0.02 s (i.e., the highestquality picks) and the lowest weight (4) was given to picks with uncertainties of >0.5 s (we exclude picks with weights of 4 from further analyses).In addition, the minimum location quality criteria (number of observations -nobs-and azimuthal gap -gap-) are set at P nobs ≥ 4 and gap < 180 • , decreasing the number of useable catalog earthquakes to 1,210 with a total of 8,135 P picks and 6,643 S picks.Although we identified additional uncataloged smaller earthquakes using an automated detection algorithm (REST, Lanza et al., 2019), in this study we opt to use only the higher quality manually-picked events to decrease the occurrence of possible misidentification of P and/or S phases.
We performed initial relocations using the seismic tomography code simul2014, extended from simulps12/simul2000 (Thurber and Eberhart-Phillips, 1999;Eberhart-Phillips and Reyners, 2012).Simul2014 uses a 3D velocity model (in this case from Lanza et al., 2020) and a linearized earthquake location method (in our relocations, we did not allow the velocity model to vary).For this initial linearized earthquake relocation, the weighted root-meansquare misfit is reduced from 0.080 s to 0.065 s, a reduction of ∼17.5% with respect to the initial AVO catalog locations.The locations obtained from simul2014 were used as initial input for hypoDD 3D (v2.1b), which uses the double-difference algorithm of Waldhauser and Ellsworth (2000) with a 3D velocity model and cross-correlation differential times to produce final high-precision hypocenter locations.The dynamically weighted double-difference inversion of the simul2014-relocated dataset uses 5,541,317 pickderived differential times and 354,002 cross-correlation-derived differential times calculated between pairs of events with a maximum hypocenter separation of 5 km and a normalized crosscorrelation coefficient >0.8.We compute cross-correlations in the time-domain using filtered (1-20 Hz) seismogram pairs beginning 0.1 s before the pick time ( P or S arrival) and ending 0.25 s after the pick time.The maximum allowed hypocenter separation of 5 km discards ∼5% of the events, further reducing the number of relocated events from 1,210 to 1,154.We exploit here the ability of hypoDD to combine differential times derived from phase-picked data and waveform cross-correlation and perform a joint inversion.
A crucial aspect of the relocation process is the weighting and reweighting of the different kinds of data (Waldhauser and Ellsworth, 2000).Here, we used six data sets with five iterations and chose an appropriate reweighting of the differential times for each set.For the first 10 iterations, we downweighted the crosscorrelation data by a factor of 100 to allow locations using only the catalog differential time data for event pairs with large separation distances (4-6 km) to define the overall swarm structure.For the following 10 iterations we downweighted the catalog-based data also by a factor of 100 and allowed the cross-correlation differential times to locate the event pairs having shorter separation distances (∼2-4 km) to refine the event locations.The last 10 iterations maintained the same weighting scheme, but we used events pairs with an even smaller separation distance of 0.5-2 km.All of the calculations were performed using the conjugate gradients method after appropriate damping.Because the uncertainties in the final locations estimated by hypoDD using the LSQR method are grossly underestimated and poorly constrained (Waldhauser and Ellsworth, 2000), we assess them independently using a bootstrap resampling method (Mesimeri et al., 2018) in which relocation is repeated for 200 samples.The resulting average median 95% confidence error ellipse semi-major axis lengths are ∼24 m, 11 m, and 14 m in the model-aligned x, y, and z directions, respectively (Fig. S4).Given the uncertainties in the seismic velocity model, we acknowledge that these are still probably underestimated; however, the relative locations are well-constrained, as they depend less on the velocity model (e.g., Waldhauser and Ellsworth, 2000;Zhang and Thurber, 2003).

High-precision earthquake relocations
Precise relative relocations of 1,154 earthquakes reveal a complex pattern.Although there are no clear alignments, the events seem to follow two predominant orientations, NW-SE and SW-NE (see cross-section profiles AA' and CC' in Fig. 4; there are no mapped surface faults for this area).The temporal evolution, with ∼48% of the events occurring in the first 2 days of the swarm, indicates no obvious earthquake migration (Fig. 5).In August, a small cluster oriented NW-SE activated in the east portion of the main swarm seismicity.A weak shallowing in depths is then observed for the December events, which locate 1-2 km shallower and east with respect to the main cluster (Figs.5c, 6).In general, we observe a consistent shift in depths with respect to the initial catalog locations, with the relative relocated events ∼1-2 km shallower than catalog depths (Fig. 6).The lack of spatio-temporal evolution is not uncommon, and it has been observed in well-recorded preeruption swarms (e.g., Roman andCashman, 2006, 2018).

Fault Plane Solutions (FPS)
For all VT earthquakes between June 15, 2020 and December 31, 2020 with an AVO catalog location (USGS, 2017) at Makushin (latitudes 57.78 • -54.01 • N, longitudes 166.9 • -166.62 • W), a catalog magnitude greater than 1.0, and a catalog azimuthal gap <180 • (281 earthquakes), we repicked P -wave first-motion polarities (when a minimum of six could be identified clearly) and combined our polarity picks with the relative relocation for the event (from Section 4.2) to calculate a double-couple FPS.FPS were calculated using FPFIT (Reasenberg and Oppenheimer, 1985) and the AVO 1D velocity model for Makushin (Dixon et al., 2013), and evaluated for quality and rejected if they did not meet all of the following criteria: unique or similar multiple solutions (e.g., two strike-slip solutions with similar orientations), misfit of no more than 15% of the P -wave polarities, uncertainty in strike, dip, and rake of ≤25 • , and pressure ( P ) and tension (T ) axis 95% confidence regions that cover ≤25% of the focal sphere.FPS that met these quality criteria have an average of 11 first-motion polarities, a maximum of 23 first-motion polarities, and a minimum of 7 first-motion polarities.While the number of first-motion polarities for individual FPS does generally correlate with magnitude, most events are M L 1-2 (as shown in Fig. 8): M L 1-1.5 events have an average of 10 first-motion polarities, and M L 1.5-2 events have an average of 11 first-motion polarities.Example FPS, along with waveforms indicating first-motion polarities, are shown in Fig. S5.FPS for 11 'background' events occurring in the swarm volume between March 2015-October 2018 show oblique strike-slip faulting (and a minor component of oblique normal faulting) consistent with NW-SE regional maximum compression (Fig. S6).
In total, we were able to calculate 109 double-couple FPS for swarm earthquakes of M L 0.5-4.2 between June and December 2020.We classify FPS as dominantly strike-slip (P-and T-dip <45 • ), thrust (P-axis dip <45 • and T-axis dip >45 • ), or normal (P-axis dip >45 • ) and further classify these events as consistent (∼ σ 1R ) or inconsistent with the regional stress field ( σ 1R ): Strike-slip and thrust FPS with P-axis azimuths within 45 • of 327 • (assumed σ 1R ), and normal FPS with T-axis azimuths within 45 • of 57 • (assumed σ 3 R ), are classed as consistent with the regional stress field; all others are classed as inconsistent.
Approximately 43% (47/109) FPS indicate oblique strike-slip, thrust, or normal mechanisms inconsistent with the orientations of regional maximum and minimum compression (Fig. 7).The remaining 62 FPS indicate a mix of oblique faulting mechanisms consistent with the regional stress field.All but one of the FPS consistent with regional stresses have M L 2.5 or greater (Fig. 8).The inconsistent FPS are not distinct in space or time from consistent FPS, but concentrate in the east and west portions of the swarm volume (Fig. 7, 8).The inconsistent FPS provide evidence for the presence of a local stress field induced by dike inflation/magma intrusion during the 2020 swarm (e.g., Roman and Cashman, 2006;Lehto et al., 2010;Roman and Gardine, 2013;Roman et al., 2021).

Coulomb stress change modeling
Using Coulomb 3.3 (Lin and Stein, 2004;Toda et al., 2005), we forward-model two possible end-member scenarios: 1) Strike-slip earthquakes inconsistent with regional stresses are the result of stress transfer (i.e., aftershocks) from the two M4+ earthquakes, 2) strike-slip earthquakes inconsistent with regional stresses are induced entirely by stresses from an inflating dike.In both scenarios, receiver faults are represented by a left-lateral pure strike-slip fault with a strike of 105 • (the average fault strike for inconsistent strike-slip FPS).For scenario (1), we use the locations and FPS for the two M4+ events obtained in this study as sources of Coulomb stress change.In scenario (2), the position of the inflating dike (the stress source) is obtained through optimization of Coulomb stress changes resulting from a model dike following the approach of Roman et al. (2011).Specifically, in scenario (2), we searched for a model dike that maximized the number of earthquakes with inconsistent FPS that experienced a Coulomb stress increase ( σ f ) based on their location, fault orientation, and sense of slip.To limit the number of free model parameters, we assumed an isotropic regional stress field, a dike orientation of 327 • (parallel to the NUVEL-1A plate convergence vector, Fig. 7 inset) and a 1-km-tall dike centered at a depth of 8.5 km BSL.We also assume 1 m of instantaneous opening based on average dike widths in Thiele et al. (2020).Tested dike lengths ranged from 0.1-1 km (in 100 m increments).Tested dike locations were within a 3 km (latitude) by 5 km (longitude) grid centered on 53.8625 • N and 166.77• W with 200 meter increments.All calculations were made in a homogeneous linearly elastic half-space with typical elastic parameters including a shear modulus of 5.0 × 10 4 MPa, Poisson's ratio of 0.25, and an effective coefficient of internal friction of 0.4.
The results of our Coulomb stress modeling indicate that swarm FPS orientations and sense of slip are consistent with emplacement and inflation of an ∼0.6-km-long model dike located just south of the swarm centroid during the second half of 2020, as shown in Fig. 9.We find that the majority of earthquakes with inconsistent strike-slip FPS (P-axes NE-SW) that are located west of 166.77 • W are equally-well modeled by Coulomb stress changes resulting from the two M4+ events or from dike inflation (i.e., they are located in regions of positive Coulomb stress change).However, earthquakes with inconsistent strike-slip FPS (P-axes NE-SW) that are located east of 166.77 • W are only consistent with Coulomb stress changes from dike inflation, providing Fig. 7. Map (upper panel) and cross-section (lower panel) of FPS, coded to indicate whether they are normal, strike-slip, or thrust, and consistent or inconsistent with regional stress, based on P-and T-axis orientation (see text for details).Inconsistent FPS are highlighted with yellow P-axis quadrants.Arrow in inset shows the NUVEL-1A (DeMets et al., 1994.)relative motion between Pacific (moving) and North American (fixed) plates at the latitude and longitude of Makushin.(For interpretation of the colors in the figure, the reader is referred to the web version of this article.)strong support for dike inflation as the dominant driving mechanism of these earthquakes.Misfits of individual events may result from variability in fault orientation and rake, spatiotemporal complexity in the inflating dike, absolute location errors, and the absence of regional stresses in the model scenarios.While for simplicity we restrict our Coulomb stress modeling to inconsistent strike-slip earthquakes, it is reasonable to posit that inconsistent thrust and normal earthquakes also result from a combination of intrusion-induced stresses and elevated pore pressures resulting from magma intrusion, as well as additional unmodeled spatiotemporal intrusion complexity.

Tectonic or magmatic?
Because the June 2020 seismic sequence initially appeared to be a mainshock-aftershock sequence (MS-AS), it led to significant debate as to its cause (tectonic vs magmatic).Here, we present multiple lines of evidence from our analysis implicating magma intrusion as the predominant driving process for the 2020 swarm, prompting us to re-examine the fundamental distinctions between tectonic and magmatic earthquake sequences.
The seismic sequence's temporal evolution coupled with its cumulative seismic moment and seismogenic volume can help discriminate between different types of earthquake sequences (Vidale et al., 2006).For the Makushin seismic activity of 2020, the temporal duration of ∼6 months associated with the relatively small mainshock is indicative of a characteristic duration that is uncoupled from the seismicity.Its cumulative seismic moment-duration relationship is incompatible with "regular" (MS-AS) tectonic earthquake sequences as seen in Fig. 10a (Holtkamp and Brudzinski, 2011), where the seismicity burst plots in the swarm-type domain.Similar observations in which swarm sequences show a deviation from the scaling behavior of ordinary earthquakes are documented in Passarelli et al. (2018b) and Passarelli et al. (2021) and also shown in Fig. 10a.We also compared the Makushin 2020 swarm to the Vidale and Shearer (2006) and Vidale et al. (2006) analyses of Southern California and Japan "bursts" (Fig. 10c-d).To ensure a fair comparison, we have followed the same selection criteria as outlined in Vidale and Shearer (2006).Thus, we consider only events within a radius (defined as the average distance to the swarm centroid) of 2 km that occurred within 28 days from the mainshock.This limits the analysis (and duration of the swarm) to July 13, 2020, and includes 860 This window was shortened to 14 days (and therefore 736 events) when comparing the Makushin swarm with the Japanese swarm sequences as stated in Vidale et al. (2006).The other two requirements in Vidale and Shearer (2006) are also valid here, as there were only few events in the prior 28 days within the same 2 km radius, and there were no more than 20% more events (∼1%) between 2 and 4 km from the initiating event in the same 28 days afterward.We note that although for the Japanese sequences the radius was increased to 4 km, we maintained the same 2 km radius in our analysis.Fig. 10c-d shows the Makushin swarm to be again clearly in the domain of swarm-like sequences.We note here that the Makushin swarm is an outlier or at the edge of each of the four plots shown in Fig. 10.This makes it a good case for comparison in future studies.
Another simple, yet robust statistical indicator to distinguish a sequence type is the skewness (S) of the temporal distribution of the moment release in a cluster (Roland and McGuire, 2009;Passarelli et al., 2018b).It has been demonstrated that MS-AS sequences have positive skewness (S >> 8), whereas swarms have a small or negative skewness (−4 < S < 4), as the moment release is often delayed from the beginning of the sequence.When coupling the analysis of S with the Kurtosis (K) of moment release history, it allows for a parabolic representation where MS-AS and swarm sequences can be clearly separated (Mesimeri et al., 2019).For the Makushin swarm, we calculated S and K consider-ing the entire duration of the swarm (∼6 months) on (1) the 888 events using Mw calculated as in section 3; (2) the relocated 1,154 events using M L ; and (3) the relocated 1,154 events using Mw derived from a quadratic fitting of the obtained Mw-M L relationship (see Fig. S2).For all three cases we obtained positive values for both S and K (4 < S < 12; 22 < K < 231; Fig. S7).These results put the Makushin swarm right at the transition between MS-AS type and swarm sequences, with a tendency to rather follow a MS-AS behavior, showing that the sequence could have initiated as a tectonic-like sequence (with an abrupt begin and a slow decay) before transitioning to a complex swarm-type event.Evidence for a tectonic contribution to the swarm in its early stage is also confirmed by the swarm decay pattern as it follows the modified Omori-law (Utsu et al., 1995), which usually well-describes the aftershock decay in a MS-AS sequence (Fig. 2b).
We next address whether the swarm is driven by magma intrusion or fluid circulation.Typically, swarms fitting a diffusivelike pattern in space are interpreted as related to fluid pressure propagation (e.g., Hill and Prejean, 2005;Shelly et al., 2013b).In this swarm, we find a lack of pronounced migration of earthquakes outward from the initial source, except for a weak upward migration of the earthquake depths in early December (Fig. 6).In addition, in Fig. 10d we also show how the seismogenic volume of the swarm does not scale with magnitude, as was observed by Vidale et al. (2006) for fluid-driven Japanese swarms.Similarly, the elevated b-values (1.56) seem to favor the hypothesis of the presence of a crust that is being more extensively cracked than elsewhere as a result of stress triggering (i.e., mag-  7).Upper-left inset illustrates the modeled receiver faults, which are idealized left-lateral strike-slip events with average fault orientation of inconsistent strike-slip events as shown in Fig. 7. Small black symbols show epicenters of all catalog earthquakes in the study period.(For interpretation of the colors in the figure, the reader is referred to the web version of this article.)matic processes such as intrusions, dike formation) rather than fluid triggering and fluid diffusion.Often elevated b-values are associated with (1) high thermal gradients, (2) high pore pressure found in the vicinity of a magma reservoir, or (3) vesiculation and fragmentation of ascending magma that can create stress-induced cracks (Wyss et al., 2001).Further, the FPS inconsistent with regional stresses in zones of positive Coulomb stress change induced by dike inflation (Fig. 9) provides additional evidence for the presence of a local stress field induced by a magmatic intrusion (e.g., Roman and Cashman, 2006;Lehto et al., 2010;Roman and Gardine, 2013).Local ∼90 • stress field rotation can be explained by failure within the walls of inflating dikes (Roman, 2005) and therefore can be indicative of magmatic processes.The largest-magnitude earthquakes all have solutions consistent with the regional compressive stress field (Fig. 8).These FPS may be attributed to regional stresses, and could explain the observed initial MS-AS characteristics.Thus, the observed heterogeneity of the FPS as well as the ambiguity of the sequence statistics likely point to a superposition and rapid modulation of background tectonic stresses and transient stress perturbations due to a magmatic intrusion.
During many intrusions, additional corroborating observations such as increases in degassing and/or ground deformation are observed, whether before or simultaneously with the occurrence of the swarm seismicity (e.g., Sigmundsson et al., 2018).No increased SO 2 or deformation during the Makushin sequence was observed.Neither GPS nor InSAR recorded ground deformation.Computed InSAR interferograms indicate that there was no significant deformation (larger than a couple of cm) during 2019-2020 as well as for a smaller time frame focused on the largest earthquake recorded on 2020.However, many intrusions, especially small-volume and/or deep (>5 km BSL) intrusions, are not accompanied by observable degassing or deformation (e.g., Shishaldin Volcano [Moran et al., 2006]).Thus, the absence of multi-parameter evidence does not rule out deep and/or smallvolume intrusion.On the other hand, we acknowledge that the lack of complementary geophysical datasets, as well as detailed knowledge of mapped faults and their geometry and simplified assumptions of the stress modeling, pose limitations to our interpretation.

Spatial relationship of earthquake activity to magmatic intrusion
VT earthquakes result from a combination of stresses induced by volcanic fluid/magma transfer and regional stresses and the presence of suitably-oriented fault planes, thus the location of intruded magma may be distinct from the resulting VT hypocenters (e.g., Roman, 2005).Our results from the Coulomb modeling  Holtkamp and Brudzinski (2011) and Ide et al. (2007).In addition to the Makushin swarm (red circle) we have plotted the Icelandic tectonic earthquakes swarm (TES) from Passarelli et al. (2018b) and the swarmgenic slow-slip events (SSE) as shown in Passarelli et al. (2021).Duration is time between the first and last events in the swarm or sequence, and magnitude is the total seismic moment of earthquakes that comprise the swarm or sequence.(b-d) Characteristics of the 2020 Makushin seismicity compared against Japanese swarms and a of Southern California sequences as measured using parameters from Vidale and Shearer (2006) and Vidale et al. (2006).Following Vidale et al. (2006) terminology, "Average" indicates those sequences that do not fit aftershock-like or swarm-like sequences and have in between characteristics.b) Comparison of the number of events in each Southern California sequence as a function of the magnitude of the largest event.Modified from Fig. 7 in Vidale and Shearer (2006).c) Same as (b) but for Japanese swarms from Fig. 2a of Vidale et al. (2006).d) Comparison of the swarm radius with the equivalent magnitude of the summed moment in the swarm for the Japanese sequences.Modified from Fig. 2b in Vidale et al. (2006).(For interpretation of the colors in the figure, the reader is referred to the web version of this article.)(Fig. 9b) point to intrusion of a dike trending NW-SE and located to the WSW of the swarm centroid.The orientation of the dikes coincides with the NW-SE-elongated pattern of earthquake epicenters with hypocenters dipping to the SE.Recent morphometric studies have also suggested that the more recent plumbing system at Makushin is dominated by N135 • -striking dikes based on the prevalence of dike occurrence in Quaternary deposits (Tibaldi and Bonali, 2017).In addition, Lerner et al. (2020) have found that low-magma-flux and young volcanic systems, such as Makushin, have less centralized magma storage and transport networks than systems with higher magma flux or longer durations of magmatic activity.In these systems, preexisting structural and upper-crustal stress heterogeneities may more significantly influence ascending magma (i.e., increased lateral transport along zones of weakness, such as faults or lithologic contacts).This is consistent with the 2020 earthquake swarm, which was offset from the main summit vent in proximity of the offset magma reservoir as imaged by tomography (Syracuse et al., 2015;Lanza et al., 2020; see Fig. 4) and exhibits hybrid and complex interactions of tectonic and magmatic signatures.

"Large" (M3+) magnitude earthquakes and dike intrusion
We have shown that the Makushin 2020 intrusion was possibly of rather small volume (no InSAR or GPS signal).However, the associated seismic activity was one of the strongest (with two M4+ earthquakes) recorded on the Island since 1996.Several recent studies aimed to relate seismic moment release to the intruded volume (e.g., White and McCausland, 2016).However, the predictive capability of these empirical relationships bears large uncertainties and systematically underpredicts volume changes observed during VT swarms near volcanoes (Meyer et al., 2021).Conversely, it has been shown that stress imposed by high rates of deformation due to rapid magma accumulation can promote largemagnitude VTs (and thus large cumulative seismic moment), despite a small intruded/erupted volume, especially when relatively large preexisting structures or planes of weakness are present (e.g., Nishimura et al., 2001;Wauthier et al., 2016).The latter two factors are clearly observed at Makushin; the swarm productivity is mostly concentrated in the first few days of the swarm, indicating rapid emplacement, and the initial MS-AS behavior can be attributed to the influence of pre-existing tectonic stresses.Moreover, all M2.5+ have FPS consistent with regional stresses and theoretically locate in zones of induced tensile stress change in Coulomb models.Effectively, the regional stress magnitude is enhanced in these zones, resulting in larger-magnitude earthquakes.Thus, we postulate that the difference in magnitudes of the seismicity attributed to the tectonic and magmatic processes can be related to the size of the pre-existing faults (relatively large and with fast rupture as shown by earthquake patterns in time and space) and the volume affected by the magmatic stress perturbation (relatively small as indicated by the lack of a detectable deformation transient).

Implications for earthquake sequence terminology and swarm classification
We have documented evidence in support of a magmatic intrusion to explain the 2020 seismic activity at Makushin Volcano.However, because the 2020 swarm started with the strongest magnitude earthquake, a typical characteristic of MS-AS sequences, in which aftershocks are the response to stress transfer from mainshock and the driving mechanisms are primarily tectonic, it poses implications for the currently accepted seismicity cluster classification scheme and casts doubts on the fundamental criteria that distinguish tectonic/MS-AS from magmatic seismicity.
Generally, following Vidale et al. (2006), the following characteristics indicate a sequence belongs in the swarm-like category: (1) the volume of activity is large compared to the sequence moment release; (2) there is a lack of correlation between the number of earthquakes and the magnitude of the largest one; (3) the occurrence of the largest event is usually well after the beginning of the sequence; (4) the rate of seismicity is steady after an initial burst; and (5) there is temporal expansion of the seismicity zone.All listed criteria except for point 3 can be applied to the 2020 seismicity burst at Makushin.On the other hand, mainshockaftershock sequences (1) begin with the largest event, the mainshock, (2) the presence of smaller foreshocks is frequent; (3) the activity is shorter in duration, and (4) the number of events scales with magnitude, as does the volume and duration of seismicity.Here, the 2020 seismicity only fits the first criterion.In the literature, it is acknowledged that hybrid sequences exist (Vidale et al., 2006).Nonetheless, these hybrid swarms still retain the characteristic of having the largest event well after the beginning of the sequence.Our analysis shows that we need to reconsider the definition of a swarm as having the 'largest magnitude in the middle of sequence' and possibly consider swarms holistically, especially in terms of duration, b-value, and focal mechanisms.

Summary and conclusions
The 2020 earthquake swarm at Makushin is a complex sequence both in terms of source processes and sequence statistics.We derived earthquake relocations, FPS, and model Coulomb stress changes to identify the triggering mechanisms of the observed seismicity.Several compelling arguments favor the hypothesis that the seismicity at Makushin is likely connected to a magmatic intrusion; they are summarized here: (1) The temporal duration as well as the scaling of the volume of seismicity and the number of events with magnitude are consistent with a swarm-type rather than a MS-AS-type sequence.Thus, a purely tectonic driving mechanism is unlikely.
(2) b-values are high.The elevated b-values in conjunction with the co-location of the swarm with low-velocity tomographic volumes point to the presence of a close-by active magma reservoir that generates a heterogeneous, intensely cracked crust, where, in turn, when critically and/or magmatically stressed, earthquakes can be triggered on the pre-existing fractures/faults.This is consistent with the observed complexity of the earthquake pattern and the absence of a clear migration.
(3) Many earthquakes have FPS with P-axes oriented are ∼90 • from the regional maximum compressive stress direction.
(4) Coulomb stress models indicate that FPS inconsistent with the regional stress field have a spatial distribution consistent with an inflating dike located WSW of the swarm and paralleling the NW-SE orientation of regional maximum compression.The occurrence of these earthquakes on pre-existing faults may also explain the complicated rupture pattern and the absence of hypocenter migration.This ensemble of seismically active fractures is in agreement with the swarm being located in an area close to the inferred magma reservoir at Makushin in a region of low Qp (150-250) and low Vp/Vs ratios (Lanza et al., 2020, point 2).
We conclude that magmatic processes are the primary driving mechanisms likely at play during the Makushin 2020 earthquake swarm, with tectonic processes playing an early and secondary role.Thus, we interpret the seismic swarm as likely the result of the incremental stresses from magmatic and tectonic processes on pre-existing faults, where the M3+ earthquakes are responsible for the initial aftershock-mainshock behavior observed in the data.The unusual characteristics of the Makushin swarm demonstrate that a re-evaluation of earthquake sequence classifications is warranted and paramount for guiding reliable interpretation and response during future swarms.

Fig. 1 .
Fig. 1.Maps (a, b) and cross-sections (c, d) showing Makushin Volcano and Alaska Volcano Observatory (AVO) catalog locations of 8,156 earthquakes from July 1996 to May 31, 2020 pre-swarm period (a) and 2,266 earthquakes from June 1, 2020 to December 21, 2021 post-swarm period (b) Note how activity is still ongoing in the region of the 2020 swarm (white dots).In b) yellow stars indicate the two M > 4 mainshocks from the 2020 seismic swarm.The dashed arrow in both cross-sections (c, d) shows the shallowing of the earthquakes from ∼15 km to 0-7 km depth closer to the volcano summit.(For interpretation of the colors in the figure, the reader is referred to the web version of this article.)

Fig. 2 .
Fig. 2. a) Magnitude (M L )-time plot of seismicity between July 1996 and December 2020.Notice that the events in June 2020 are the strongest recorded since 1996.The lowering of the magnitude threshold apparent in the last 5 years is due to a combination of network improvements together with improved data processing tools (Dixon et al., 2005; Power et al., 2020).b) Cumulative number of events for January 1, 2020 -December 2020.The red curve indicates the number of events according to an Omori law function, K /(c + t) with K ∼ 160 and c ∼ 0.1 and t time in days after mainshock.c) Helicorder record of the first 48 hours of the swarm recorded from station MGOD.Red lines indicate that the clip threshold was exceeded.(For interpretation of the colors in the figure, the reader is referred to the web version of this article.)

Fig. 3 .
Fig. 3. Frequency-magnitude distribution and associated variations in b-value versus time.a) Event moment magnitude versus time for 888 of the relocated catalog events (black dots).The blue line shows cumulative moment release with time for all events.b) b-value versus time calculated using the method described in Roberts et al. (2016).5000 random windows of sizes 50 to 500 are used to obtain the temporal b-value variation.PDFs are derived by stacking every 10 samples.c) Frequency-magnitude relationships for the 888 events from June 15, 2020 to December 7, 2020 for which we obtained a reliable moment magnitude Mw (blue curve, black dots) and for all the 1,154 relocated events using the AVO-calculated M L (red curve, red dots).Green circles are discrete values for the 888 events with reliable Mw.Magnitude of completeness (Mc) obtained with Mw is 1.56 and it is indicated by the blue dashed line, whereas the magnitude of completeness (Mc) obtained with M L is 0.36 and it is indicated by the red dashed line.(For interpretation of the colors in the figure, the reader is referred to the web version of this article.)

Fig. 4 .
Fig. 4. Map view and cross sections of earthquake relocations using hypoDD (red-filled circles).The AVO background seismicity from July 1996 -December 2020 is shown in gray.The white-filled star shows the location of the modeled inflating point source at ∼7 km depth from Lu et al. (2002), whereas the white-filled square and diamond indicate the centroids of the seismic tomography low velocity anomalies from Syracuse et al. (2015) and Lanza et al. (2020) respectively at depths of 8 km.The summit of Makushin is indicated by the white triangle in the upper left panel.The black hexagrams indicate the two largest events (M > 4.0).(For interpretation of the colors in the figure, the reader is referred to the web version of this article.)

Fig. 8 .
Fig. 8. Magnitude vs time vs FPS orientation.a) Entire swarm period, b) first month of swarm.as in Fig. 7. Horizontal black line indicates an apparent cutoff magnitude (M L 2.5) above which all but one earthquake have FPS consistent with regional stresses.(For interpretation of the colors in the figure, the reader is referred to the web version of this article.)

Fig. 9 .
Fig. 9. Coulomb stress change models for two end-member scenarios for the cause of swarm earthquakes.(a) Swarm earthquakes are the result of stress transfer (i.e., aftershocks) from the two M4+ earthquakes, and (b) swarm earthquakes are the result of stresses from dike inflation.In (a) modeled sources of Coulomb stress changes are the two M4+ earthquakes, represented by red boxes (fault geometry) and FPS showing the location and mechanism for the two earthquakes obtained in this study.In (b) the modeled source of Coulomb stress change is an inflating dike with position indicated by the green line.Green triangles represent the epicenter locations of earthquakes with inconsistent strike-slip FPS (from Fig.7).Upper-left inset illustrates the modeled receiver faults, which are idealized left-lateral strike-slip events with average fault orientation of inconsistent strike-slip events as shown in Fig.7.Small black symbols show epicenters of all catalog earthquakes in the study period.(For interpretation of the colors in the figure, the reader is referred to the web version of this article.)

Fig. 10 .
Fig. 10.a) Magnitude-duration relationships for the 2020 Makushin swarm modified afterHoltkamp and Brudzinski (2011) andIde et al. (2007).In addition to the Makushin swarm (red circle) we have plotted the Icelandic tectonic earthquakes swarm (TES) fromPassarelli et al. (2018b) and the swarmgenic slow-slip events (SSE) as shown inPassarelli et al. (2021).Duration is time between the first and last events in the swarm or sequence, and magnitude is the total seismic moment of earthquakes that comprise the swarm or sequence.(b-d) Characteristics of the 2020 Makushin seismicity compared against Japanese swarms and a of Southern California sequences as measured using parameters fromVidale and Shearer (2006)  andVidale et al. (2006).FollowingVidale et al. (2006)  terminology, "Average" indicates those sequences that do not fit aftershock-like or swarm-like sequences and have in between characteristics.b) Comparison of the number of events in each Southern California sequence as a function of the magnitude of the largest event.Modified from Fig.7inVidale and Shearer (2006).c) Same as (b) but for Japanese swarms from Fig.2aofVidale et al. (2006).d) Comparison of the swarm radius with the equivalent magnitude of the summed moment in the swarm for the Japanese sequences.Modified from Fig.2binVidale et al.  (2006).(For interpretation of the colors in the figure, the reader is referred to the web version of this article.)