Seismic tremor reveals active trans-crustal magmatic system beneath Kamchatka volcanoes.

The occurrence and the style of volcanic eruptions are largely controlled by the ways in which magma is stored and transported from the mantle to the surface through the crust. Nevertheless, our understanding of the deep roots of volcano-magmatic systems remains very limited. Here, we use the sources of seismovolcanic tremor to delineate the active part of the magmatic system beneath the Klyuchevskoy Volcanic Group in Kamchatka, Russia. The tremor sources are distributed in a wide spatial region over the whole range of crustal depths connecting different volcanoes of the group. The tremor activity is characterized by rapid vertical and lateral migrations explained by fast pressure transients and dynamic permeability. Our results support the conceptual model of extended and highly dynamic trans-crustal magmatic systems.

INTRODUCTION

While most of volcanoes are fed by basaltic magmas originated in the upper mantle, the style and frequency of eruptions and composition of erupted rocks vary very strongly between different volcanic systems. This variability is due to the transport of magma from its generating region in the mantle toward the surface. This transport is particularly complex across a thick and relatively light continental crust where large volumes of magma can be stored at intermediate depths and exhibit chemical transformations via crystallization, degassing, differentiation, and interaction with surrounding rocks. However, our knowledge about the deep structure of the plumbing systems of active volcanoes remains limited, and, in many publications, their representation remains at a conceptual level. A traditional concept is that magma storage occurs mainly in relatively shallow and long-lived reservoirs. Geological evidence accumulated during the past three decades, however, spawned the image of trans-crustal systems where mixtures of melt, crystals, and volatiles are heterogeneously distributed with intermittent vertical connectivity (, ). The interaction between shallow reservoirs and deeper storage might affect the size and evolution of eruptions (). How these elaborate systems get built over time and how they contribute in individual eruptive events remain poorly understood.

Seismic () and sometimes electromagnetic () imaging are widely used to study the structure of volcanic systems and understand their dynamics. These approaches are mostly successful either in imaging shallow magma chambers in the upper crust or large-scale magma and melt containing regions in the uppermost mantle. Volcano-related structures within the midlower crust remain a difficult target for geophysical imaging and could only be reported beneath some of the largest or well-instrumented volcanic systems (–). This difficulty may reflect the distributed nature of magma storage in the midlower crust, resulting in much weaker contrasts in media properties compared with localized magma chambers. Here, we suggest that the deep parts of active volcano plumbing systems can be delineated on the basis of their seismic activity and apply this idea to one of the World’s largest volcanic groups located in Kamchatka, Russia.

Different physical processes occurring in volcanic systems can generate seismic signals (). First, volcanic seismicity can originate from faulting within the volcanic edifice and the underlying crust that generate so-called volcano-tectonic (VT) earthquakes (). This type of seismicity mostly occurs in the upper crust and often delineates geologic structures (rift zone, crustal detachment, etc.) existing in the vicinity of volcanoes rather than magma/melt-containing regions (, ). Also, rheological conditions leading to the brittle faulting are not often realized in the on average ductile midlower crust resulting in much fewer VT earthquakes occurring at depths below 10 to 15 km.

One of our main interests here is to study the deep roots of the magmatic system. Therefore, we focus on seismic waves generated by pressure perturbations within volcano plumbing systems. Resulting signals can be very variable including relatively impulsive events and long-duration tremors. Their distinctive feature is the depletion in high frequencies in comparison with VT earthquakes. Therefore, they form a large class known as volcanic long-period (LP) seismicity (, ). While other mechanisms have been hypothesized to generate very shallow volcanic LP earthquakes (, ), we will not consider them here because of our focus on the plumbing system beneath the upper crust.

Therefore, our main working hypothesis is that if we detect a LP seismovolcanic event and locate its source, then this points to an activated part of the plumbing system. Therefore, by building a representative catalog of locations of LP sources, we can delineate the extent of the active plumbing system. In addition, this approach can be used to study the internal dynamics of the volcanic system over time. While LP seismicity sometimes takes form of short impulsive signals (LP earthquakes), most often long-duration tremors are observed. LP seismicity is often observed close to the surface; however, it has also been reported to originate below the upper crust (–). This deep LP (DLP) seismicity might be suitable for illuminating the deep roots of plumbing systems.

We study the Klyuchevskoy Volcanic Group (KVG) located in the Russian Kamchatka peninsula. This is one of the World’s largest and most active clusters of subduction volcanoes. During recent decades, three volcanoes erupted: Klyuchevskoy, Bezymianny, and Tolbachik (Fig. 1A). Seismic tomography (, ) highlights the existence of a ∼30-km-deep magmatic reservoir that is possibly connected to active volcanoes through a rifting zone developed after a recent subduction reconfiguration (, ). This large-scale faulting structure can channel fluids and transfer pressure resulting in LP seismicity observed beneath the KVG at a broad range of depths (, ).

Fig. 1.

Spatial density of the tremor sources delineating the active trans-crustal plumbing system beneath the KVG.

(A) Map of the studied region. The inset shows the location of Kamchatka peninsula. Seismic stations are shown with white and black inverted triangles for permanent and temporary networks, respectively. Black dashed rectangle indicates limits of the area used in the grid search. Red dashed line shows a profile for the cross section shown in (C). (B) Sum along vertical lines projected on the horizontal plane. (C) Sum along horizontal lines with fixed depths and along-fault distances projected on a vertical plane corresponding to the profile shown in (B). (D) Similar to (C) but projection on a perpendicular profile shown in (B). Tremor hypocenters are shown with open circles. Light magenta crosses in (B), (C), and (D) show the DLP swarms that occurred on 23 November 2015, 26 November 2015, and 25 February 2016.

KVG LP seismicity is dominated by tremors that are difficult to analyze with standard methods and require using data from seismic networks. Therefore, we used an approach based on the network covariance matrix (a Fourier domain equivalent of an ensemble of interstation cross-correlations) to detect coherent seismic signals and locate their sources (–). The parameters of the analysis were tuned to better detect deep sources of LP tremors. Therefore, we analyzed frequencies between 0.5 and 5 Hz and selected rather long time windows (20 min). Despite this, the method sometimes detected relatively strong impulsive events, some of which could not be associated with LP activity, such as the low-frequency part of radiation of VT earthquakes (). Therefore, we use a simple classification scheme to identify such impulsive sources and to verify their possible bias in the final results.

We applied our method to 45 stations in the vicinity of the KVG (Fig. 1A). We selected 13 stations from the permanent network run by the Kamchatka Branch of the Russian Geophysical Survey (KBGS) and 32 stations from a temporary network installed between August 2015 and July 2016 in the framework of the “Klyuchevskoy Investigation - Seismic Structure of an Extraordinary Volcanic System” experiment (KISS) (). The use of this dense seismic network allowed us to obtain a catalog of LP tremors with a three-dimensional (3D) spatial and temporal resolution that has not been previously achieved in the KVG.

RESULTS

We analyzed 44,020 20-min-long windows overlapping at 50% and computed three network-based characteristic functions (figs. S1 and S2) that have been used to detect and classify coherent signals (figs. S3 and S4). Processes and criteria selection are summarized in the flowchart (fig. S5). As a result, we selected 13,027 windows identified as tremors or swarms of LP earthquakes. For every window, we backprojected the cross-correlations in a 1D model (fig. S6) () to estimate a 3D spatial likelihood function of the source location (fig. S7).

Spatial density of tremor sources delineating the active trans-crustal plumbing system beneath the KVG

To characterize the overall spatial distribution of the seismic sources corresponding to active parts of the plumbing system, we stacked all these likelihood functions. The result shown in Fig. 1 (B to D) demonstrates that the highest source density is observed beneath the main part of the KVG including Klyuchevskoy, Bezymianny, and Ushkovsky. Also, some notable amounts of sources are located beneath Tolbachik and Udina. This suggests that the active plumbing system extends approximately 50 km from southwest of Tolbachik volcano to northeast of Klyuchevskoy volcano, i.e., along the rift structure (, ).

The vertical cross section shown in Fig. 1C suggests that the active plumbing system extends through the whole crust with a “teapot” shape where the Moho-level magmatic reservoir is connected to the Klyuchevskoy volcano via a main conduit and to the Tolbachik volcano and recently active Udina volcano () via a branch deviating at ∼25 to 30 km depth (i.e., near the crust-mantle boundary).

The analysis of location uncertainties (figs. S8 and S9) shows that the errors increase with depth and that they are larger in the vertical direction than horizontally (which is due to the method based solely on S-waves). Also, a comparison with some known LP earthquakes (fig. S9) indicates that depths might be systematically overestimated. With taking into account these uncertainties, the overall southwest-northeast ∼50 km horizontal extend of the tremor source area is well above the location error. At the same time, we cannot exclude that the apparent widening of the “main conduit” beneath Klyuchevskoy with depth might be due to the increase of the errors.

In the vertical direction, this main conduit is clearly separated in two groups. First, many very shallow tremor sources are detected above 8 km. The second large group of sources extends through depths between 10 and 50 km with density maximizing around 30 to 35 km (i.e., near the crust-mantle boundary), where the DLP earthquake activity has been previously reported (). The penetration of this area toward depths significantly larger than 35 km is most certainly caused by the location errors and bias. At the same time, its overall vertical extent over more than 20 km and existence of tremor sources above the previously known DLP activity cannot be explained only by vertical location errors, which systematically decrease from ∼10 to ∼5 km between depths of 35 and 10 km.

The “gap” in detected LP sources around 10 km depth may be either caused by the influence of the 1D velocity model on tremor source locations or due to some real difference of rheological properties such as accumulation of magma or saturation in gas as suggested by previously imaged seismic structures at these depths ().

We also verified that the imaged distribution of source density does not significantly depend on the selection of detected signals. Figure S10 shows the distribution in which the sources classified as “earthquakes” (i.e., impulsive individual events) were added. The only notable difference is the appearance of very shallow events northwest of Ushkovsky volcano that might be related to the presence of actively moving glaciers in this region (). At the same time, the imaged active plumbing system in the midlower crust remains essentially unchanged.

Spatiotemporal patterns of KVG tremors

The space-time distribution of the tremor hypocenters (Fig. 2 and movie S1) is not fully random and is mainly composed of short bursts. During many bursts, the plumbing system is activated at a large range of depths spanning the whole crust.

Fig. 2.

Spatiotemporal patterns of the KVG tremors.

(A) Tremor depth as function of time. Daily activity level determined by the KBGS () for three main volcanoes (KLU, Klyuchevskoy; BEZ, Bezymianny; TOL, Tolbachik) are indicated at the top with colored circles. Orange color code corresponds to ongoing eruption. (B) Along fault distance as function of time. Orange and green dashed lines indicate positions of Klyuchevskoy and Tolbachik volcanoes, respectively. Green points correspond to Tolbachik tremor source location for which the difference in temporal tremor stability between Tolbachik and Klyuchevskoy is greater than 0.05. Light magenta crosses in (A) and (B) show the DLP swarms occurred on 23 November 2015, 26 November 2015, and 25 February 2016. (C and D) Detailed view on depth-time tremor patterns during periods indicated with orange rectangles in (A). Source positions are represented with open circles plotted on top of spatial likelihood function smoothed in time. Yellow and gray dashed lines represent visually identified upward and downward migration episodes, respectively. Corresponding speed estimations are written in gray rectangles.

A detailed view reveals migration-like patterns when the activity moves very rapidly from depth toward the surface (Fig. 2, C and D). We visually identified different upward migration episodes and approximately represented them with simple linear trends (yellow dashed lines in Fig. 2, C and D). The resulting migration velocities were approximately estimated between 3 and 10 km/hour. Sometimes, weaker patterns with smaller number of sources and larger errors might indicate downward migrations. We represent these possible episodes with gray dashed lines in Fig. 2 (C and D) and show that the estimated velocities are in the range of those found for vertical upward migrations.

While, most of the time, the tremor activity occurs approximately beneath the Klyuchevskoy volcano, during a few episodes it migrates laterally toward Tolbachik (Fig. 2B). Sources classified as earthquakes (shown as blue circles in fig. S4) mostly correspond to additional episodes of DLP earthquakes and also to the previously mentioned shallow earthquakes in the vicinity of Ushkovsky. The sources of this type do not affect the observed tremor migration patterns.

DISCUSSION

The obtained spatial distribution of the LP sources allows us to image an extended trans-crustal plumbing system beneath the KVG, and the temporal evolution provides useful information on the processes of destabilization and activation of the volcanic system before eruptions. The regional tectonic stresses dominated by the east-west extension (, ) facilitate the vertically efficient magmatic system. From August to November 2015, no observable volcanic activity was reported in the KVG and most of the tremors were observed at depth, while the shallow part of the system remained relatively quiet. On November 23 and 26, two particularly strong DLP swarms occurred (indicated with magenta crosses in Fig. 2, A and B). A recent study () has demonstrated that this type of seismic activity may be due to deep degassing caused by an influx of fresh basaltic magma from the mantle. This, in turn, may act to destabilize the whole plumbing system. In the beginning of December 2015, the DLP swarms were followed by tremor activity over the whole range of crustal depths and numerous upward (and possibly downward) migrations (Fig. 2, C and D) between the Moho and the surface. Simultaneously, the shallow part of the plumbing system beneath Klyuchevskoy became very active, as shown by the increase in number of detected tremors and the shifting to yellow daily activity level (Fig. 2A). While this first activation episode did not end with an eruption, it continued until mid-February 2016 when the activity partially migrated toward Tolbachik. Then, after a short quiescence, another DLP swarm likely marking the injection of a fresh magma pulse occurred on 25 February 2016 followed by a subsequent destabilization phase that led to an eruption in April.

These results reveal some key characteristics of the KVG plumbing system and help us to understand how it evolves before eruptions. First, we obtained a nearly continuous spatial distribution of tremor sources over a large depth range (Fig. 1, B to D). Second, the spatiotemporal patterns of tremor sources constrain the physical mechanisms controlling the evolution of the activity. Rapid migrations (Fig. 2, C and D) cannot be explained by magma upwellings through dikes because the observed migration speeds (up to 10 km/hour) are faster than those associated with dyke propagation or other volcanic fluid migrations that have been reported of the order of 0.1 km/hour or smaller (–). Possible downward migrations indicated by some source patterns also could not be explained by the dyke propagation.

At the same time, these fast speeds are compatible with the very short delays between DLP swarms and the onsets of shallow activity that were observed during the build-up of the 1991 Mount Pinatubo eruption (). These migrations look rather like random diffusions that may be related to the physics of pressure transients in a fluid-saturated permeable medium. The implication, then, is that the KVG lies above a magma-filled network, which is in agreement with the conceptual model of extended trans-crustal plumbing systems ().

Pressure diffusion in a hydraulically connected magmatic conduit has been previously evoked by Shapiro et al. () to explain the link between DLP and shallow LP earthquakes observed at the KVG when the activity migrated over ∼30 km from the Moho level up to the crust during a few months. This corresponds to an average migration rate of ∼1 km/day or an average conduit diffusivity of ∼100 m2 s−1. With our new high-resolution tremor source catalog, we observe much faster migration episodes (up to 10 km/hour; Fig. 2, C and D). These rapid migrations can arise on top of slower diffusion processes as “pressure waves” in a system with pressure-dependent permeability (). This leads to a nonlinear diffusion equation that, for example, has been solved to model rapid migrations of tectonic tremors in a subduction zone ().

Trans-crustal volcano-plumbing systems are formed over long periods of time by multiple episodes of magmatic intrusions and remobilizations (, ), resulting in a strong compositional and mechanical heterogeneity. Thus, transport properties are likely to be highly variable in both space and time. In these systems, magma conduits likely contain many barriers that impede magma motion and lead to the buildup of large pressure gradients. These barriers act as valves, i.e., open and close under certain thresholds (–), resulting in rapid decompression events that are known as sources of LP seismic radiation (, ). In crystal-laden trans-crustal plumbing systems (), rapid changes of transport properties can occur because of sintering mechanisms reducing permeability in fluid pathways. As a result, low-permeability barriers can emerge. Solid-state sintering causes large reductions in porosity and permeability over days to months (). While the precipitation of mineral phases can induce sealing with a time scale ranging from hours to years (, ), viscous sintering can reduce the connected porosity within seconds to hours (). Intermittent flow can also arise when a magma with strongly temperature-dependent viscosity is transported through a network of reservoirs and conduits ().

In a conduit with multiple low permeability barriers, the overall permeability is proportional to the number of open versus closed valves (), which, on average, increases with overall pressure gradient across the system. A recent study of dynamics of such a system () has shown that closely located valves can interact, which results in cascades forming very rapidly migrating pressure transients (analogous to nonlinear pressure waves) able to propagate in both upward and downward directions, similar to what we observe for the KVG tremors (Fig. 2, C and D). The modeling shows that in a conduit with background diffusivity of ∼100 m2 s−1 [as previously inferred for the KVG ()], rapid migration velocities can reach ∼10 km/hour (), which is again close to the observations presented here.

We conclude, therefore, that the large trans-crustal KVG plumbing system remains in a hydraulically open state and is subject to very frequent fast fluid pressure transients emerging from a local valving behavior (). Because pressure controls most key volcanic processes such as fluid flow rates, mineralogical reactions, and degassing, its fast migrations provide a very efficient mechanism of destabilization/activation of the system before eruptions.

Several important implications for the structure and behavior of the KVG magma supply system can be drawn from the data of this paper. One is that this system stretches across the whole crust and that Klyuchevskoy is fed by dominantly vertical flow over the whole depth range. Another implication is that the supply system extends horizontally over a large distance encompassing both Tolbachik and Klyuchevskoy. Strong activity bursts propagate horizontally away from Klyuchevskoy but do not necessarily lead to a Tolbachik eruption.

While Tolbachik remained mostly inactive during the campaign operation period, our detailed tremor analysis suggests the existence of a secondary conduit that connects this volcano to the main plumbing system beneath Klyuchevskoy. This secondary conduit may play an important role during the large eruptions of Tolbachik (the last one in 2012–2013) when its average permeability rises allowing a discharge of a substantial part of fluid pressure from the main KVG channel. This large-scale pressure migration may explain the perturbation of “regular” cycles of activity of Klyuchevskoy and Bezymianny volcanoes revealed on the basis of analysis of seismological and satellite data ().

Tremor activity beneath the old Tolbachik (∼300 ka) is much weaker and scattered than beneath the much younger Klyuchevskoy [8 ka; ()]. Two factors may be at play. One is that the older Tolbachik rests over a mature network of magma pathways that are mostly open, whereas the much younger Klyuchevskoy depends on a supply system that is still under construction. Another factor is that the deep magma source has migrated laterally away from Tolbachik, such that the largest amount of magma is now generated beneath Klyuchevskoy. In this case, contrasting levels of tremor activity may be attributed to differences of magma supply rate. In both cases, however, it is clear that the magmatic system extends over large distances horizontally and vertically.

Our findings show that analysis of tremors based on seismic networks with appropriate coverage and density can be used to delineate the active parts of volcano-plumbing systems including those in the midlower crust and to follow their hydrodynamic evolution in time. Results of this approach with the KISS data reveal a large trans-crustal plumbing system beneath the KVG that connects the main active volcanoes of this group whose interactions () can, therefore, be explained by pressure transients at different time scales. Similar approach could be applied to other large volcanic complexes [or to different types of hydraulically controlled natural systems ()], either to those with very sustained seismic activity or when long time series of seismic network data are available.

MATERIALS AND METHODS

Data preparation

We use the data of 45 seismic stations shown in Fig. 1 from August 2015 to June 2016. The preparation and preprocessing of the data are done with the CovSeisNet software (). First, seismograms are aligned in time and stored in 24-hour-long segments. After down-sampling them at 25 samples/s, we apply demeaning, linear detrending, and band-pass filtering between 0.1 and 10 Hz. Instrument responses are considered being stable in time and not corrected for.

Network-based characteristic functions

We use three network-based characteristic functions described below to detect signals corresponding to LP tremors and earthquakes. The first two functions are sensitive to the coherence of the long-duration signals, i.e., mostly to tremors. The third function has been designed to be more sensitive to the presence of impulsive signals such as individual earthquakes and swarms.

Covariance matrix spectral width

For this characteristic function, we preprocessed the time series according to () with the application of a spectral whitening followed by a temporal normalization where the whitened signal is divided by a 0.25-s long running average of its absolute value. Following (), daily traces are cut into a number of overlapping time windows called averaging windows. Then, each averaging window is cut into M overlapping subwindows of length δt such that Δt = Mrδt, where Δt is the length of the averaging time window and r is the overlapping ratio that is 0.5 on our case. The array data vector is obtained from the N seismic stations u(f) as u(f) = [u1(f), u2(f), …, u(f)] and we then compute the outer product u(f)u†(f) within each subwindow, where † stands for the Hermitian transpose. The covariance matrix (f) is lastly estimated on any averaging window from M consecutive subwindowswith u(f) the data Fourier spectra vector in the subwindow m. We perform calculation on consecutive averaging windows over the entire time period of our dataset using M = 48 and δt = 50 s. We thus obtain estimation of frequency-dependent covariance matrices in 20-min-long windows that overlap with a 50% factor.

The covariance matrix is Hermitian and positive semidefinite, so it can be decomposed on the basis of its complex eigenvectors v associated with real positive eigenvalues λ

Then, eigenvalues are arranged in decreasing order, and the covariance matrix spectral width is computed as a function of frequency as

The time evolution of function σ(f) is shown in fig. S1A. Last, we average it between 0.5 and 5 Hz to obtain a simple 1D function characterizing the level of wavefield coherence across the network as shown in fig. S1B.

Temporal stability of tremors estimated from single-station intercomponent cross-correlations

If we consider a seismic source remaining at the same location with a constant mechanism, then the cross-correlation of the generated wavefield between two receivers or two components of the same receiver will remain stable in time, and it can be used as a “fingerprint” of this particular source (). The idea of this method is to compute correlation coefficient as a measure of similarity between intercomponent cross-correlation waveforms over time and to compute a mean value for a few consecutive windows to estimate the level of stability of these intercomponent cross-correlation. We filter the waveforms in a dominant tremor spectral band: 0.5 to 5 Hz and preprocess them in a similar way as for computation of the covariance matrix spectral width. Then, after computing intercomponent cross-correlations CC (i, j = E, N, Z) in moving windows of length Δt and shifted by Δt/2, we compute the correlation coefficients between cross-correlation waveforms from consecutive time windows and take the mean value of five consecutive computations. Using moving windows of length Δt = 200 s, the obtained functions provides us with an estimation of the stability of the cross-correlation waveforms over 600 s that, in turn, is related to the existence of a stable seismovolcanic tremor source. Using the three pair components at a single station, we therefore end up with three functions per station. We can then average these functions to obtain a single function per station ts(t) (n is the stations index). By averaging these functions over the whole network, we obtain a third network-based characteristic function ts(t) that we call as temporal tremor stability shown in fig. S1C.

Comparing values of ts at different stations can be useful in a region where different volcanoes can be active. Following this idea, we define the difference in temporal tremor stability as difference between the averages of ts for stations IR1, IR2, and SV6 close to Tolbachik volcano and for stations IR11, IR18, and SV13 close to Klyuchevskoy (Fig. 1A). We indicate with green points in Fig. 2B the sources for which this difference is greater than 0.05 that systematically correspond to locations close to Tolbachik, confirming the activation of the plumbing system beneath this volcano.

Network response function

This characteristic function is designed to be more sensitive to the presence of impulsive signals such as individual earthquakes and swarms. Therefore, for its computation, we apply a temporal normalization that does not fully suppress the amplitude information. Instead of applying the running absolute mean normalization, we divide the traces (once spectral whitening is applied) by their mean absolute deviation ().

We then compute the covariance matrix and its inverse Fourier transform that results in an ensemble of interstation temporal cross-correlations between all pairs of available stations (we use the dominant tremor spectral band: 0.5 to 5 Hz). Their smoothed envelopes is then computed by convolving the amplitude of the analytic signal derived from Hilbert transform with a 1.5-s-long Gaussian window, an optimal value suggested in ().

We assume that tremors are dominated by S-waves and compute the likelihood of their source positions with a grid-search approach (, ). For every tested source position in a 3D spatial grid, we compute travel times for all stations using a 1D S-wave velocity model slightly modified from () as shown in fig. S6.

Then, for every grid point r, cross-correlation envelops S are shifted by the time differences needed for the wave to travel from the tested source to the two stations i and j, and we lastly compute the spatial likelihood function R(r) by stacking at zero lag time the value of the shifted envelopes for all stations pairs aswhere N is the number of stations and t(r) is the travel time between the tested source and the station i. The 3D grid covers the region shown in Fig. 1 with a 650 m–by–650 m horizontal resolution and a 500-m vertical resolution down to 50 km depth.

We then evaluate the level of focusing of this spatial likelihood function to design another characteristic network-based function. For this, we first normalize R(r) by its sum over all grid pointsand then add a temporal normalization by the maximum over the total time period

The maximum value of is directly dependent on the level of focusing of the spatial likelihood function and is expected to be high for well-localized sources and low when a random incoherent wavefield is recorded. We therefore define the “network response function” NRF(t) as maximum of computed in every of 20-min-long 50% overlapping windows. Its evolution is shown in fig. S1D.

Source location

Hypocenters of LP sources are estimated as positions of the maximum of the spatial likelihood functions described in the previous section. For better isolating the dominant source, we apply a spatial filtering (, ). Namely, we compute a modified spatial likelihood function that is estimated in the dominating tremor frequency band (0.5 to 5 Hz) from cross-correlations obtained via an inverse Fourier transform of the filtered covariance matrix computed from the complex outer product of the first eigenvector v1(f) with itself

This filtering is aimed at denoising the spatial likelihood function to enhance the dominant sources. Examples of located sources are shown in fig. S7. Positions of their maxima give us estimations of source hypocenters for corresponding windows.

Signal detection and classification

On the basis of methods described above, for every 44,020 20-min analyzed window, we compute three spatial coordinates of the hypocenter and the values of the three network-based characteristic functions. We then use them simultaneously to select windows with “detected” signals and to approximately classify them into different groups such as LP tremors, individual earthquakes, and swarms.

First, we removed all windows for which the hypocenter was located on the limit of the 3D grid to reject possible detections associated with teleseismic and regional (nonvolcanic) earthquakes. On the basis of this criterion, from a total number of 44,020 windows, we retained 14,179.

Then, based on a visual inspection of the three computed network-based functions (figs. S2 and S4), we established respective threshold values separating “detections” from noise (windows with all three functions not passing the thresholds): 0.8 for the spectral width, 0.17 for the temporal tremor stability, and 0.45 for the NRF. On the basis of these thresholds, the number of detections was reduced to 13,709.

Figure S2 (B, D, and F) shows the 2D diagrams comparing pairs of network-based characteristic functions for the 13,709 retained detections. Similar to results of (), the largest concentration of points (more than 85%) is in the area with relatively small spectral width and small NRF that corresponds to tremors. Another “branch” with large NRF and large spectral width corresponds to earthquakes. The separation between these two branches is less clear than at the Piton de la Fournaise volcano (), which is due to the much more sustained and variable seismic activity in the KVG. We define lines (1), (2), (3), and (4) to classify the retained detections in three main groups: tremors, earthquakes, and “intermediate.” Line expressions are (1) ; (2) ; (3) tstremor = 1.5 × (NRF − 0.33) − 0.1; and (4) tsearthquake = − 0.03 × (NRF − 0.05) + 0.164.

A visual inspection of respective signals (fig. S3) shows that intermediate detections mostly correspond to swarms of LP earthquakes. Therefore, for our final analysis, we selected 12,384 “tremor” windows and 643 intermediate windows whose temporal evolution is shown in fig. S4. We also checked how excluding earthquake detections from the analysis could affect our results and inferences, as shown in figs. S4 and S10, and concluded that this criteria selection was not critical.

Spatial density of tremor sources

We stack for all selected time windows to obtain a spatial density function representing spatial distribution of seismovolcanic sources. To avoid contaminations from location artifacts, we focus on most probable source positions and we keep for stacking only points for which . The results of this procedure are shown in Fig. 1.

We also use the spatial likelihood function to visualize the uncertainties associated with tremor migration episodes shown in Fig. 2 (C and D). For every source, we keep the values of the likelihood exceeding 98%. Then, we smooth the retained values in time in 2-hour-long windows. Last, the source density function is renormalized between 0 and 1 in every time window to better highlight the migration patterns.

Covariance matrix parameters choice

Two main parameters controlling the computation of the covariance matrix are the number of subwindows M and their length δt. M = 48 was selected to slightly exceed the number of stations used in the computation (45 stations in our case) to ensure the covariance matrix estimation to be statistically robust.

The length δt should be long enough so that the wave has time to propagate through the entire seismic network. The largest interstation distance in our network is close to 100 km, and the slowest seismic velocity in the media is ∼2 km/s (). This implies that δt should not be significantly smaller than 50 s.

We then test a range of δt spanning from 40 to 70 s to find a parameter that optimizes detection of sources in the mid and lower crust (>15 km). Longer windows would prone detecting very long duration sources like shallow tremors rather than shorter and more impulsive events more present at depth. We show in fig. S11 the influence of the length of subwindow δt on the result of tremor source location in December 2015, where δt is progressively increased from 40 to 70 s with a 5-s step. It can be seen that with δt < 50 s, only a very few shallow sources are detected and that deep sources start to be missed with δt > 55 s.

As another test, we decided to compare our detections with some well-known sources of LP seismic radiation in the KVG, i.e., with the DLP earthquakes (). We selected 45 sufficiently strong (M ≥ 1.3) DLP earthquakes that occurred during the campaign period and were correctly located by the KBGS (we established a simple threshold of minimum of 25 stations used for the location). The selected events are listed in table S1. We show in fig. S12 the resulting percentage of detection for shallow and DLP sources during December 2015 and the percentage of DLP detection in table S1. The detection of deep sources clearly maximizes with δt = 50 s, and therefore, we chose this subwindow length for our analysis.

Uncertainties of source location

Uncertainties of the source location can be estimated from the modified spatial likelihood function computed from the dominant eigenvector of the covariance matrix (Eq. 7). For every detected source, we measure the size of the area where this likelihood exceeds 98% laterally and vertically to estimate its respective location errors. We then compute histograms of these errors in different depth ranges from all detected sources and obtain 2D distributions shown in fig. S8 and use their mean values at every depth to estimate average errors. As a result, we estimate that while both horizontal and vertical errors are about 4 km near the surface, the latter increase more rapidly with depth reaching its maximum value of ∼14 km beneath the crust-mantle boundary (the apparent decrease of this error beneath 40 km is an artifact of the limited depth range used in the location procedure). As an independent test, we compare hypocenters for 37 detected DLP earthquakes (table S1) obtained with our method and based on standard approach minimizing S and P travel time residuals (fig. S9). We represent the difference of location obtained from the classical manual method and from our automatic method as blue points in fig. S8 and find them well in agreement with the errors estimated from the spatial likelihood function.