Persistent shallow micro-seismicity at Llaima volcano, Chile, 11 with implications for long-term monitoring

21 Identifying the source mechanisms of low-frequency earthquakes at ice-covered volca- 22 noes can be challenging due to overlapping characteristics of glacially and magmatically 23 derived seismicity. Here we present an analysis of two months of seismic data from 24 Llaima volcano, Chile, recorded by the permanent monitoring network in 2019. We ﬁnd 25 over 2,000 repeating low-frequency events split across 82 families, the largest of which 26 contains over 200 events. Estimated locations for the largest families indicate shallow 27 sources directly beneath or near the edge of glaciers around the summit vent. These 28 low-frequency earthquakes are part of an annual cycle in activity at the volcano that is 29 strongly correlated with variations in atmospheric temperature, leading us to conclude 30 that meltwater from ice and snow strongly aﬀects the seismic source mechanisms related to 31 glacier dynamics and shallow volcanic processes. The results presented here should inform 32 future assessments of eruptive potential at Llaima volcano, as well as other ice-covered 33 volcanoes in Chile and worldwide.


Introduction
network. Window lengths of 0.7 and 8 s were used for the short-and long-term windows, 140 respectively, with a ratio threshold of 3.5 used to define a detected event at each station. 141 These parameters were decided using manual inspection of events detected over 24 hours 142 of seismic data recorded at station AGU and differ from those used in 2015 (minimum 2 143 stations, 1 and 9 s for short-and long-term windows, and a ratio threshold of 5). Seismic 144 data during this step were pre-filtered with a 1-10 Hz bandpass filter to improve the   Each event is cross-correlated with all other events within each day, using a minimum 156 cross-correlation coefficient of 0.8 to define 2 events as closely matching; this is higher 157 than the 0.7 used for repeating events in 2015. The first 5 s of each event waveform is 158 used, which is sufficient to maximise the SNR of each event. As station AGU had the 159 highest number of detected events, waveforms from this station were used to build the 160 catalogue of families. Families of repeating waveforms were defined using a hierarchical 161 clustering method similar to that used by Buurman and West (2010) and Lamb et al.

162
(2017); the scipy.cluster.hierarchy Python package is used for this step (For more 163 details, see https://docs.scipy.org/doc/scipy/reference/, last accessed October 2021). In 164 this approach, branches within the hierarchy are joined at nodes whose height is the mean 165 cross-correlation value between each event pairs spanning the two groups. These nodes 166 may join individual events or between clusters of events, depending on which linkage has 167 the highest mean cross-correlation, and families are defined by nodes whose values are 168 higher than a threshold. 169 Next, for each day a median waveform stack is computed for each family of 2 or 170 more events, which are then compared with all other stacks across the whole time period 171 to find larger, multi-day families. The last step ensures a more complete repeating 172 event catalogue by using a frequency-domain approach with an overlap-add method, (1) where x grid is a target source position, C ij is the normalized cross-correlogram between

211
The above calculation was performed within the 'enveloc' Python package (Wech 212 and Creager, 2008). This optimization problem was performed on a grid with 0.005 • 213 lateral and 100 m depth intervals (down to 5 km) centered on the summit of the volcano.

214
However, the 'enveloc' package was originally developed for searching for locations across 215 a much larger area than that which is used here so topography was previously not taken  Table S1 in 249 Supplementary Materials.

250
To estimate source amplitudes for each of the regional events, we first remove the site 251 amplification and instrument response, before filtering to 1 -10 Hz. The maximum of 252 the smoothed Hilbert Transform was used as maximum amplitude A i at each station i.

253
The amplitude at the source, A 0 , was then calculated at each station based on (Battaglia 254 and Aki, 2003) where r is the source-to-receiver distance, f is the central frequency, β is wave velocity   , and 73; Figs. S11, S13, S17, S20 and S21, respectively). The remaining 12 families 290 are all located at shallow depths (<500 m) around the summit of the volcano (Fig. 5).

291
Note that 5 of these families are located at or within close approximation to each other  could be calculated (Fig. 5). Local magnitudes for all located repeating events fell 298 within the 0.9 -1.5 range, with little distinct difference between families (Fig. 6). As we 299 Figure 3: Catalogue of family occurrence in our dataset. Each plotted point represents the time of an event, and lines join events from the same family. The total number of events in each family is noted with grey numbers before the first event. Families containing 30 or more events are highlighted using blue diamonds for the individual events.
assumed straight wave propagation between source and receiver, we underestimate A 0 for 300 regional events which, in turn, implies that local magnitudes for all repeating events are 301 overestimated; therefore, the values calculated here represent a maximum feasible value. to be continuing in 2019, albeit with higher numbers of detected events (Fig. 2a, b).  and depths for most families support the hypothesis that these were generated by glacial 341 activity rather than volcanic.

342
Qualitative event magnitudes for each located family suggest little distinct differences 343 in energies between families (Fig. 6). The lack of obvious differences between each family 344 suggests they may be generated by a similar source mechanism, but at different locations.

345
A review of glacial seismicity has suggested that source mechanisms may be identified

406
Only temperature has a peak (0.32) which exceeds the maximum correlation value that 407 could be obtained randomly (0.3). There is also a strong diurnal cycle in event rates 408 from 2013 to 2020 with lower detection rates from 1000 to 1800 local time (Fig. 7d).

409
Visual comparison with hourly wind and temperature data suggests this is possibly due 410 to increased seismic noise from wind drowning out smaller magnitude events (Fig. 7d). A