Edinburgh Research Explorer Repetitive fracturing during spine extrusion at Unzen volcano, Japan

. Rhythmic seismicity associated with spine extrusion is a well-documented phenomenon at a number of dome-forming volcanic systems. At Unzen volcano, Japan, a 4-year dome-forming eruption concluded with the emplacement of a spine from October 1994 to February 1995, offering a valuable opportunity to further investigate seismo-genic processes at dome-forming volcanoes. Using continuous data recorded at a seismic station located close to the dome, this study explores trends in the seismic activity during the extrusion of the spine. We identify a total of 12 208 volcano-seismic events in the period between October 1994 and February 1995. Hourly event counts indicate cyclic activity with periods of ∼ 40 to ∼ 100 h, attributed to pulsatory ascent deﬁned by strain localisation and faulting at the conduit margins. Waveform correlation revealed two strong clusters (a.k.a. multiplets, families) which are attributed to fracturing along the margins of the shallow, ascending spine. Further analysis indicates variable seismic velocities during the spine extrusion as well as migration of the cluster sources along the spine margins. Our interpretation of the results from seismic data analyses is supported by previously published ﬁeld and experimental observations, suggesting that the spine was extruded along an inclined conduit with brittle and ductile deformation occurring along the margins. We infer that changes in stress conditions acting on the upper and lower spine margins led to deepening and shallowing of the faulting sources, respectively. We demonstrate that the combination of geophysical, ﬁeld and experimental evidence can help improve physical models of shallow conduit processes.


Introduction
Lava dome growth and collapse represent a major volcanic hazard globally (Sparks, 1997).The transition from effusive to explosive activity of a dome may be rapid and can generate destructive pyroclastic flows (Calder et al., 2002).This range of behaviour presents a significant challenge for forecasting and hazard mitigation, for example, the 1991-1995 eruption at Unzen, Japan (see Sect. 1.1).Continuous monitoring of growing and unstable lava domes is thus central to risk mitigation.The movement of magma during domebuilding eruptions is a rich source of geophysical signals, particularly earthquakes.Analysis of shallow (< 3 km) volcanic earthquakes has highlighted the importance of discerning their characteristics for the purpose of hazard assessment (Neuberg, 2000;Chouet and Matoza, 2013).
Volcanic earthquakes are conventionally classified as either high-frequency or low-frequency.High-frequency signals, or volcano-tectonic events, are typical earthquakes within a volcanic edifice or in the crust beneath it (Chouet et al., 1994;Lahr et al., 1994;Chouet and Matoza, 2013).Their occurrence may be related to the movement of magma (e.g.Wright et al., 2006) or gas (e.g.Hill, 1996), or can be the result of local-regional stress interactions (e.g.Roman et al., 2008).Despite extensive research into low-frequency volcanic seismicity, the interpretation of its source mechanisms and propagation effects is still controversial.Typical low-frequency events (1-4 Hz) are ascribed to volumetric sources involving motion of fluids within resonating volcanic conduits and dykes (Chouet, 1996).Earthquakes characterised by high-frequency and impulsive onsets with a low-O.D. Lamb et al.: Cyclic fracturing during spine extrusion at Unzen volcano, Japan frequency coda, known as hybrids, are often considered part of the low-frequency family of events (Green and Neuberg, 2006).Neuberg et al. (2006) suggested that the long ringing coda of these signals are generated by conduit resonance triggered by brittle fracturing of magma, a process that can be preserved in a conduit margin (Tuffen and Dingwell, 2005;Kennedy et al., 2009;Kendrick et al., 2012Kendrick et al., , 2014a)).Repeating low-frequency events have been interpreted as the result of stick-slip motion along the margins of a viscous magma plug forced upwards by ascending magma (Iverson et al., 2006;Lensky et al., 2008;Iverson, 2008;De Angelis, 2009;Power and Lalla, 2010;Varley et al., 2010;De Angelis and Henton, 2011), and possibly regulated by frictional melting (Kendrick et al., 2014b).Harrington and Brodsky (2007) noted that hybrid events recorded at different volcanic systems all exhibit corner frequencies and magnitudes consistent with brittle-failure sources.A recent study by Bean et al. (2013), employing a more proximal instrument network (< 500 m from volcanic vent), proposed that slow-rupture failure of weak volcanic materials can produce low-frequency events.
In this study, we present a detailed characterisation of the seismicity associated with the extrusion of a spine during the late stage of the 1991-1995 dome-building eruption at Unzen volcano, Japan.We use the results from our methods, together with previously published field and experimental constraints, to develop a conceptual model of the processes occurring during the magma ascent and spine extrusion.

Unzen eruption and spine extrusion
Following almost 200 years of dormancy at Unzen volcano, Japan, a new period of activity began with a swarm of volcano-tectonic earthquakes to the west of the volcano in November 1989 (Nakada et al., 1999;Umakoshi et al., 2001).The hypocentral distribution of volcano-tectonic events subsequently progressively evolved eastwards, towards the volcano at the surface (Umakoshi et al., 2001).Eventually, eruptive activity at the surface began in November 1990, with the first months dominated by phreatic and phreatomagmatic activity as magma approached the surface (Nakada et al., 1999).In May 1991 exogenic dome growth began, with a brief phase of spine extrusion, followed by a series of 13 viscous, dacitic lava lobes (Nakada et al., 1999).Vulcanian explosions occurred during the initial dome eruption stages in June and August 1991, when the effusion rate was at its highest (Nakada et al., 1999).The growing lava dome was structurally unstable and experienced repeated partial collapses that generated pyroclastic density currents, forcing the evacuation of over 10 000 local residents.This activity continued until February 1995, by which point effusion had ceased and seismicity returned to background levels (Nakada et al., 1999;Umakoshi et al., 2008).
Protracted spine growth defined the last 5 months of the eruption, beginning in mid-October 1994 and growing at an average rate of 0.8 m day −1 in November (Nakada et al., 1999;Yamashina et al., 1999;Hornby et al., 2015).Extrusion proceeded obliquely, with the spine inclined at an angle of 45 • towards the east-south-east (Kohno et al., 2008).The final dimensions of the spine were 150 m long, 30 m wide, and 60 m high (Nakada et al., 1999), with a record of ductile and brittle deformation preserved along the spine margins (Smith et al., 2001).The exposed shear zone bounding the upper margins of the inclined spine was up to 10 m thick and comprised of sheared dacite, cataclasite and breccias (Hornby et al., 2015), whereas the underside has less pervasive damage and shows slickensides on the surface (Calder et al., 2015).Seismic activity during the spine extrusion was relatively low and shallow (< 0.5 km; Umakoshi et al., 2008), but periodic increases in event output were observed synchronous with tilt oscillation on the volcanic edifice (Yamashina et al., 1999), constrained to originate from a fluctuating pressure source at 0.7-1.3km depth (Kohno et al., 2008)., 6, 1277-1293, 2015 www.solid-earth.net/6/1277/2015/ 2 Data and methods

Data source
The seismic data used for this investigation was collected by the Shimabara Earthquake and Volcano Observatory (SEVO1 ; Shimizu et al., 1992).Three stations were installed < 1 km from where the dome would eventually emerge, and all were equipped with a 1 Hz vertical component seismometer (see Fig. 3 in Umakoshi et al., 2008).Signals were telemetered to SEVO and recorded continuously with a sampling rate of 100 Hz.For the last 5 months of the eruption, the most complete data set was recorded by station FG1.The station also lies close to the spine (∼ 600 m) and was positioned in line with the direction of spine inclination.Here we analysed the complete high-resolution data set from station FG1 for our investigation.

Automatic event detection
Single-station detection was performed from 1 October 1994 to 28 February 1995 to extract all types of seismic events from continuous vertical channel data using a conventional short-term average/long-term average (STA/LTA) algorithm (Allen, 1978).Using a single station improves the chronology of the progression of seismic events as it includes those not large enough to be detected at multiple stations.However, it is well known that the STA/LTA technique is prone to false alarms or missed events, depending on the choice of parameters.To help alleviate this issue, we tested different STA/LTA parameters on a 24 h period of the seismic record.
Once triggers were identified, we used a multi-parametric filter to remove false triggers from the record.Details of the parameters for the STA/LTA algorithm and the trigger filter can be found in Appendix A. The STA/LTA calculation and waveform correlation below were implemented within the framework of a Matlab toolbox tailored for large data sets (GISMO; Reyes and West, 2011).The STA/LTA algorithm is widely used for automatic detection of volcanic earthquakes from a large data set (e.g.Matoza et al., 2009;Thelen et al., 2010;De Angelis and Henton, 2011;Ketner and Power, 2013).Amongst others, the single-station detection approach was successfully employed for the analysis of continuous seismic data recorded during the 2009 eruption of Redoubt volcano (USA) where it allowed a fine-scale view of variations in seismic behaviour as potential indicators of impending explosions (Ketner and Power, 2013).

Statistical analysis
Volcanic time series are inherently non-linear and can show cyclicity over a range of timescales (e.g.Odbert et al., 2014).Superposition of multiple cycles within a data set can make it difficult to observe and interpret individual signals.The Fast Fourier transform (FFT) offers an efficient means of examining the characteristics of a time series (Danielson and Lanczos, 1942) via the power spectral density (PSD) estimate, which highlights the strength of periodic components in the signal (Percival and Walden, 1993).Here the PSD is estimated using the multitaper method (MTM), demonstrated to be the most robust method when there is no prior knowledge of the signal-generating source (Thomson, 1982).
The SSA-MTM Toolkit presented by Ghil et al. (2002) was used to perform the spectral analyses.A detrending correction was used to prepare the data by rendering the time series approximately stationary, then either padded with zeroes at either end or truncated to a length of n 2 , for integer n, as required for FFT.The significance of spectral peaks was assessed against a statistical red noise model (Mann and Lees, 1996), considered appropriate for systems with background noise heavily influenced by long-term processes.It is impossible to completely characterise the nature of the noise without prior knowledge of the source; therefore, the red noise model acts only as a guide to interpretation.Here only peaks above the 99 % confidence threshold are considered significant for discussion.MTM analysis requires statistical stationarity over the whole data window, which is not a common feature of many geophysical systems and can result in spectra which are difficult to interpret.Short-term Fourier Transform (STFT) analysis calculates a series of PSD estimates using a moving window of specified length with results illustrated using spectrograms.An assumption of stationarity is only required within an individual sample window; therefore, spectrograms are useful for tracking changes in the spectral content of a time series (Odbert and Wadge, 2009).In other words, the STFT spectrogram is a series of PSD estimates overlapping each other, instead of a single PSD estimate as produced by MTM analysis.The choice of parameters (window length and overlap) is critical and has been optimised for our analysis.A window length of 336 h (14 days) with 99 % overlap provided the best compromise between achieving sufficient temporal resolution and maintaining a long-enough window for robust analysis.The frequency distribution of each window was normalised to unity in order to remove the influence of changes to absolute spectral power, thereby allowing direct comparison of the relative frequency distributions between contiguous windows.For the analysis, a high-pass Butterworth filter with a cut-off of 672 h was applied to the time series to enhance the clarity of shorter-period cycles of interest in the final spectrogram.Comparison between spectrograms generated from both non-filtered (Fig. S1 in the Supplement) and filtered data indicated that the use of the filter did not affect either the timing or the frequency of spectral peaks.This combination of methods has not been previously used at Unzen volcano but has been successfully applied to geophysical data sets from other volcanic systems (Odbert and Wadge, 2009;Nicholson et al., 2013;Lamb et al., 2014).
At Soufrière Hills volcano (Montserrat), this combination of methods identified cyclic patterns in tiltmeter displacement (Odbert and Wadge, 2009), SO 2 flux (Nicholson et al., 2013) and seismic event rates (Lamb et al., 2014).Lamb et al. (2014) also described similar sub-annual cyclic patterns in a data set from Volcán de Colima (Mexico), attributed to cyclic motion of magma within the system and expansion and contraction of a feeder dyke.

Waveform correlation
Assessment of waveform similarity is useful for investigating trends in earthquake activity within large data sets.Repeating waveforms are significant because they represent the product of earthquakes generated by the same source in the same location.Small, impulsive earthquakes that accompanied spine extrusion at Mount St. Helens from 2004 to 2008 were dubbed "drumbeats" due to their rhythmic occurrence (e.g.Iverson et al., 2006;Iverson, 2008).These waveforms were hypothesised to derive from "stick-slip" motion of the magma plug as it was forced upwards (Iverson et al., 2006;Iverson, 2008).The margins of this magma plug experienced strain localisation, and the surface textures (deformation features) of the extruded spine can provide the key to subsurface deformation that may help regulate magma ascent (Kennedy et al., 2009;Kendrick et al., 2012).Here we explore the role of repeating earthquakes as an indicator of the conditions along the margins of the spine at Unzen.
We use cross-correlation to measure the similarity of waveforms in the event catalogue constructed using the single-station detection method (Sect.2.2).To calculate waveform cross-correlation, we extract a 5 s window of data beginning at the picked arrival time.All waveforms are filtered with a 0.5-20 Hz bandpass Butterworth filter to minimise background noise and increase the signal-to-noise ratio.Five seconds of data following the arrival is sufficient to include the largest amplitude section of most waveforms while minimising the effect of background noise.Changes in window length of a few seconds have a minor influence on the calculated correlation coefficient (not shown here).In the chosen period of analysis, each waveform is compared to every other waveform and the pair is designated a correlation value between 0 and 1 (i.e. a value of similarity where 0 is completely different and 1 is identical).We identify clusters (repeating earthquakes, a.k.a.families, multiplets) using a hierarchial clustering method similar to that used by Buurman and West (2010).Based on visual inspection of the cluster tree, we choose a correlation threshold of 0.8 and a minimum number of five events to define clusters.The value is equal to or higher than in comparable studies (e.g. Green and Neuberg, 2006;Buurman and West, 2010;Rodgers et al., 2013); further discussion of the choice of parameters is found in Appendix B. Previously, waveform correlation has detected low-frequency clusters during swarms at Soufrière Hills volcano which were linked to a fixed shal-low conduit source in the volcanic system (Green and Neuberg, 2006).Buurman and West (2010) found clusters of waveforms closely tied to explosions at Augustine volcano (USA) and proposed that the method has considerable potential as an automated real-time volcanic explosion monitoring tool.

Singular value decomposition
Since the original data set is derived from only three seismic stations, it is not possible to accurately calculate locations and, thus, absolute magnitudes and earthquake energies.However, we can assess the relative magnitudes between events of a cluster because the entire cluster must derive from the same source (e.g. Green and Neuberg, 2006).To determine the relative earthquake amplitudes and relative time lags within the two largest clusters we apply a method based on standard matrix factorisation: the singular value decomposition (SVD).The method was proposed by Rubinstein and Ellsworth (2010) and used on clusters located around Parkfield, northern California.SVD estimates relative earthquake size for repeating earthquakes by taking advantage of the waveform similarity of all the events within the sequence.The method uses the entire waveform, which reduces the uncertainty in earthquake size and produces a measurement which is consistent from event to event within a cluster.The technique robustly computes earthquake size while exploiting the information from the whole waveform, unlike the parametric information used in most standard magnitude estimation methods.The SVD can be used to describe a matrix M (i.e. a group of horizontal vectors) in terms of two sets of basic vectors (input U and output V) and set of nonnegative singular values ascribed to the weight of the output-basis vector, as in the following relationship: where U is a matrix that maps the weight applied to each output-basis vector for each data vector.Each column − → U i of U gives the relative amplitude of the associated output-basis vector ( − → V i ) for each data vector.When the input data consists of repeating events, the similarity of their waveforms means that most of the data's power is explained by the first output-basis vector (a.k.a. the first principal component).If we assume that all variables except size and arrival times are exactly the same, the other output-basis vectors will only be describing background noise.This means that the relative amplitude, and therefore moment, of the seismograms is described by the first input-basis vector.However, this can only be true when the SVD is applied to seismograms filtered such that the frequency content is below the corner frequency of the event that is being examined.Rubinstein and Ellsworth (2010) demonstrate that estimating the moment using SVD has much less uncertainty than the whole-network estimates of moment based on coda duration.A more detailed method-ology, as well as a discussion of the associated errors and pitfalls of the method, is presented in Rubinstein and Ellsworth (2010).
To find the relative time lags of the earthquakes in each cluster, we employed a cross-correlation method using the first output-basis vector calculated by SVD as a reference.That is, all waveforms in each cluster are cross-correlated with the first output-basis vector, under the assumption that the first output-basis vector is effectively the mean waveform and representative of an unperturbed background state.The results of this method should reflect temporal variations in velocity relative to such unperturbed background levels.

Coda wave interferometry
While it is not possible to calculate source locations, coda wave interferometry (CWI) has the advantage in using repeating waveforms from a single station to provide an estimate of source separation.By exploiting the crosscorrelation properties of the repeating waveforms with the use of CWI, it is possible to separate temporal changes of the material and source properties from the spatial variation of the earthquake source.The migration of a seismic source results in a change in distance to scatterers in the surrounding medium, affecting the arrival times of the waves that interfere to produce the seismic coda.De Angelis (2009) used this method to infer a source migration of 235 m over 24 h for repetitive low-frequency earthquakes at Soufrière Hills volcano.The correlation coefficient, R, is related to the variance of the travel-time perturbation, σ τ , and to the frequency, ω 2 , according to the following relationship (Snieder, 2003): (2) The frequency can be calculated from the seismogram data, u(t): where the integral is performed over a window of length 2T centred at time t.The relationship between the variance of the travel-time perturbation and inferred source migration depends on the source mechanism.Snieder and Vrijlandt (2005) have demonstrated this relationship for different types of source such as a explosive, point or fault plane.If displacement occurs, for instance, along the fault plane, the source dislocation, δ, is given by where v p and v s are P -and S-wave velocities in the medium.We applied Eq. ( 4) to 36 h waveform stacks from clusters identified by waveform correlation (Sect.2.4), to maximise temporal resolution and reduce the computing power required.Waveform stacking also improves the signal-to-noise ratio and improves correlation.The correlation coefficient was calculated for different frequency bands (1-5, 1-10 Hz) and time windows (1-5, 7-11 s) between each waveform stack and the first stack, and converted to source displacement.The calculations were based on a P -wave velocity of 4.5 km s −1 and S-wave velocity of 2.6 km s −1 , consistent with measured velocities during the Unzen Scientific Drilling Project (Ikeda et al., 2008).These velocities are smaller than these measured in laboratory samples from Unzen (Scheu et al., 2006) and passive seismic source experiments at Unzen (Ohmi and Lees, 1995).They are also higher than those inferred during an active seismic source experiment deployed a few kilometres to the west of Unzen (Matsumoto et al., 2012); our chosen values represent a compromise between these measurements.However, it must be noted that these velocities do not account for macroscale fractures and fracture density along the seismicity pathways.Therefore our velocity values represent the upper limit of the range of actual velocities.To test the sensitivity of the equations, we used the higher velocities of Ohmi and Lees (1995) and find they do not change the results significantly (Fig. S2).Therefore Eq. ( 4) is sensitive to the seismic velocities employed, but uncertainties do not appear to affect the results dramatically.
The velocities represent a second-order effect in comparison to the scarce knowledge of the source dynamics.The CWI method relies on a series of simplifying assumptions: the velocity field is considered isotropic and the medium locally homogenous, mode conversion at the surface is ignored and the source separation must be less than a wavelength; a full discussion of the limitations of the method is presented by Snieder and Vrijlandt (2005).However, results from its application have been validated by comparison with other techniques such as the double-difference relocation method (Snieder and Vrijlandt, 2005).

Single-station detection
After filtering out false triggers, a total of 12 208 seismic events were detected; hourly event counts are plotted in Fig. 1.Two key observations of the data set are the following: 1.The event rate for all types of events begins with peaks of 10-20 events per hour in October 1994, and decreases to less than 5 events per hour by the end of the eruption in February 1995.

Fast Fourier transform
MTM analysis was carried out on the complete record of hourly event counts from single-station detection analysis to provide a first-pass assessment of the cyclic character of the data set (Fig. 2a).The PSD estimate reveals two peaks that appear significant above the 99 % noise confidence threshold: 40 and 24 h.However, analysis of the first and last 4 weeks of the data set (Fig. 2b, c) reveals the cyclic components are not wholly persistent throughout the data set.The 40 h cyclic component does not appear in the PSD estimate for the last 4 weeks of the data set (Fig. 2c).Although this cyclic instability brings into question the assumption of statistical stationarity (required by definition of MTM spectral analysis), STFT analysis only requires the time series to remain stationary within each window (described in Sect.2.3).We therefore used the STFT method to explore the temporal variability of the hourly event count in more detail.Inspection of the spectrogram from the hourly event count (Fig. 1b) highlights the following key observations: 1. Mid-October to mid-November 1994 is dominated by a strong ∼ 40 h cyclicity.This behaviour does not appear again in any other part of the time series.
2. From December to the end of January the cycle period "glides" from ∼ 50 to ∼ 100 h.

3.
The 24 h periodicity noted in the MTM analysis (Fig. 2) is present in the spectrogram, but is not clear against the background noise.

Waveform correlation
For waveform correlation, we focused on 5738 events detected during the period of 1 October to 15 November 1994 to explore the emergence of the spine and its strong 40 h pulsatory extrusion (Fig. 3a).The key observations are the following: 1. Two large clusters dominate the results of the waveform correlation.The clusters, henceforth known as cluster 1 and 2, are highlighted and their waveforms contain a strong low-frequency component (Fig. 3a, b and c).
3. Cluster 1 (487 events) and 2 (181 events) are significantly larger than the other 27 clusters detected during this period; by comparison the third largest cluster has 32 events.
4. The occurrence of both clusters also shows strong cyclicity, corresponding to the approximately 40 h cycles identified in the STFT analysis of the hourly event counts (Sect.3.2).
5. A total of 84 % of events in this period (1 October-15 November 1994) did not fall under our criteria of being part of a cluster; thus the clusters represent a fraction of the seismic record (Fig. 3d).

Cluster characteristics
For each cluster, the frequency index for each waveform (calculated for single-station detection filtering; Sect.2.2) is plotted over time (Fig. 4a).There is a clear division in the frequency component of clusters 1 and 2, with the latter containing waveforms with more low-frequency components.SVDs of the waveforms in cluster 1 and 2 show that, during their lifetimes, the relative amplitudes (and by extension, relative moment) remain substantially unchanged (Fig. 4b).In contrast, relative time differences for each cluster show opposite trends (Fig. 4c).Cluster 1 has a positive trend, indicating that the events are arriving sooner relative to the mean cluster waveform; whereas for cluster 2, a negative trend indicates that the events are arriving later relative to the mean waveform (Fig. 4c).CWI was applied to the bandpass-filtered stacked waveforms from each cluster.Displacements calculated in the 1-5 s window of the waveforms, filtered at 1-5 Hz, are relatively stable during the early stages of both clusters (Fig. 5b, e).However, the last week of each cluster is defined by displacements of over 100 m.The displacements calculated after the 5-10 Hz filter show greater scattering but smaller displacements (Fig. 5c, f).The codas (5-9 s) in both bandpass filters do not reflect these changes, with relatively stable movement throughout each cluster's lifetime (Fig 5b,  c, e, and f).This indicates little or no motion of the scatterers in the source-receiver path.

Cyclic variations in seismicity
Single-station detection with a STA/LTA algorithm detected a total of 12 208 volcano-seismic events between 1 October 1994 and 28 February 1995 at Unzen volcano (Fig. 1a).Our count is in good agreement with previous results published by Nakada et al. (1999) and Umakoshi et al. (2008).Furthermore, good visual correlation between the temporal pattern of earthquake rates and real-time seismic amplitude measurements reinforces our confidence in the results of the detection algorithm (Fig. S3).However, our rates are much higher than those of Yamashina et al. (1999), who used data recorded 5 km south-west of the vent and which likely did not detect many smaller events during the activity.A combination of MTM and STFT analysis on hourly event counts revealed a strong ∼ 40 h cyclicity during the extrusion of a spine (Fig. 1b).We also see a "gliding" in the cycle period from ∼ 40 h in mid-November up to ∼ 100 h at the end of January 1995.The ∼ 40 h periodicity is very similar to that previously described in seismicity and tilt data (Nakada et al., 1999;Yamashina et al., 1999).However, the ∼ 100 h cycle at the end of January 1995 is longer than the ∼ 60 h period described by Nakada et al. (1999).

O. D. Lamb et al.: Cyclic fracturing during spine extrusion at Unzen volcano, Japan
Cluster 1; 487 events Cluster 2; 181 events (d) Hourly event rates for all events (grey) and events in clusters 1 and 2 (black) during the period of 1 October to 15 November 1994.
Numerous studies of seismicity during dome-forming eruptions have described characteristic sequences of earthquakes that correlate with changes in surface activity (e.g.Chouet, 1996;Iverson et al., 2006;Neuberg et al., 2006;Buurman and West, 2010;Ketner and Power, 2013;Johnson et al., 2014).This has been modelled to reflect frictional processes along conduit walls as the ascending magma releases accumulated energy by faulting and/or pulsatory ascent (Denlinger and Hoblitt, 1999;Voight et al., 1999;Iverson et al., 2006;Neuberg et al., 2006;Iverson, 2008;Lensky et al., 2008;Moore et al., 2008;Kennedy et al., 2009;Kendrick et al., 2012;Pallister et al., 2013;Kendrick et al., 2014b;Scharff et al., 2014).Yamashina et al. (1999) proposed a cyclic upward movement of the Unzen spine due to periodic variations in pressure in the conduit, and Umakoshi et al. ( 2011) invoked a similar model to explain a short period 1-2 h cycle during an earlier phase of the eruption in May 1991.The increase in observed cycle periods between May 1991 and November 1994 is likely due to a decrease in the effusion rate from > 2 to < 0.3 m 3 d −1 through the same period (Nakada et al., 1999).This is a similar mechanism to that inferred by Voight et al. (1999) to explain cycles of a similar periodicity at Soufrière Hills volcano, and later successfully modelled by Costa et al. (2012).The gradual decline and eventual halting of effusion by February 1995 would also explain the observed "gliding" of the cycle periodicity.The short-term cyclicity at Unzen could also be explained by the presence of an elastic-walled dyke (e.g.Costa et al., 2007); indeed, a dyke was emplaced beneath the dome in the early stages of the eruption (Yamashina and Shimizu, 1999;Umakoshi et al., 2011).However, the calculated thickness of the dyke (13 m; Yamashina and Shimizu, 1999), coupled with the low extrusion rate in November 1994, would not be able to account for the cycle period observed here (see Fig. 6 in Costa et al., 2007).
Given the high ambient temperatures in the shallow magma column, local temperature increases during slip events are likely to produce pseudotachylyte (e.g.Kendrick et al., 2012Kendrick et al., , 2014a, b), b), greatly altering fault properties (Otsuki, 2003;Di Toro et al., 2006;Nielsen et al., 2010;Hornby et al., 2015;Lavallée et al., 2015).A recent study has shown that slip velocities as little as 0.1 m s −1 are sufficient to induce frictional melting (Kendrick et al., 2014a) -a threshold regularly met during faulting in the seismic events related to the spine extrusion at Unzen volcano (Hornby et al., 2015).However, it is difficult to infer the importance of mechanical contributions to the cyclicity detected above, as our description of these signal's trigger mechanisms remains incomplete.Here, using a combination of waveform correlation, SVD and CWI we have attempted to get a better understanding of this problem.

Repeating waveforms
Waveform correlation was carried out on the period around the extrusion of the spine instead of the entire period to  focus on repetitive events related to its upward movement.A total of 29 clusters were detected between 1 October and 15 November 1994 (Fig. 3b).The two largest clusters, which emerge coincidentally with the first observation of the spine (Yamashina et al., 1999), also feature a strong ∼ 40 h periodicity similar to that detected in our FFT analysis (Sect.3.2).The similar cluster waveforms (Fig. 3b, c) and the cyclic, non-destructive source character of the seismicity suggest a common process occurring within a similar setting.In addition, the opposing arrival polarities of each stacked cluster waveform (Fig. 3b, c) indicates opposite failure directions.Therefore, we propose that each cluster of low-frequency events derive from brittle shear failure on opposite sides of the ascending magma plug at Unzen volcano.It is the rheology of magma that helps drive this process; during ascent magma becomes increasingly brittle as it both vesiculates and outgasses (Caricchi et al., 2007).This behaviour, on short timescales in the upper conduit, provides exceptionally dynamic rheological conditions that favour strain localisation and failure (Lavallée et al., 2008), providing a hypothesised source for repeating volcanic earthquakes during the delicate interplay of magma across the glass transition (e.g.Neuberg et al., 2006).As magma undergoes a similar pressure temperature path through time, the conditions of failure would remain constant and failure would be achieved at a similar depth in the conduit, thereby generating seismicity from a recurring source (Thomas and Neuberg, 2012).Field examination of eruptive products elsewhere have exposed the importance of multiple fault processes (Tuffen et al., 2003;Tuffen and Dingwell, 2005;Kendrick et al., 2012Kendrick et al., , 2014b;;Plail et al., 2014), where magma fracture is followed by friction, inducing comminution, brecciation, cataclasis (Kennedy et al., 2009), frictional melting and viscous remobilisation near the conduit margin (Kendrick et al., 2014a).These findings support the view that the bulk of the magma is able to ascend as a plug, where brittle fracturing along the margins of the ascending spine is responsible for seismicity, including the repetitive "drumbeats" observed during the 2004-2008 eruption at Mount St. Helens (Iverson et al., 2006;Iverson, 2008;Kendrick et al., 2014b).These mechanisms have been observed or, in some instances, inferred at other volcanic systems such as Augustine (Power and Lalla, 2010), Volcán de Colima (Varley et al., 2010) and Soufrière Hills volcano (Green and Neuberg, 2006;Neuberg et al., 2006;Lensky et al., 2008;Hammer and Neuberg, 2009;De Angelis, 2009;De Angelis and Henton, 2011;Kendrick et al., 2014a).
At Unzen, the gas flux scaled relatively consistently with magma extrusion rate throughout the eruption, with no marked increase during spine growth (e.g.1995); therefore, we exclude the movement of hydrothermal fluids as a source for the repetitive seismicity at Unzen (see e.g.Waite et al., 2008;Matoza and Chouet, 2010).In contrast, evidence of brittle failure along the spine margins at Unzen is present in field observations and laboratory measurements.Parallel-plate viscometry experiments have indicated that, under the temperature and stress conditions in the conduit at Unzen, the magma's rheology would have induced brittle failure in regions of high strain rate (Goto, 1999;Cordonnier et al., 2009).Indeed, the brittle failure of magma at high temperature is integral to magma discharge (Lavallée et al., 2008;Tuffen et al., 2008;Smith et al., 2011) and occurs when the strain rate surpasses the melt's timescale of relaxation (Webb and Dingwell, 1990).In the crystalrich Unzen magma (56 % crystal content; Cordonnier et al., 2009) there is a wide transitional area between viscous flow and brittle failure, where cracks nucleate, propagate and coalesce as catastrophic failure is approached (Cordonnier et al., 2012).In the conduit this shear failure occurs at the margins (Lavallée et al., 2012a), creating a marginal damage zone that enhances permeability and hence the outgassing capability of the magma column (Watts et al., 2002;Gaunt et al., 2014;Kendrick et al., 2014b;Plail et al., 2014).It is thought that lava domes may evolve from endogenous to exogenous growth through the development of such damage zones (Hale and Wadge, 2008;Cordonnier et al., 2009); indeed, sheet-like layers of different porosities were described in rocks within the Unzen spine (Kueppers et al., 2005) along with a variety of cataclastic and brecciated fault rocks that preserve a record of brittle and ductile deformation along the margins of the spine (Smith et al., 2001;Calder et al., 2015;Hornby et al., 2015).

Source variations
Since our data set only includes data from three stations, it is not possible to calculate the location of each earthquake in the record, which would aid our understanding of the source processes for the clusters.Instead, we employed the SVD method (Rubinstein and Ellsworth, 2010) to calculate relative amplitudes (and by extension, relative moments) and relative time differences (Fig. 4).While relative amplitudes for both clusters remain the same throughout each lifespan (Fig. 4b), the relative time differences have positive and negative trends for cluster 1 and 2, respectively (Fig. 4c).This suggests changes in either seismic velocity or sourcereceiver distance.The presence of a cluster generally implies a relatively stable source volume for the earthquakes, allowing for displacement within a quarter wavelength (Geller and Mueller, 1980).However, De Angelis (2009) used CWI to estimate a source migration of 235 m for a cluster at Soufrière Hills volcano, Montserrat.Determining relative relocation for the cluster events at Unzen volcano is difficult due to the low number of stations available, but published locations have shown that the vast majority of earthquakes during our period of analysis occurred in the shallow subsurface (< 1 km below summit; Nakada et al., 1999;Umakoshi et al., 2008).
Using CWI, we analysed the stacked waveforms from cluster 1 and 2 to assess changes in relative source locations (Fig. 5).We measured displacements of up to 100 m from the original position towards the end of each cluster's lifetime.However, displacements measured in the codas are relatively small, remaining close to zero.This shows that the position of scatterers remained relatively fixed during this period.Additionally, the calculated displacements are sensitive to which bandpass filter is used (Fig. 5b, c, e and f).However, as each waveform stack indicates a low-frequency nature (1-5 Hz, Fig. 3b, c) we do not consider displacements in the higher bandpass filter as relevant.The results presented in Fig. 5 do not account for errors that are due to the choice of source models and seismic velocities (see Sect. 2.6).For the former, it is reasonable to expect the displacement occurs along a fault plane rather than an explosive or point source; no explosions were reported at Unzen during the spine formation (Nakada et al., 1999).Using seismic velocities consistent with those measured during the Unzen Scientific Drilling Project (see Sect. 2.6 here; Ikeda et al., 2008), we estimate that these displacements cannot explain the entire relative time difference trend seen in Fig. 4c.Therefore, we ascribe the relative time differences in both clusters (Fig. 4c) to source displacements and an increase and decrease in seismic velocities for cluster 1 and 2, respectively.Such characteristic opposite trends must derive from a common phenomenon with different conditions, as described by our conceptual model (Sect.4.4).

A conceptual model
The displacement in the cluster sources and changes in seismic velocity at the source requires a mechanistic understanding.There are multiple processes by which the seismic velocity in a volcanic edifice can change: damage in and close to the volcanic conduit, the presence of hydrothermal fluids, or changes in stress and temperature conditions (e.g.Scheu et al., 2006;Heap et al., 2015).Damage around the margins of spines is reported at Mount St. Helens (Kennedy et al., 2009;Kendrick et al., 2012;Pallister et al., 2013;Gaunt et al., 2014) and Unzen volcano (Smith et al., 2001;Hornby et al., 2015;Calder et al., 2015).It is highly unlikely that the source location for cluster 1 was becoming less damaged during the cluster lifespan, but variable loading of damaged rocks may affect material response; for example, as fractures perpendicular to the principal stress close, seismic velocity can be increased (e.g.O'Connell and Budiansky, 1974;Heap et al., 2014).On the other hand, it is likely that the source www.solid-earth.net/6/1277/2015/Solid Earth, 6, 1277-1293, 2015 location for cluster 2 was becoming increasingly damaged by both the thermal and cyclic stressing associated with the eruption (Kendrick et al., 2013;Schaefer et al., 2015;Heap et al., 2014).The frequency content of each cluster suggests the presence of fluids around the source regions, with higher fluid content around the source for cluster 2 (Fig. 4a), despite no reports of significant changes in the flow of hydrothermal fluids towards the end of the eruption (Hirabayashi et al., 1995).However, we do not believe it can explain the observed opposite trends in time lag for the events in clusters 1 and 2, though hydrothermal fluids may contribute towards mechanical weakening of the rocks (e.g.Heap et al., 2015).
We therefore propose a model where the extrusion of an inclined spine can lead to both variations in seismic velocities at the cluster sources and migration of each source (Fig. 6).
In this model the source regions for each low-frequency cluster are located on the eastern and western margins of the spine, for cluster 1 and 2 respectively.The spine, rising in a pulsatory manner, undergoes volatile exsolution leading to densification and a loss in buoyancy.The dense, eastwardinclined spine then exerts an increasing amount of normal stress on the eastern, lower margin, and a decreasing amount on the western, upper margin.The inclination of the spine is evident from the inclined distribution of seismic hypocentres (Umakoshi et al., 2001) and pressure sources beneath the volcano (Kohno et al., 2008).Experimental work on Unzen dome samples by Scheu et al. (2006) has shown that, under isothermal conditions, an increase in normal stress increases the P -and S-wave velocities.Each seismic event in the cluster is the result of a shear stress increase above the material strength or the frictional resistance on the surface of the spine, which are equally affected by normal stress (Byerlee, 1978).On the eastern, lower margin of the spine, an increase in the normal stress occurred throughout extrusion and cooling and densification resulted in a shallowing of the cluster source as the failure criteria were met at the same stress conditions at increasingly shallower depth.On the western, upper margin the opposite effect is seen, the source deepened to where failure conditions were met as the spine sunk away from the contact, unloading it from below and decreasing normal stress.Furthermore, a reduction in normal stress across the western margin could introduce a degree of aseismic slip, causing cluster 2 to cease first.Therefore, the response of the magma subjected to differing stress-strain conditions on the upper and lower spine margin could in part explain the difference in occurrence times and rates between clusters 1 and 2. Additionally, the locations of each cluster source on opposite sides of the plug explains the opposite arrival polarities seen in the cluster waveform stacks (Fig. 3b,  c).The lower-frequency content in the waveforms from cluster 2 (Fig. 4a) could result from greater fluid saturation due to degassing of the magma through the upper surface of the spine as the failure point propagates ever deeper towards the base of the magma plug.These stress variations form a minor contribution to the variations in seismic velocity in the source region for cluster 2, dominated by the decreasing normal stress due to the inclined spine.

Conclusions
We have characterised the seismicity of spine extrusion during the last stage of the 1991-1995 dome-forming eruption at Unzen volcano, Japan.Using a single-station detection approach, we identified 12 208 seismic events between October 1994 and February 1995.Two Fast Fourier transform analysis methods (multitaper method and Short-term Fourier transform) revealed strong ∼40 h cycles in hourly counts during the emergence of the spine.By the end of January 1995, the cycle "glided" to a ∼ 100 h period.The cycles are linked to pulsatory extrusion of the spine, in part governed by the frictional properties of the magma plug and driven by magma supply from below.Waveform correlation of the data set revealed two strong cyclic clusters of low-frequency volcanoseismic events during the spine extrusion, interpreted as repeated fracturing events at the spine margins.A combination of singular value decomposition and coda wave interferometry revealed changes in seismic velocities in the Unzen edifice during the extrusion.Our conceptual model proposes that late-stage densification of the inclined magma spine resulted in increasing normal stresses on the lower margin, serving to close preferentially aligned fractures, and the opposite effect on the upper margin, where pervasive damage resulted from unloading from below as the spine slumped.This results in increasing and decreasing seismic velocities respectively, as well as a migration of the source regions along the conduit walls as failure criteria are met.This study highlights the application of single-station detection, statistical analysis, waveform correlation, singular value decomposition, and coda wave interferometry to deduce subtle variations during eruptions.We anticipate that the application of such mechanism-based interpretation of geophysical signals, combined with field and experimental evidence, will help improve physical models needed to inform future assessment of hazards during lava dome eruptions.

Figure 1 .
Figure 1.(a)The hourly event record of events detected from continuous data recorded at station FG1 at Unzen from 1 October 1994 to 28 February 1995.(b) The STFT spectrogram for the hourly event count time series plotted in (a).Spectrograms are plotted from 20 to 672 h (4 weeks) cycle periods; the maximum defined by the high-pass Butterworth filter.The power spectral density of each window has been normalised to unity.The last 2 weeks of November 1994 is dominated by a peak cyclicity whose period (> 100 h) lies close to the size of the moving window (336 h); it is therefore considered a methodological artefact.

Figure 2 .
Figure 2. MTM spectra showing time-series power spectral density of the daily event counts of (a) the whole hourly event count time series (1 October 1994-28 February 1995), (b) the first 4 weeks (1-29 October 1994), and (c) the last 4 weeks of the time series (31 January-28 February 1995).The power spectral density is plotted against various confidence levels of the red noise model.Peaks exceeding at least the 99 % confidence level are annotated with corresponding cycle period in hours.

Figure 3 .
Figure 3. (a) Catalogue of cluster occurrence in our data set from 1 October to 15 November 1994.Each plotted point represents an earthquake, and earthquakes on the same line are part of the same cluster.Only clusters with five or more events are shown.The two largest clusters, clusters 1 and 2, are plotted with red squares and blue triangles respectively.(b)The waveform stack of all events in cluster 1, overlain on a frequency spectrogram of the stack.(c) The waveform stack of all events in cluster 2, overlain on a frequency spectrogram of the stack.(d) Hourly event rates for all events (grey) and events in clusters 1 and 2 (black) during the period of 1 October to 15 November 1994.
Figure 4. (a)The frequency index of earthquakes in clusters 1, 2, and 3-29.(b) Relative amplitudes of earthquakes in clusters 1 and 2, as measured using singular value decomposition.(c) The relative time difference (i.e.time lag) of each earthquake in clusters 1 and 2, as measured using cross-correlation against the first principal component.In each panel the linear best fit for each cluster is plotted.Relative amplitudes and time differences were not calculated for clusters 3-29.

Figure 5 .
Figure5.Stacked waveforms (left column) and source displacements between the first stack and all other stacks for data bandpassed between 1-5 Hz (middle column) and 5-10 Hz (right column).The top and bottom row are for cluster 1 and 2 respectively.

Figure 6 .
Figure 6.(a) An E-W cross section of the conceptual model for spine extrusion and cluster source locations at Unzen volcano in October/November 1994, looking north.The inset shows a map view of the spine at the surface of the lava dome.Squares indicate the locations of (b) and (c).(b) Detailed view of the source region for cluster 2, showing how it moves progressively deeper over time (t 1 , t 2 , t 3 ) in response to the decrease in normal stress.(c) Detailed view of the source region for cluster 1, showing how it moves progressively shallower over time (t 1 , t 2 , t 3 ) in response to the increase in normal stress.Please see Sect.4.4 for more details and explanations.