Investigating potential icequakes at Llaima volcano, Chile

2 Glaciallyand magmatically-derived seismic events have been noted to heavily overlap in 3 characteristics, thus there exists the potential for false-alarms or missed warnings at ice-covered 4 volcanoes. Here we present the first study to specifically investigate icequakes at an ice-covered 5 volcano in Southern Chile. Two months of broadband seismic data collected at Llaima volcano 6 in 2015 were analyzed in order to quantify, characterize, and locate potential glacially-derived 7 seismic events at one of the most active ice-covered volcanoes in the region. We find over 8 1,000 repeating seismic events across 11 families, the largest of which contains 397 events. 9 Approximate locations and characteristics of the largest families lead us to conclude that these 10 events were derived from persistent stick-slip motion along the ice-rock interface at the base of a 11 glacier near the volcano summit. These results have implications for future seismic monitoring at 12 Llaima volcano and other ice-covered active volcanoes in the region. 13

. Map of Llaima volcano with the locations of the 2015 seismic stations used in this study marked with green diamonds (5 of the 26 stations are not visible). Names marked with asterisks (*) were those used in the STA/LTA method described in Section 3. Also marked are the mapped summit glacial areas marked as 'clear' or 'debris-covered' ice (white area). Thick and thin contours mark 500 and 100 m altitude intervals, respectively. Inset: Map of Southern Chile with the location of Llaima volcano (red triangle) and Santiago (red star, SG) marked. Also plotted are the locations of other ice-covered volcanoes within the Southern Volcanic Zone of Chile that have displayed eruptive activity in last 200 years (white triangles; Venzke, 2013).
The glacial area presented in this study (white area in Fig. 1) was calculated by using high-resolution 77 (0.5 m pixel size) Digital Globe panchromatic satellite image taken on March 6 2016. This image was 78 georeferenced to geographical coordinate system using WGS1984 datum yielding an estimated horizontal 79 accuracy of 1-2 pixels. The image was manually analyzed following the methods presented by Paul et al. 'unclear'. These categories were defined by characterizing the surface patterns where clear ice polygons 82 are bright snow or ice surfaces with or without crevasses and little debris cover material. The albedo is 83 high for snow and lower when bare ice is present. Sometimes the ground area is slightly brown due to 84 minor surface material. The 'debris-covered ice' class was assigned when surface debris form patterns 85 indicating that underlying ice was present. Such patterns could be formed by the presence of crevasses 86 ice-cliff backwasting, undulations, or morrenic-like alignments. Finally, unclear areas were selected when 87 surface patterns were similar to inactive rock glaciers but without crevasses or other features indicating ice 88 dynamic were observed. The spatial resolution of the utilized image (0.5 m) was very good for detecting 89 very small features like erratic blocks, small crevasses and other glacier origin forms. This resolution 90 allowed estimating the total glacial area on Llaima volcano with higher detail than previous work. This 91 explains why our estimate is significantly larger than 5.5 km 2 estimated by Reinthaler et al. (2019), that 92 used Landsat 8 OLI images (15 m pixel size for band 8) in which debris-covered ice would not be clear. 93 Nevertheless, it is clear from satellite images that the glacial area at Llaima volcano has been significantly 94 reduced in recent decades due to eruptive activity and global climate change (Reinthaler et al., 2019). 95 To provide a degree of security for nearby population centers, OVDAS (Observatorio Vulcanológico 96 de los Andes Sur 1 ) has deployed a network of stations around the volcano to continuously monitor its 97 activity. OVDAS use the criteria described in Lahr et al. (1994), Chouet (1996), and Chouet and Matoza 98 (2013) to identify and classify the earthquakes recorded by the seismic network around the volcano. Arrival 99 times and waveform amplitudes are used to differentiate between volcanic and tectonic events. Using a 100 reference station within the network, the volcanic earthquakes are classified as volcano-tectonic, long-period 101 (including hybrid), or tremor events. Each type of earthquake has been associated with multiple distinct 102 source mechanisms and relative temporal trends of each type has important implications for assessing the 103 activity state of a volcano (see Chouet and Matoza, 2013, and references therein). OVDAS also manually 104 distinguishes other non-volcanic or non-tectonic events such as cryogenic earthquakes, but have no mandate 105 to track these events therefore little is known about their prevalence in the seismic record (Mora-Stock et al.,  University and Southern Andes Volcano Observatory (OVDAS) collaboration (Fig. 1). Application of 118 receiver function analysis to this seismic data revealed a low-velocity zone at 8-13 km depth beneath the 119 volcano that could be interpreted as a magmatic body (Bishop et al., 2018). The network used a variety of 120 broadband seismometers that used various digitizers recording the data at 100 samples per second, see 121 Table S1 for specific details of what each station used.
To detect candidate seismic events at Llaima volcano, we applied a multistation detection algorithm on 123 seismic data collected from 1 February to 31 March 2015. Trigger times were extracted from multiple 124 stations using a short-term average/long-term average algorithm (STA/LTA), on condition that an event 125 was detected coincidentally in time at ≥2 stations. We used lengths of 1 and 9 s for the short and long 126 time windows, respectively, and a ratio of 5 was used to define a detected event; these parameters were 127 decided by manual inspection of events detected over 24 hours of seismic data recorded at station GEO.

128
Considering the low magnitude and strong attenuation noted for icequakes at other volcanoes (e.g. Allstadt 129 and Malone, 2014), we used only the eight closest stations to the summit for this step (marked by asterisks 130 in Fig. 1 and Table S1). Seismic data were preprocessed with a bandpass filter of 0.5-10 Hz to improve the 131 signal-to-noise ratio (SNR).

132
From the catalog of candidate triggers compiled by the multi-station detection algorithm, our next step 133 was to find seismic events that were repetitive over the period of study. In order to reduce the computing 134 load, we followed a similar methodology to that detailed by Allstadt and Malone (2014) who used an 135 algorithm modified from Carmichael (2013). The method uses unsupervised clustering of seismic events so 136 the user does not need to define templates in order to detect repeating events. First, we cross-correlate every 137 event with all other events within each day and group them into families, using a minimum cross-correlation 138 coefficient of 0.7 to define two events as a match; we use the scipy.cluster.hierarchy Python package to carry 139 out this step (see https://docs.scipy.org/doc/scipy/reference/ for more details, last accessed 7 February 2020).

140
For each event, we used the first 5 s of the waveform, sufficient to include the largest wave amplitudes   165 The earthquakes allocated to the largest family of repeating events (henceforth called Family 1) are 166 generally small, with magnitudes of less than 1, and appear to be of an emergent low-frequency nature ( Fig.   167 3a). However, the low-frequency and emergent nature of these events were likely the result of path effects 168 as the waves were strongly altered and attenuated as they traveled away from the volcano (Fig. 4a). To

178
Calculating the source locations for each of the families is crucial for understanding the source 179 mechanism(s) involved. However, locating individual events within each family detected at Llaima 180 volcano without unacceptable error margins is impossible due to the emergent and low SNR nature of 181 each waveform, as well as the rapid attenuation of the signal as it moves away from volcano (Fig. 4a).

182
Nonetheless, following Allstadt and Malone (2014), we can take advantage of the repeating nature of 183 these families and calculate median stacks for each family at each station. If there are enough events in the 184 family, clearer signals with relatively high SNR can be acquired on at least 3 stations in the network (Fig.   185 4b). The improvement in the SNR is such that relative P-wave arrival times across the station can be used 186 for a grid-search location algorithm. In addition, we can also determine the direction of first motions at 9 187 of the closest stations to the volcano summit for Family 1 (Fig. 5). The first motions for Family 1 in the 188 vertical component show mixed polarities across these stations. However, the stacking method only applied 189 for three families, as the SNR did not improve enough for clear P-wave arrivals in the remaining families.

190
Once the P-wave arrival times were picked, we used a brute-force 3D grid-search algorithm to estimate 191 source locations. This method uses the relative arrival times between the first recorded arrival and all In other words, artificial relative arrival times are calculated for each point in the grid, and compared to the 194 real relative arrival times. The location is that which most closely matches the real relative arrival times. 195 We defined the grid of source nodes using a 29 m horizontal and 37 m vertical resolution. A previous study where the mean-squared frequency, ω 2 , can be calculated from the seismogram data, u(t): 231 Lamb et al.

Llaima icequakes
where the integral is performed over a window of length 2T centered at time t andu is the time derivative 232 of the waveform, u. We also apply a correcting factor to R to account for bias due to noise in the waveforms 233 (Douma and Snieder, 2006). The relationship between the variance of the travel-time pertubation and 234 inferred source migration depends on the source mechanism, such as explosive, point, or fault-plane 235 (Snieder and Vrijlandt, 2005). Evidence from the mixed first-motion polarities (Fig. 5) suggest that it is 236 reasonable to assume, for the purposes of this calculation, that the source is dominated by shear motion 237 along a fault-plane. Therefore, if displacement occurs along a fault-plane, the source dislocation between 238 waveforms, δ, is given by: where v p and v s are the P-and S-wave velocities in the medium. Note that using different seismic volcanoes that approximately range from 1.5 to 2.5. Here we calculate displacements using P-wave 243 velocities ranging from 1 -4 kms -1 , with a v p /v s ratio of 2. As the individual waveforms within Family 1 244 have relatively low SNR, we instead apply CWI to stacked subsets of the family in order to improve the 245 SNR. Family 1, featuring 397 events, was divided up into 13 subsets of 30 or 31 events, and median stacks 246 were calculated from each stack. R was calculated using 8 second windows starting 5 s after the start of the 247 stacked waveform, bandpass filtered at 1-10 Hz, for each stack relative to the first stack, and converted to δ.

248
For Family 1, the calculated displacements from waveforms recorded at two different stations (HRD, 249 GEO) are <1 m/day (Fig. 7). The largest displacements appear to occur during the first part of the recorded  Métaxian et al., 2003;Jónsdóttir et al., 2009;Allstadt and Malone, 2014). Indeed, during our study period,

260
OVDAS officially cataloged no icequakes as it is not within their mandate to do so (Fig. S1). While we 261 study a relatively small time period, from our observations described above we would argue that glacially 262 derived seismic events may be far more prevalent in the seismic record than previously thought. 263 We conclude that the low-frequency and repetitive seismic activity detailed here is caused by glacial  stations, but first motions observed here are mixed (Fig. 5), and 2) a volumetric source would be expected 282 to generate low-frequency resonance and thus common spectral peaks would be seen at all stations (e.g.

283
Waite et al., 2008) and this characteristic is not observed here (Fig. 4c). It would be appropriate to carry 284 out a full moment tensor inversion to quantify the source mechanism but the lack of an accurate shallow 285 velocity model prevents us from doing so; further work is needed to calculate a reliable velocity model for suggests that failure in poorly consolidated materials such as ash tuffs is unlikely to generate repetitive 297 low-frequency seismicity (Heap et al., 2015). However, with the evidence presented here we cannot completely rule out shallow slow-slip as a potential source of minor seismic activity on other regions of the 299 volcano.

300
For glacial sources of seismicity, there are multiple different mechanisms that have been documented 301 (Podolskiy and Walter, 2016). We can disregard mechanisms involving hydraulic resonance in or below the 302 ice (e.g. Lawrence and Qamar, 1979;Métaxian et al., 2003) because there are no consistent spectral peaks 303 between stations or evidence of harmonics (Fig. 4c), though the resonant character of the signal could 304 be lost due signal alteration in the heterogeneous medium at shallow depths. Furthermore, we observe 305 mixed polarity first motions (Fig. 5) , 1979;Thelen et al., 2013). Again, the mixed polarity first motions present a strong 319 argument against crevassing as it is a volumetric source and should generate isotropic first motions. It is 320 worth noting that our analytical workflow made a key assumption that most of the icequakes that could 321 be occurring at Llaima are of a repetitive and persistent nature. It is possible that there were also many 322 small, non-repetitive seismic events of a glacial origin that were not automatically detected here. Outside 323 of manually and inefficiently picking these possible events from the seismic record, it is not yet feasible to 324 build a catalog of these events.

325
Of all the candidate source mechanisms, basal stick-slip sliding close to or at the interface between ice 326 and rock is the most likely. Repetitive, low-frequency seismicity generated by discrete glacial movements 327 along the base has been well documented (e.g. Weaver andMalone, 1976, 1979;Ekstrom et al., 2003;328 Caplan -Auerbach and Huggel, 2007;Zoet et al., 2012;Thelen et al., 2013;Allstadt and Malone, 2014). The 329 repetitive, persistent families observed at Llaima volcano ( Fig. 2c) require non-destructive and repeatable 330 sources, which can be provided by stick-slip motion over a stationary asperity at the ice-rock interface.

331
Alternatively, stick-slip motion can also be generated by rocks embedded in the ice (i.e. 'dirty patch'; e.g.

332
Allstadt and Malone, 2014) but the low or stationary motion of the source calculated from CWI (Fig. 7) 333 suggests the former is more likely. The mixed polarity first motions for Family 1 (Fig. 5) are also consistent 334 with shear failure at the source, in agreement with what is inferred to occur during stick-slip motion.

335
Stick-slip behavior requires two conditions be met: 1) friction must decrease with slip velocity, so that the 336 associated acceleration can be sustained, and 2) healing (i.e. strengthening) must occur at the slip interface, 337 so that static stress can be recharged (Zoet and Iverson, 2018). With the latter condition, one effect is 338 that longer time periods without slip would lead to bigger stress build-up and bigger subsequent seismic 339 events, a behavior that is hinted at for Family 1 (Fig. 3b) ice-covered volcanoes due to eruptive products such as tephra. Lastly, it is important to note that any linear 346 relationships between repose time and seismic event energies (e.g. Fig. 3b) could be explained by other 347 physical mechanisms. For example, it could also be indicative of repeated pressurization and failure of 348 a fluid-driven crack (e.g. Matoza and Chouet, 2010;Matoza et al., 2015). Therefore, while the tangible 349 relationship presented in Fig. 3b  These studies also noted that the base of a glacier is a dynamic environment with some time periods more 366 favorable for basal stick-slip behavior than other time periods. Thus, there is a good chance that the seismic 367 station network deployed in early 2015 were coincidentally in the right place at the right time to detect the 368 icequakes at Llaima volcano. As this study only looks at a relatively short two month period at the volcano,

369
it is clear there is a need to expand the analysis to a multi-year scale so that seasonal changes in glacial 370 seismic activity can be constrained. Furthermore, the locations calculated here would be of an unacceptably 371 low quality for the needs of continuous monitoring and risk assessment. Therefore, future deployments 372 at Llaima will need to explore new deployment configurations around the glaciers to help constrain the 373 source locations for such low energy events.

374
The findings detailed in this study have important implications for continuous monitoring at Llaima