Clustering of Long-Period Earthquakes Beneath Gorely Volcano (Kamchatka) during a Degassing Episode in 2013

: Gorely is one of the most active volcanoes in Kamchatka with a rich magmatic and eruptive history reﬂected in its composite structure. In 2013–2014, a temporary network of 20 seismic stations was installed on Gorely for one year. During the four months of its high degassing rate, seismic activity was mostly expressed in the form of a long-period (LP) seismic tremor. In this study, we have developed a workﬂow based on the combination of back-projection (BP), cluster analysis, and matched-ﬁlter (MF) methods. By applying it to continuous seismic records for the study period, we were able to identify discrete LP events within the tremor sequence automatically and individually investigate their properties. A catalog obtained using the BP detection algorithm consist of 1741 high-energy events. Cluster analysis revealed that the entire variety of LP earthquakes in this catalog could be grouped into ﬁve families, which are sequentially organized in time. Utilizing templates of these families in the MF search resulted in the complementary catalog of 80,615 low-energy events. The long-term occurrence of highly repetitive LP events in the same location may correspond to resonating conduits behaving in response to the high-pressure gases ﬂowing from the decompressed magma chamber up to the volcano’s crater.


Introduction
Processes in active magma systems can set in motion different types of seismic sources occurring either through abrupt fractures of rocks, or oscillations of magma containing reservoirs, or as a combination of these two processes [1,2]. These sources generate seismic waves that can be recorded by seismic stations and used to monitor volcano activity and to diagnose the state of the magma plumbing system. Unlike purely tectonic earthquakes in non-volcanic areas, the volcano-related seismicity has a broad range of types starting from volcano-tectonic (VT) earthquakes with clear arrivals of the P Eastern Volcanic Front, 25 km from the Pacific coast, it is related to the ongoing subduction of the Pacific Plate, which is located at a depth of~130 km below Gorely [22,23]. Morphologically Gorely is a compound shield-like stratovolcano with an altitude of~1800 m above sea level and a relative elevation of~850 m. Its upper part forms a linear northwest striking ridge of three merged primary cones and 11 superimposed summit craters complicated by more than 40 flank cones [24]. The modern Gorely edifice is located inside an ancient elliptic caldera with a size of 9 × 13 km, which is apparent on the topography map ( Figure 1).
Geosciences 2020, 10, x FOR PEER 3 of 19 Morphologically Gorely is a compound shield-like stratovolcano with an altitude of ~1800 m above sea level and a relative elevation of ~850 m. Its upper part forms a linear northwest striking ridge of three merged primary cones and 11 superimposed summit craters complicated by more than 40 flank cones [24]. The modern Gorely edifice is located inside an ancient elliptic caldera with a size of 9 × 13 km, which is apparent on the topography map ( Figure 1). The contemporary Gorely volcano represents the evolutionary development of an older volcanic center, followed by a radical transformation of its magma-feeding system [23]. Based on the age and composition of the erupted rocks, one can define three major stages of its formation [25,26]. The first (pre-caldera) stage is associated with the development of Pra-Gorely (also referred to as "Old Gorely"), which was an extensive (approximately 12 × 15 km in size) Middle-Pleistocene shield volcano stretched in the northeastern direction. Nowadays, the remnants of Pra-Gorely are mainly represented by peripheral parts of massive lava flows at the edges of the caldera and some relicts in the surrounding plateau [25]. The second stage led to the formation of a large caldera and massive felsic pyroclastic deposits in the surrounding area of 600 km 2 . There is debate about whether this thick ignimbrite and pumice complex with the total volume >100 km 3 has been deposited during a single [25] or multiple caldera-forming eruptions ranging in age from 361 ka to 38 ka [27,28]. Regardless of the eruptions number, such depletion of a large magma chamber embedded in the Earth's crust below Pra-Gorely caused its roof to collapse. Limited by steeply dipping arc faults, the Gorely volcano caldera is a typical collapse structure of the Krakatau type that is confirmed by magmatic permeability of individual sections in the caldera boundary [25]. The last (post-caldera) major stage started toward the end of Late Pleistocene with monogenetic volcanism on the weakened zone of the caldera rim. It was continued by the formation and development of the modern edifice in the central part of the caldera [26]. During Holocene Gorely's activity mainly consists of a cyclic The contemporary Gorely volcano represents the evolutionary development of an older volcanic center, followed by a radical transformation of its magma-feeding system [23]. Based on the age and composition of the erupted rocks, one can define three major stages of its formation [25,26]. The first (pre-caldera) stage is associated with the development of Pra-Gorely (also referred to as "Old Gorely"), which was an extensive (approximately 12 × 15 km in size) Middle-Pleistocene shield volcano stretched in the northeastern direction. Nowadays, the remnants of Pra-Gorely are mainly represented by peripheral parts of massive lava flows at the edges of the caldera and some relicts in the surrounding plateau [25]. The second stage led to the formation of a large caldera and massive felsic pyroclastic deposits in the surrounding area of 600 km 2 . There is debate about whether this thick ignimbrite and pumice complex with the total volume >100 km 3 has been deposited during a single [25] or multiple caldera-forming eruptions ranging in age from 361 ka to 38 ka [27,28]. Regardless of the eruptions number, such depletion of a large magma chamber embedded in the Earth's crust below Pra-Gorely caused its roof to collapse. Limited by steeply dipping arc faults, the Gorely volcano caldera is a typical collapse structure of the Krakatau type that is confirmed by magmatic permeability of individual sections in the caldera boundary [25]. The last (post-caldera) major stage started toward the end of Late Pleistocene with monogenetic volcanism on the weakened zone of the caldera rim. It was continued by the formation and development of the modern edifice in the central part of the caldera [26]. During Holocene Gorely's activity mainly consists of a cyclic alternation between phases Vulcanian-style explosive eruptions, voluminous (>0.1 km 3 ) lava flow eruptions, and intense degassing [25].
Having relatively high explosive eruption potential [29], Gorely may represent a significant hazard for aviation [30], tourists, and nearby infrastructures such as Mutnovsky Geothermal Power Plant (MGPP on Figure 1b) with a capacity of 50 MWt which provides a significant part of the electrical energy to the Petropavlovsk-Kamchatsky city and its surroundings. Therefore, the volcano has been thoroughly investigated by specialists in different disciplines of geosciences. The first robust data on the geological structure and development of Gorely volcano were published in [31,32], with the descriptions of the caldera, the associated pumice-ignimbrite deposits, the structural and material composition of the pre-caldera complex and the modern edifice. The comprehensive reconstruction of its Holocene activity via tephrochronological analysis was presented in [33]. In 1974-1977 a geological survey on a scale of 1:50,000 had resulted in a detailed geological map of the Gorely volcano [25]. In 20th century all eruptions (1921-1931; 1959-1960; 1980-1981; 1984-1986) were moderately explosive (VEI < 3) and occurred through the central summit with emission of basaltic-andesitic ash [24]. After the most recent one in 1986, a large fumarole was formed in the crater, through which an active emission of gases ensued. In the period of strongest degassing activity in 2010, the mass of gases emitted through this fumarole was estimated at 11,000 tons per day, with the outlet temperature reaching 900 • C. It is determined that these gases were composed of water (93.5%), CO 2 (2.6%), SO 2 (2.2%), HCl (1.1%), HF (0.3%), H 2 (0.2%), as well as some bromine and iodine compounds. It was estimated that under this regime, Gorely emitted about 0.3% and 1.6% of the total global fluxes from arc volcanism for CO 2 and HCl, respectively [34].

Seismic Data
In 1980, the Kamchatka Branch of the Geophysical Survey (KBGS) installed one telemetric seismic station GRL on the eastern slope of Gorely volcano. Two more permanent stations were installed in the summer of 2008 on the neighboring Mutnovsky and Asacha volcanoes. All stations were equipped with three-component sets of short-period channels based on SM-3 seismometers for recording the ground displacement velocity in the frequency band of 0.8-20 Hz. These permanent stations were used to investigate the seismicity beneath Gorely since 1984 [20,35]. Note, however, that these studies could only provide count and energy estimates for the events, but not the information about their locations.
A dense temporal seismic network ( Figure 1) consisting of 20 three-component broadband seismographs was deployed on the Gorely volcanic edifice and its surroundings in August 2013 by joint efforts of scientists from Trofimuk Institute of Petroleum Geology and Geophysics SB RAS (initiated this project), Department of Geology and Geophysics of Novosibirsk State University (provided seismic instruments) and Institute of Volcanology and Seismology FEB RAS (provided logistical support of the fieldwork). The network was removed in August 2014 and provided~350 Gb of continuous seismic records in total. Each of the temporal stations consisted of a CME-4311 (R-sensors, Moscow, Russia) three-component broadband sensor and a digital recorded Baikal-ACN-87/88 (R-sensors, Novosibirsk, Russia) with power supply provided by one box of 10 high-capacity power batteries Baken VTs-1 (UralElement, Verchniy Ufaley, Russia), external GPS antenna and necessary ventilated protection against dust and moisture. Baikal-ACN series recorders are three-channel autonomous seismic stations of an extended frequency range with an internal or external GPS module, a USB 2.0 channel for communication with a laptop, and a memory slot for SD card supporting volumes up to 32 GB. The CME-4311 three-component broadband velocimeter is built of three orthogonally oriented molecular-electronic transducer, and an electronic board, placed in a protective outer casing. The manufacturer stated flat instrument response in the frequency band of 0.016 (60 s) to 50 Hz.
We use the seismic records from a single permanent station GRL located on the volcanic edifice and maintained by the Kamchatka Branch of the Geophysical Survey as a reference. Even though ground waters eventually flooded some stations of the temporal array, overall data coverage remained consistent enough for more than eight months with maximal spatial density during the four starting months. It allows us to analyze the seismicity on Gorely during a significant period containing an episode of volcano's intense degassing. In this study, we consider the period from the 28th August till the 17th December, when the maximum seismic volcano-related activity occurred beneath Gorely. During this period, 18 stations of the temporary network were functioning, providing dense observation system on the volcano.
Preliminary analysis revealed that the dominant part of seismic energy was emitted by numerous LP earthquakes, which occurred on average twice per minute during most active phases of the degassing episode. Corresponding signals have a duration of about 10-15 s long, with an energy peak at 3 Hz. The strong similarity of these waveforms for consecutive events is the most notable feature of the dataset. Figure 2a presents an example of a five-minute seismogram of vertical components recorded by all available stations that clearly shows the LP swarm beneath Gorely. In this interval, at least eight events can be visually identified, and all of them have almost identical waveforms, as seen in an example in Figure 2b,c. four starting months. It allows us to analyze the seismicity on Gorely during a significant period containing an episode of volcano's intense degassing. In this study, we consider the period from the 28th August till the 17th December, when the maximum seismic volcano-related activity occurred beneath Gorely. During this period, 18 stations of the temporary network were functioning, providing dense observation system on the volcano. Preliminary analysis revealed that the dominant part of seismic energy was emitted by numerous LP earthquakes, which occurred on average twice per minute during most active phases of the degassing episode. Corresponding signals have a duration of about 10-15 s long, with an energy peak at 3 Hz. The strong similarity of these waveforms for consecutive events is the most notable feature of the dataset. Figure 2a presents an example of a five-minute seismogram of vertical components recorded by all available stations that clearly shows the LP swarm beneath Gorely. In this interval, at least eight events can be visually identified, and all of them have almost identical waveforms, as seen in an example in Figure 2b,c. A large number of events to be identified with a massive amount of the continuous seismic records urged us to develop a particular automated approach for compiling and analyzing the LP earthquakes. A robust catalog of specific type earthquakes is something highly desirable in seismology because it can give us valuable insights about the underlying mechanism and its evolution in time. We were trying to construct one by using three methods that successively built, investigate, verify, and enhance a catalog of the LP earthquakes. To build the first catalog with the most energetic events, we used an automatic detection algorithm based on back-projection (BP) technique. We then investigated the acquired catalog via cluster analysis that gave us a set of templates for several LPs clusters and served as a verification tool. Finally, applying a matched-filter (MF) technique, we searched for less-energetic events with waveforms similar to the obtained templates, thus compiling the complementary extended catalog of low-energy LP events. This A large number of events to be identified with a massive amount of the continuous seismic records urged us to develop a particular automated approach for compiling and analyzing the LP earthquakes. A robust catalog of specific type earthquakes is something highly desirable in seismology because it can give us valuable insights about the underlying mechanism and its evolution in time. We were trying to construct one by using three methods that successively built, investigate, verify, and enhance a catalog of the LP earthquakes. To build the first catalog with the most energetic events, we used an automatic detection algorithm based on back-projection (BP) technique. We then investigated the acquired catalog via cluster analysis that gave us a set of templates for several LPs clusters and served as a verification tool. Finally, applying a matched-filter (MF) technique, we searched for less-energetic events with waveforms similar to the obtained templates, thus compiling the complementary extended catalog of low-energy LP events. This catalog, however, should be used as a supplementary one in the analysis since it could be noise-contaminated.

Back-Projection Detection and Location Algorithm
The back-projection (BP) method is a practical approach to detect and locate seismicity by taking advantages of a seismic network or array. The core idea of of this method is a stacking of seismic records shifted by precomputed travel times to the theoretical origin points followed by a grid-search for the local maximum in space and time. Several shift-and-stack methodologies have been described in literature starting from "semblance analysis" [36] and "source scanning algorithm" [37,38]. In this paper, we use "beamforming" version of the BP method suggested for studying tectonic low-frequency earthquakes, which share many signal features with long-period volcanic earthquakes [39].
A seismic event originated in the location → x * at the time moment t * is recorded by a set of receivers located in → r i (i = 1, 2, . . . , N) as a set of waveforms u i (t). Assuming velocity model v for a study region, we can compute theoretical travel-times τ( x ) between each receiver and some virtual source location → x . The BP technique, which general concept is schematically demonstrated in Figure 3, is based on the stacking of normalized signal envelopes shifted in accordance to these precomputed travel-times. Thus, the recorded wavefields are kinematically projected back to the point → x . We are using a bending algorithm from the Local Tomography Software (LOTOS) [40] for raytracing and calculating of travel times, which gives us the potential to improve results confidence by using a more realistic velocity model. In this study, our primary goal is the detection and only then relative location of LP seismicity, thus a simplified model with a constant velocity value equal to 2 km/s was used. For the case of Gorely, this appears to be suitable because the LP seismicity is generated at shallow depths within the volcanic edifice. catalog, however, should be used as a supplementary one in the analysis since it could be noisecontaminated.

Back-Projection Detection and Location Algorithm
The back-projection (BP) method is a practical approach to detect and locate seismicity by taking advantages of a seismic network or array. The core idea of of this method is a stacking of seismic records shifted by precomputed travel times to the theoretical origin points followed by a grid-search for the local maximum in space and time. Several shift-and-stack methodologies have been described in literature starting from "semblance analysis" [36] and "source scanning algorithm" [37,38]. In this paper, we use "beamforming" version of the BP method suggested for studying tectonic lowfrequency earthquakes, which share many signal features with long-period volcanic earthquakes [39].
A seismic event originated in the location * ⃗⃗⃗⃗ at the time moment * is recorded by a set of receivers located in (i = 1, 2, …, N) as a set of waveforms ( ). Assuming velocity model v for a study region, we can compute theoretical travel-times ( , ) between each receiver and some virtual source location . The BP technique, which general concept is schematically demonstrated in Figure  3, is based on the stacking of normalized signal envelopes shifted in accordance to these precomputed travel-times. Thus, the recorded wavefields are kinematically projected back to the point . We are using a bending algorithm from the Local Tomography Software (LOTOS) [40] for raytracing and calculating of travel times, which gives us the potential to improve results confidence by using a more realistic velocity model. In this study, our primary goal is the detection and only then relative location of LP seismicity, thus a simplified model with a constant velocity value equal to 2 km/s was used. For the case of Gorely, this appears to be suitable because the LP seismicity is generated at shallow depths within the volcanic edifice. As described in the following paragraphs, for each point in the 3D space and time, we calculate a likelihood function that shows the coherency of the shifted waveforms across all receivers. Due to the small-scale media heterogeneity, the waveforms from an earthquake recorded at different stations are not coherent. That is why we ignore the signal phase and use a characteristic function (CF), ( ), instead of raw records. Common examples of CF are high-order statistics of the seismic signal As described in the following paragraphs, for each point in the 3D space and time, we calculate a likelihood function that shows the coherency of the shifted waveforms across all receivers. Due to the small-scale media heterogeneity, the waveforms from an earthquake recorded at different stations are not coherent. That is why we ignore the signal phase and use a characteristic function (CF), f i (t), instead of raw records. Common examples of CF are high-order statistics of the seismic signal (Kurtosis, Skewness) [41], short-term average to long-term average (STA/LTA) ratio [42], and signal envelope [43]. After many trials with different types of CF, we found that in our case, the most optimal form is an energy envelope or an absolute values of seismogram smoothed in a moving window: where τ h = 3 s is a half-size of the moving window for smoothing, and u i (t) is a horizontal component record of the i-th station. We define back-projection intensity (BPI) function p x and shifted according to travel-times between this point and each receiver: where N is the number of used receivers. The geometrical spreading A → r i , → x for a virtual source point → x and a receiver located in → r i is calculated as follows: where r 0 is the source size, which is approximated in our case by a unit sphere and d x is the length of the ray path in the reference model between the source and receiver.
For an actual source position (Figure 3b), the CFs calculated in Equation (2) are correctly shifted back in time and are stacked constructively, forming a maximum at the event's origin time. For any other points (Figure 3c), the same procedure will result in lower values of BPI. Computing BPI for a grid of virtual-source points → x j ∈ → X allows us to obtain an array of time-dependent functions p → x , t . Each of these functions represents a transformed wavefield kinematically projected back to a specific location inside the study region. They form a spatio-temporal distribution of BPI-p → x , t , which thereby may be considered as a time series of spatial images (snapshots) defined on the grid → X. Each snapshot p → x j , t depicts the likelihood of finding seismic source inside the study region at the specific time moment. In Figure 4, we present an example of a BPI snapshot corresponding to an average event based on the Gorely experimental data.
Incorporation of the BPI procedure inside a grid-search strategy is a core part of the BP detection technique. From full spatio-temporal BPI distribution p → x , t we construct a compressed BPI: and use it as a detecting function. By scanning through p c (t) in time for local maxima larger than a threshold value p d , we can effectively obtain time moments t max that correspond to the local maxima of the full BPI distribution. Then, we use a snapshot of BPI distribution at the t max moment for estimation of source location in 3-D space. Computing the time interval τ w between the absolute maximum to the nearest local minimum after the detected event gives us an approximate value of the signal duration. In the case illustrated in Figure 4, resulting length of the event's signal was approximately 16 s. For earthquakes located outside the study region, the snapshot maxima are usually observed on edges of the grid. We use this criterion in conjunction with the limitation on a signal duration τ w to exclude teleseismic and slab-related earthquakes out of the catalog. The BP detection technique was implemented to build a catalog of LP earthquakes beneath Gorely using two subsets of the whole data. Initially, with event detection threshold = 2 µ / , we performed detection for nine months of data available for five most consistent seismic stations (Figure 5a). This result showed that the significant of detections is condensed in the four starting months, with minor activity after middle December 2013. For this period of intense degassing, we were able to use data from 18 stations and obtain 9691 detections (Figure 5b). We then performed the BPI procedure with more conservative conditions presuming a more substantial value of the threshold equal to 4 µ / , which led to decreasing the catalog size to only 1741 most energetic events (Figure 5c). Spatial distributions of the detected events (third column in Figure 5) show that the significant part of the detected LP seismicity is located right beneath the volcano edifice. These results point out similar characteristics of seismicity in time for both low and high energy parts of the catalog. It is also apparent that after December 7th, the activity on Gorely is rapidly decaying. The BP detection technique was implemented to build a catalog of LP earthquakes beneath Gorely using two subsets of the whole data. Initially, with event detection threshold p d = 2 µm/s, we performed detection for nine months of data available for five most consistent seismic stations (Figure 5a). This result showed that the significant of detections is condensed in the four starting months, with minor activity after middle December 2013. For this period of intense degassing, we were able to use data from 18 stations and obtain 9691 detections (Figure 5b). We then performed the BPI procedure with more conservative conditions presuming a more substantial value of the threshold equal to 4 µm/s, which led to decreasing the catalog size to only 1741 most energetic events (Figure 5c). Spatial distributions of the detected events (third column in Figure 5) show that the significant part of the detected LP seismicity is located right beneath the volcano edifice. These results point out similar characteristics of seismicity in time for both low and high energy parts of the catalog. It is also apparent that after December 7th, the activity on Gorely is rapidly decaying.

Figure 5.
Results of applying the BP detection algorithm to Gorely data: (a) detections acquired using records of five temporary stations which functioned nine months; (b) same for 18 stations working during the four starting months, when degassing activity was exceptionally strong; (c) filtered detections obtained with 18 stations, which follow Gutenberg-Richter law and constitute BP-based catalog.

Cluster Analysis
The catalog of most-energetic events ( Figure 6) derived from the implementation of the BP technique can be further investigated in detail via cluster analysis. By dividing earthquakes into groups, we can effectively reveal the overall seismicity structure that can be related to the characteristics of underlying processes. To do so, one needs a well-defined principle to group 'similar' earthquakes together. Our approach for similarity quantification between pairs of events, in general, resembles the one used for swarms of repeating long-period earthquakes at Shishaldin Volcano in Alaska [44].

Cluster Analysis
The catalog of most-energetic events ( Figure 6) derived from the implementation of the BP technique can be further investigated in detail via cluster analysis. By dividing earthquakes into groups, we can effectively reveal the overall seismicity structure that can be related to the characteristics of underlying processes. To do so, one needs a well-defined principle to group 'similar' earthquakes together. Our approach for similarity quantification between pairs of events, in general, resembles the one used for swarms of repeating long-period earthquakes at Shishaldin Volcano in Alaska [44].
For each pair of detected earthquakes, the corresponding waveforms of vertical components are cropped in a time window τ w after the respective origin time (τ w = 16 s, in our case). We define the similarity of two earthquakes l and m as a correlation coefficient (CC) computed between the waveforms and mean averaged for all stations. To see the relationships between all events, one may plot these coefficients as a matrix, where each row or column reflects how similar the selected earthquake to the other ones in the catalog. Figure 7a shows the calculated CC matrix for the set of 1741 events identified for the Gorely volcano at the BP step with the higher threshold. For each pair of detected earthquakes, the corresponding waveforms of vertical components are cropped in a time window after the respective origin time ( = 16 s, in our case). We define the similarity of two earthquakes and as a correlation coefficient (CC) computed between the waveforms and mean averaged for all stations. To see the relationships between all events, one may plot these coefficients as a matrix, where each row or column reflects how similar the selected earthquake to the other ones in the catalog. Figure 7a shows the calculated CC matrix for the set of 1741 events identified for the Gorely volcano at the BP step with the higher threshold. Despite the general similarity of all LP events in the BP-based catalog, some subgroups of events can be distinguished from the visual analysis of the CC matrix in Figure 7a. The diagonal of this matrix represents auto-correlations. We can see from four to six square patterns of high CCs along diagonal, that indicates a group of similar events localized in time. Some of the 'squares' are prolonged to off-diagonal part, suggesting the existence of subgroups inside. The number of groups of events with different properties can be approximately identified by estimating the rank of the CC matrix, which can be done by computing its eigenvalues. The resulted distribution of eigenvalues ranged in the decreasing order (Figure 7b) shows that only a small part of the first eigenvalues was large enough, while the rest is close to zero. We decided to take into account only five eigenvalues that are larger than 5% of the maximum one, thus estimating N c = 5, the number of meaningful earthquake groups in the catalog. Figure 7c presents the mean average CC of every detected LP earthquake with all other 1740 events in the catalog. The majority of them have the CC between 0.3 and 0.4 showing relatively high similarity of all events. At the same time, there are a few events having the correlation of around 0.1. We manually inspected all these events and identified that they are relatively short-duration VT earthquakes having completely different waveforms compared to the LP events. It does not mean that no other VT events occurred during the studied period because the used parameters of the BP method were specially adapted for searching the LP events (reference velocity, time window, frequency of filtering). Thus, after excluding these four VT events, a verified initial catalog of the LP earthquakes contained 1737 events with the highest energy.
To separate all detected LP earthquakes in five groups, we followed the iterative approach described in [45]. To find out a reasonable partition of initial clusters, we consequently excluded groups of similar earthquakes from the catalog. First, we calculate a mean average CC for every single earthquake with all other events across the whole catalog, as shown in Figure 7c. Next, an earthquake with the maximum average CC is taken as a master event for the first group. All events having the CC with the master event larger than a particular threshold h are excluded from the catalog to form the first cluster. The procedure is repeated for the rest of the catalog: on each step earthquake with maximum average correlation is taken to form the next initial cluster that is excluded from the catalog until we get all N c clusters. After such selection, each master event is close to others in its group while staying far from the other master events. Depending on the chosen threshold, the entire catalog may be completely divided into N c clusters, or some earthquakes may stay ungrouped. In the second step, the defined clusters iteratively resorted in accordance with the CC matrix. On each iteration, we first check every grouped event and place it in the group where it has the highest similarity to the corresponding master event. After that, the new master event for each cluster is determined by computing a new vector of the mean similarity from a subsection of the CC matrix. As a result, stable cluster distributions are organized after several iterations. Since this method converges to a local minimum, the final result depends on the cluster's starting 'centers' (master events). Reasonable choice of the starting cluster composition via excluding them from the catalog helps us form highly diverse clusters. Despite the general similarity of all LP events in the BP-based catalog, some subgroups of events can be distinguished from the visual analysis of the CC matrix in Figure 7a. The diagonal of this matrix represents auto-correlations. We can see from four to six square patterns of high CCs along diagonal, that indicates a group of similar events localized in time. Some of the 'squares' are We performed described cluster analysis with the CC threshold of 0.3 and found that it converges to the stable distribution after the 3rd iteration. The CC values within the five groups and their time distributions are shown in Figure 8a,b. It can be seen that the CC values within distinct groups are higher than CCs with all events shown in Figure 7c, which demonstrate the adequacy of such classification of events. Final clusters are sequentially arranged in time that may indicate the possible evolution of seismic source properties or changes in the seismic velocity structure of the volcano.

Matched-Filter Detection Algorithm
The results of the cluster analysis allow us to create a set of templates that reflect common waveform features for all earthquakes in a certain group of the catalog. We create the cluster templates by stacking waveforms with the weights equal to their correlation coefficients. Thus, for the i-th receiver, the resulting template waveform ( ) of the k-th cluster is computed as follow: where ( ) is a waveform of the j-th event inside the cluster, denotes event-to-master CC, and Nk stands for the cluster size. Stacking increases signal-to-noise ratio so that we may treat the template as a fingerprint of a composite event with the common source mechanism for all events in the

Matched-Filter Detection Algorithm
The results of the cluster analysis allow us to create a set of templates that reflect common waveform features for all earthquakes in a certain group of the catalog. We create the cluster templates by stacking waveforms with the weights equal to their correlation coefficients. Thus, for the i-th receiver, the resulting template waveform u k i (t) of the k-th cluster is computed as follow: where u k i (t) is a waveform of the j-th event inside the cluster, C k n denotes event-to-master CC, and N k stands for the cluster size. Stacking increases signal-to-noise ratio so that we may treat the template as a fingerprint of a composite event with the common source mechanism for all events in the corresponding cluster. Since a template represents a generalized image of a cluster, we can compare it to another one visually and numerically by calculating CCs between them. In Figure 9, we show examples of the composite waveforms in some stations corresponding to the selected five groups of events.
Geosciences 2020, 10, x FOR PEER 14 of 19 Figure 9. Template waveforms for five closest to summit stations (same colors as Figure 8).
We can use the constructed templates to search for other earthquakes, similar to ones selected for the corresponding clusters. The matched-filter (MF) search, which we use to detect multiplets [42], has become the standard way to identify families of tectonic low-frequency earthquakes [46][47][48][49][50] and volcanic LP earthquakes [17]. The MF algorithm consists of a matching template to continuous seismograms by computing CCs between the template and the waveforms in a sliding window. Comparing the derived time-dependent CC with a predefined threshold cd, we can identify all events with the waveforms similar to the template, which creates a more extensive catalog than initially selected. Depending on the template-to-template similarity and chosen threshold, the same event can be detected with different templates. In these cases, we attribute such events to the template group having the highest CC. Unlike the BP technique, the MF approach is sensitive to the shape of seismic signals rather than its amplitude so that we can form a complementary catalog of weak LPearthquakes. For the Gorely case, we were able to identify 80,615 LP earthquakes divided into five clusters. The time distributions of these events can be compared with high-energy events of the initial catalog in Figure 10. We can use the constructed templates to search for other earthquakes, similar to ones selected for the corresponding clusters. The matched-filter (MF) search, which we use to detect multiplets [42], has become the standard way to identify families of tectonic low-frequency earthquakes [46][47][48][49][50] and volcanic LP earthquakes [17]. The MF algorithm consists of a matching template to continuous seismograms by computing CCs between the template and the waveforms in a sliding window. Comparing the derived time-dependent CC with a predefined threshold c d , we can identify all events with the waveforms similar to the template, which creates a more extensive catalog than initially selected. Depending on the template-to-template similarity and chosen threshold, the same event can be detected with different templates. In these cases, we attribute such events to the template group having the highest CC. Unlike the BP technique, the MF approach is sensitive to the shape of seismic signals rather than its amplitude so that we can form a complementary catalog of weak LP-earthquakes. For the Gorely case, we were able to identify 80,615 LP earthquakes divided into five clusters. The time distributions of these events can be compared with high-energy events of the initial catalog in Figure 10.

Discussion
Two implemented detection techniques gave us the possibility to carefully explore LP seismicity beneath Gorely volcano during a period of intense degassing. Our implementation of the BP approach is sensitive to the event energy, while the MF method is responsive to the waveform shape itself. By identifying signals of individual LP earthquakes in continuous seismic records of 18 temporary stations, we obtained the BP-based catalog containing 1737 high-energy events and the extensive MF-based list of 80,615 detections. It is important to note that the latter being complementary to the initial BP catalog is not entirely independent, as we are using templates constructed from the waveforms of identified events to obtain MF detections.
Cluster analysis of the BP-based catalog has demonstrated the limited variety and high repetitiveness of LP seismicity taking place beneath Gorely in the observation period. In total, we have identified five distinguishable families of LP earthquakes, which were sequentially arranged in time. In Figure 7, one can see the final distribution of high-energy LP earthquakes that gives us insight into the development of the conduit structure over time. At the beginning of the observation on August 28th, we see that the "red" family was dominating. By September 9th, its intensity has incrementally decreased, as it was gradually replaced by the "blue" family, which in turn reached the maximum on September 10th and then weakened during a couple of weeks in the same manner.

Discussion
Two implemented detection techniques gave us the possibility to carefully explore LP seismicity beneath Gorely volcano during a period of intense degassing. Our implementation of the BP approach is sensitive to the event energy, while the MF method is responsive to the waveform shape itself. By identifying signals of individual LP earthquakes in continuous seismic records of 18 temporary stations, we obtained the BP-based catalog containing 1737 high-energy events and the extensive MF-based list of 80,615 detections. It is important to note that the latter being complementary to the initial BP catalog is not entirely independent, as we are using templates constructed from the waveforms of identified events to obtain MF detections.
Cluster analysis of the BP-based catalog has demonstrated the limited variety and high repetitiveness of LP seismicity taking place beneath Gorely in the observation period. In total, we have identified five distinguishable families of LP earthquakes, which were sequentially arranged in time. In Figure 7, one can see the final distribution of high-energy LP earthquakes that gives us insight into the development of the conduit structure over time. At the beginning of the observation on August 28th, we see that the "red" family was dominating. By September 9th, its intensity has incrementally decreased, as it was gradually replaced by the "blue" family, which in turn reached the maximum on September 10th and then weakened during a couple of weeks in the same manner. Between October 3rd and 7th, there is a gap in the LP seismic activity. Between October 9th and 16th, there was a short-lived peak of the "green" family activity followed by another gap on October 16-18. On October 20th, the "orange" family abruptly started and produced the largest group of the LP events, which gradually increased until November 9th. After another short gap, the last "violet" family started on November 13th, reached its peak on November 22nd, and then completely decayed to December 6th. After this moment, LP activity only occurred infrequently.
The extensive MF-based catalog provides additional information about the evolution of Gorely's degassing system. As one may see in Figure 10, lower-energy events are repeating the same pattern as their high-energy counterpart. However, there are fewer gaps in the activity, and the transitions between LP earthquake families are much smoother. Other differences are the domination of the "blue" family over the "red" one, considerably prolonged duration of the "green" period and the minor burst of detections on December 16th.
Each LP earthquake family produced a series of nearly similar signals emitted from a localized source region. These source regions were located beneath the volcano summit at depths of less than 1 km below the surface. In the recent tomography study of Gorely [21], the top of the prominent anomaly, representing the magma chamber, was located at~2.5 km below the surface. The interface between very high Vp/Vs in the magma chamber and low Vp/Vs in the overlying carapace is interpreted as a level of the phase transition in the molten magma. We suppose that dissolved fluids were degassing at this stage due to lower pressure in the shallower part of the magma reservoir. The LP earthquakes identified in our study may indicate the following pathway of ascending high-pressure gases. The rapid degassing and fast dynamic propagation of gas bubbles through the conduit root may lead to self-sustained oscillations within the magmatic channel [51], generating the LP radiation recorded by seismic stations. It is possible that at some moment, the conduit structure changes, resulting in new oscillation parameters and characteristics of repeated LP events. Therefore, the observed evolution of LP earthquakes may reflect the structural changes in the shallow part of the volcano-magmatic system.
Alternatively, given generally shallow levels of LP radiation on Gorely, groundwaters may be involved in two-phase 'steam and water' resonator system in a similar way that was proposed for Ngauruhoe volcano in New Zealand [52,53]. The climate of Kamchatka with a heavy snow cover of volcanoes over a half-year and the glaciation of Pra-Gorely caldera both support this version. However, the chemical content of the and the rate of Gorely degassing implies that the proposed 'bubble-dynamic' mechanism of LP is primarily caused by dissolving magmatic fluids in the conduit root.
We see that only at the beginning of the observation period, two of LP earthquake families (the "red" and "blue" one) functioned at the same moments with about similar intensity of the high-energy events. In other periods, only a sole family is dominating at the time, which is apparent in the MF-based detection distribution ( Figure 10). It probably means that the preferable degassing regime of Gorely requires only one conduit acting at a time. Smooth transitions in the number of detected events between families ( Figure 10) and the structure of CC matrix (Figure 7a) imply that each dominant family gradually evolves into the next one. This may be interpreted as the slow migration of the source along a constricted pathway of magmatic gasses ascent. Close likelihood of template waveforms for master events of "blue", "green", "orange", and "purple" families also support this point, while "red" one represents alternative explanation. The co-existence of two LP families acting in the same period (but not simultaneously) may be interpreted as the balancing stage of the gas ascent process. In such conditions, the pressure gradient allows only a portion of gas bubbles to overcome constriction in the "blue" family origin point, while the rest of the gasses have to proceed laterally (Figure 8c) to the "red" family origin point.