Real-time and in-situ assessment of conduit permeability through diverse long-period tremors beneath Aso volcano, Japan

Abstract Long-period signal (LPs, 0.2–2 s) and very long-period signal (VLP, 2–100 s) observed in the shallow volcanic plumbing system are typically repetitive and time-invariant in their location and source mechanism, offering in-situ probes of hot fluid transport over the eruption cycles. While the amplitude and activity of volcanic-tectonic earthquakes and LP events have been commonly used to infer overpressure within their source region, one missing link is an observable that may permit inference on the change in the permeability of the conduit plug/wall, which can regulate the degree of pressurization, affect the mechanical strength of the surrounding rock, and consequently the likelihood of an upcoming eruption. Here we show that during the 2011–2016 eruption cycle at Aso volcano in Japan, long-period tremor events, a VLP of ~15 s period, with opposite waveform polarity can be systematically detected and categorized as pressurization and depressurization events in the same crack-like conduit. We suggest that, depending on the strength of the surrounding rock and the permeability of the crack-like conduit wall/plug, pressurization due to magmatic heat and vaporization is more likely to occur when a less permeable conduit plug/wall can effectively keep the gas inside the crack-like conduit. On the other hand, depressurization is prone to occur if the conduit wall/plug permeability is sufficiently high to allow gas to escape from the conduit. Considering the amplitude of LPT proportional to the conduit overpressure, contrasting energetics of these diverse LPT events allows us to define whether the conduit is prone to pressurization or depressurization, providing a framework to infer how the permeability of the conduit wall/plug may evolve over an eruption cycle.


Introduction
Understanding how a volcano works and foreseeing future eruptions are the ultimate goals in volcanology. Recent geophysical, geodetic, geochemical and petrological evidence has pointed to a modern picture of the trans-crustal magmatic system with a lens of crystalrich mush and crystal-poor melt where destabilization and remobilization of such a system lead to an upcoming eruption (e.g., Cashman et al., 2017;Sparks and Cashman, 2017). As an eruption occurs, when the fluid/gas pressure overcomes the strength conduit plug or/and wall, it is of great importance to make a direct assessment of how the conduit strength and pressure vary over eruption cycles, either in a period of apparent quiescence or during unrest before an impending eruption. Conduit pressure is dictated by accumulated gas pressure, which largely depends on magma flow rate, magma viscosity, volatile solubility, and permeability of magma, conduit wall and the formation of conduit plug (Heiken et al., 1988;Jaupart and Allègre, 1991;Woods and Koyaguchi, 1994;Jaupart, 1998;Gonnermann and Manga, 2003;Sparks, 2003;Diller et al., 2006;Burton et al., 2007;Edmonds and Herd, 2007;Okumura et al., 2013;Gonnermann et al., 2017;Heap et al., 2018a). On the one hand, magma degassing can result in bubble nucleation, growth, and coalescence, ultimately forming an interconnected network that governs magma permeability and dictates outgassing. On the other hand, the strength and permeability of a conduit wall/plug could vary against loading size, rate, duration, and temperature (e.g., Benson et al., 2012), crack/fracture network (or porosity), overpressure and heating (e.g., Lavallée et al., 2008Lavallée et al., , 2013Okumura et al., 2010;Plail et al., 2014;Gaunt et al., 2016;Heap and Wadsworth, 2016;Heap et al., 2017Heap et al., , 2018bKushnir et al., 2017;Bain et al., 2019).
Along with the conduit overpressure, magma and conduit wall/plug permeability seem to be one of the most critical parameters that can help evaluate whether the shallow plumbing system is prone to failure and eruption (e.g., Voight, 1988;Kilburn, 2003Kilburn, , 2018Sparks, 2003;Hammer and Neuberg, 2009). Not only it governs the ability to withhold gas inside the conduit (e.g., Edmonds and Herd, 2007;Collinson and Neuberg, 2012), it also dictates the mechanical properties (or strength) of the conduit wall/plug and regulates conduit pressure (e.g., Heap and Wadsworth, 2016). On the one hand, the permeability of conduit wall, the formation of conduit plug, and the conduit overpressure are intimately linked with magma ascent and outgassing efficiency (e.g., Diller et al., 2006). An increased conduit overpressure may be due to fast gas accumulation, low magma permeability, a stronger and more impermeable conduit wall/plug, resulting in pressurization. Similarly, a decreased conduit overpressure may be due to slow ascent of gassed magma, high magma permeability, or/and a weaker, more permeable conduit wall/plug, resulting in depressurization. On the other hand, high (low) pore pressure can effectively enhance (reduce) permeability (Zoback and Byerlee, 1975;Walsh, 1981;Berryman, 1992;Nara et al., 2011; see also Paterson and Wong, 2005) and such nonlinear feedback or two-way coupling between permeability and overpressure (or effective stress) is probably crucial, but not necessarily included in the simulation (e.g., Diller et al., 2006;Girona et al., 2019). Nevertheless, Understanding such critical parameters and their temporal evolution is needed to help realize physics-based eruption prediction models (e.g., Melnik and Sparks, 1999;Segall, 2013).
In a way, it is important to know how close the overpressure to the strength of the surrounding rock is. If the permeability is low (or the mechanical strength remains high), high overpressure does not necessarily result in the rock failure. Probably, it's a combination of high overpressure and high permeability/low mechanical strength that most likely leads to rock failure and an imminent eruption. While the pressure inside the conduit, to some extent, may be inferred from the activities of volcanic-tectonic earthquakes (VT), volcanic tremors, long-period signals (LP) (e.g., Minakami, 1974;Chouet et al., 1994;Chouet, 1996;McNutt, 1996;Roman and Cashman, 2006) and modeling of ground deformation (Poland et al., 2006;Segall, 2013;Pinel et al., 2014;Fernández et al., 2017;Magee et al., 2018), the permeability of conduit wall/plug and magma are transient in nature and it is nontrivial to infer changes in magma and conduit wall/plug permeability in situ over an eruption cycle. In particular, despite laboratory studies of limited sample size often underestimate the permeability of rock mass in the field, it has demonstrated that the permeability can be modulated by numerous factors, including stress loading cycle, transient variation in pore pressure (Heap et al., 2015;Farquharson et al., 2016), thermal perturbation (Heap et al., 2018b), chemical reaction , and hydrothermal alteration Mordensky et al., 2019). Nevertheless, recent efforts suggest that the amplitude of VLP may be linked to the amount of SO 2 emissions (e.g., Kazahaya et al., 2011;Nadeau et al., 2011;Waite et al., 2013;Zuccarello et al., 2013), which are potentially modulated by conduit permeability (e.g., Edmonds et al., 2003).
In the following sections, we will first review seismic activities observed near shallow volcanic plumbing systems and highlight signals that are potentially viable for monitoring conduit pressure, focusing on long-period signal (LP) and very long-period signal (VLP). Subsequently, we elaborate our effort in the analysis of VLP events at one of the most active volcanoes in Japan, Aso volcano. In particular, we characterize the variability and diversity of VLP events and construct a catalog documenting how such diverse VLP event families may be utilized to infer conduit permeability and pressure over the 2011-2016 eruption cycle.
While the diversity of these VLP events provides key information on processes inside the shallow magmatic or/and hydrothermal systems beneath volcanoes, there is generally a lack of systematic analysis examining how the variability or diversity of VLP in a single volcanic system changes over one or multiple eruption cycles. At Aso volcano, VLP waveforms with opposite polarity are understood as a result of depressurization and pressurization of the shallow conduit near sea level, and they are linked to outgassing and magmatic heating in the same source region (e.g., Kaneshima et al., 1996). Here we first briefly outline the volcanic activity at Aso volcano and the historical observation of VLP events.
Typically, surface volcanic activities follow a cycle starting with a period of quiescence and the presence of the crater lake. During the initial unrest, minor phreatic eruptions or/and mud-eruptions can occur. The disappearance of the crater lake is superseded by red glows at the crater bottom, followed by ash emission and ultimately Strombolian eruption, whereas phreatic, phreatomagmatic or/and explosive eruptions can also occur (e.g., Miyabuchi et al., 2006;Ikebe et al., 2008), often during the recovery of the crater lake or/and the end of the eruption cycle (e.g., Sudo et al., 2006). Presumably, the magmatic heat and input of hot volcanic fluids (e.g., Terada et al., 2012;Shinohara, 2013) originate from a deep conduit or/and a magma chamber at 6-7 km depth (e.g., Sudo and Kong, 2001;Sudo et al., 2006) and possibly even deeper below 10 km depth (e.g., Abe et al., 2010Abe et al., , 2017Unglert et al., 2011). The fact that the level of the crater lake fluctuates drastically is evidence of continuous and substantial thermal energy feeding from below.
In the case of the 2011-2016 eruption cycle, except minor phreatic eruptions in May 2011 and late 2013, Aso volcano was generally in quiescence. However, the temperature measured on the south wall of the crater shows a remarkable increase from about 60°C to about 240°C between mid-2011 to mid-2012, steadily rising to about 300°C by late 2012. The total magnetic intensity measured near the 1st crater indicates demagnetization between December 2010 and October 2012, also supporting a temperature rise underneath the crater (see also Tanaka, 1993). In the meantime, the crater lake level decreased from mid-2010 and dropped to a very low level in mid-2011 and mid-2012 before recovering in the next few months, respectively. The crater lake level again decreased to a very low level from March 2013 to late 2013. In general, the SO 2 emission remained mostly at the background level of 200-400 tons/day during 2011-2013.
In January-February 2014, a few small ash emissions occurred and the SO 2 emission increased from the background level to 1000-3000 tons/day in early 2014. The crater lake began to dry up in July 2014, and the temperature measured at the crater bottom increased rapidly from 300°C to about 600°C. At the end of August 2014, ash eruptions occurred, and intermittent Strombolian eruptions took place from the 25th November 2014 until the beginning of May 2015 Miyabuchi and Hara, 2019), when the southern part of the pyroclastic cone collapsed. Since then, the crater lake has recovered and remained at a very low level, and the SO 2 emission remained high throughout the rest of 2015 and most of 2016. Besides several minor phreatic/ash eruptions, a phreatomagmatic eruption with a small-scale pyroclastic density current occurred on the 15th September 2015 (Miyabuchi et al., 2018), and another phreatomagmatic explosion, preceded by elevated SO 2 emission (N10,000 tons/day), occurred on the 8th Oct 2016, where the plume height reached about 12 km above sea level (Ishii, 2018;Ishii et al., 2018;Sato et al., 2018). After the October 2016 eruption, the crater lake level gradually recovered, effectively marking the end of the 2011-2016 eruption cycle.  (Tanada et al., 2017), respectively. The star marks the location of the 1st crater, Naka-dake volcano. The rectangle in cyan marks the location of LPT source region, i.e. the crack-like conduit (e.g., Legrand et al., 2000;Yamamoto et al., 1999).

VLP observed in Aso volcano
VLP was first reported by Sassa in 1935(Sassa, 1935, and it was categorized as volcanic micro-tremors of the second kind, with a period of 5-8 s. More recently, Kawakatsu and co-workers (e.g., Kaneshima et al., 1996;Yamamoto et al., 1999;Kawakatsu et al., 2000;Legrand et al., 2000) analyzed broadband seismic data of several months in 1994 and systematically characterized VLP spectra and source properties. Since long-period tremor, or LPT, was first termed by Kawakatsu and coauthors to discuss VLP signal in Aso volcano and it has been widely recognized in the community, we use the same nomenclature hereafter for consistency.
As documented by Kawakatsu et al. (2000), LPT typically has a duration of 60 s and a dominant period of~15 s, with secondary spectral peaks at 7.5 s, 5 s, and 3 s, respectively. Moment tensor inversion indicates that LPT is located at~200 m southwest of the first crater near sea level. It corresponds to a low seismic velocity (Huang et al., 2018) and a modestly high conductivity zone (e.g., Hase et al., 2005;Hata et al., 2016Hata et al., , 2018Kanda et al., 2019), most likely a region of a hydrothermal reservoir. The source of LPT is predominantly isotropic (N90%) with a minor tensile crack component Legrand et al., 2000) and their excitation is generally considered as a result of fluid-solid interaction triggered by a transient pressure change in a crack-like conduit (e.g., Yamamoto et al., 1999;Kawakatsu et al., 2000;Kawakatsu and Yamamoto, 2007).
The crack-like conduit is often regarded as a pathway transporting hot gas toward a shallower hydrothermal system below the crater bottom (e.g., Hase et al., 2005;Takagi et al., 2006;Kawakatsu and Yamamoto, 2007;Kanda et al., 2008Kanda et al., , 2019Minami et al., 2018). Notably, high-frequency signals or/and LP-type events, at least in some instances, also occur in sync with LPT and fumarole activities (e.g., Sassa, 1935;Kikuchi, 1974;Churei, 1985), suggesting a causal fluid/pressure connectivity between the crack-like conduit at sea level and surface volcanic activities (e.g., Ichimura et al., 2018;Kawakatsu and Yamamoto, 2007;Mori et al., 2008;Takagi et al., 2006Takagi et al., , 2009. In general, the LPT waveform starts with a negative initial polarity and it was generally interpreted as a result of the rapid release of the liquid-gas mixture through fractures (Kaneshima et al., 1996;Kawakatsu et al., 2000). LPT waveform with a positive initial polarity was occasionally recognized during phreatic eruptions in 1993-1994, and it was typically regarded as a result of rapid pressurization due to vaporization induced by magmatic heat and supply of hot fluid/gas from below (Kaneshima et al., 1996;Kawakatsu et al., 2000).
As detailed above, while the early investigation of Kaneshima et al. (1996) and Kawakatsu et al. (2000) documented LPT waveforms of opposite polarity, it was limited to a quiet period of volcanic activity where occasional phreatic eruptions occurred in the presence of the crater lake. The diversity of LPT families during active periods, however, was unknown. In a way, it is unclear how the variability of LPT waveforms, polarity, and their activities can be understood concerning the entire eruption cycle of Aso volcano. Given the latest development of the fundamental volcano observation network, or V-net (Tanada et al., 2017), as well as the existing volcanic network from Japan Meteorological Agency (JMA), the 2011-2016 eruption cycle offers an ideal framework to systematically explore the diversity of LPT and their temporal evolution against the entire eruption cycle. In particular, not only does it allow an assessment of the state of the shallow plumbing system during low surface activities in the early stage of the eruption cycle, but it also presents an unprecedented opportunity to examine and explore how diverse LPT families and their respective activities vary from a period of volcanic quiescence to active episodes of Strombolian or/and phreatomagmatic eruptions.
In the following sections, we document the systematic detection and characterization of the diverse LPTs. After briefly outlining the data and seismic network in Section 2, we describe the methodology to detect and construct a catalog of diverse LPT families in Section 3. In particular, we report the identifications of several LPT families using the continuous wavelet transform (CWT) and the matched-filter (MF) analysis. In Sections 4 and 5, we highlight the temporal variation in the activities of these diverse LPT families and elaborate on how these new observations can be used to infer the temporal change of the state of crack-like conduit and surface volcanic activities. In particular, we contrast the energetics of these LPT families and suggest a simple metric to assess when the conduit is prone to pressurization or depressurization, allowing us to infer variations in conduit wall/plug permeability over the eruption cycle.

Data
Since late 2010, National Research Institute for Earth Science and Disaster Prevention (NIED) in Japan began establishing the fundamental volcano observational network, V-net, across active volcanoes (Tanada et al., 2017). By 2017, 16 active volcanoes are equipped with V-net. At Aso volcano, V-net includes four collocated surface broadband seismometers, borehole short-period seismometers, and tilt-meters, typically a few kilometers away from the 1st crater (Fig. 1). The stations ASIV and ASHV have been available since 2011, whereas the stations ASTV and ASNV have become available since 2016 (see also Fig. S1). While V-net typically includes collocated GNSS stations, these data are not yet made publicly available. Nevertheless, coupled with the five existing short-period seismometers from JMA Volcanic Seismometer Network, it constitutes an ideal seismic network to conduct a systematic analysis of the diversity of LPT and their activities over the 2011-2016 eruption cycle (Fig. 1). Broadband seismometers of V-net are equipped with a Nanometrics-240 sensor (e.g.,~250 s natural period), which is comparable to the bandwidth of the broadband seismic network F-net deployed by NIED across Japan. On the other hand, the natural period of short-period seismometers of JMA is about 1 s.
While data from broadband stations ASHV and ASIV are ideal to detect LPT of a 15 s dominant period, the capability of the short-period sensor in recovering long-period signals depends largely on the quality of the data and signal-to-noise ratio. Among the short-period stations, except ASO2, other short-period stations are noisy below 5 s (Fig. S2). Furthermore, we calibrate the sensor misorientation against longperiod P wave particle-motion, following the approach of Lim et al. (2018) and examining the amplitude of tangential receiver functions at zero time. Typically, the estimated sensor misorientation is less than ±3.0°(see also Supplementary material). After correcting sensor misorientation, we carefully examine the amplitudes of low-frequency modes (Davis et al., 2005) and long-period surface waves (Ekstrom et al., 2006) from large earthquakes. Among V-net broadband seismometers and nearby F-net seismometers, the amplitudes are consistent within ±3%, indicating the consistency of their sensor gains (see also Supplementary material).
As detailed below in Section 3, after considering data quality, data gap, and continuity, three high-quality vertical channels from ASHV, ASIV, and ASO2 are used for initial LPT detection. Once the initial detection is made, three-component waveforms of LPT are stacked to form 9channel waveform templates for the final detection. Since data from ASO2 are not available between December 2013 and April 2014 (Fig. S1), the detection capability is notably compromised in this period. In the following section, we detail the method of LPT detection and report the discovery of several LPT families and the pattern of their activities.

Method and detection of diverse LPT families
To conduct a systematic analysis of the diversity of LPTs at Aso, our approach involves three major steps and the workflow is presented in Fig. 2. First, we conduct a visual inspection and construct a robust LPT waveform template. Secondly, we perform an automatic detection in the continuous wavelet transform (CWT) domain, capable of identifying diverse LPTs. These initial LPT detections are evaluated, classified, and subsequently stacked to form LPT waveform templates. Thirdly, with the LPT waveform templates, we apply the matched filter (MF) against continuous-waveform data to construct the catalog of such diverse LPT families during 2011-2016. To appreciate the performance of the MF method and evaluate potential false/missed LPT identification, we perform the MF against a synthetic LPT catalog and examine the false/ missed pick rate against signal-to-noise ratio (SNR) and cross-correlation (CC) coefficient thresholds, which provides an objective selection criterion to assess the robustness of LPT detections and their activities. Hereafter, we document the data processing procedures.

Initial LPT template
To build a robust LPT template for automatic detection, we first perform the CWT to inspect LPT activities in Oct 2014, during which surface volcanic activity is prominent. Compared with the short-time Fourier transform, CWT does not require a pre-defined window size and can perform local analysis on a non-stationary signal (Grossmann and Morlet, 1984;Mallat, 2009). The continuous time series is first sliced into ten-minute long windows and the instrument response removed. Velocity seismograms in each slice are simultaneously examined in the time domain, frequency domain, and CWT time-frequency domain in the frequency band of 0.02-2 Hz (e.g., Fig. 3). After visually inspecting over 9000 windows, we recognize hundreds of isolated wave packets with strong amplitudes in the vertical and radial components, oscillating b60 s in the time domain and displaying two dominant spectra peaks at the periods of about 15 s and 7.5 s, all diagnostic features of LPTs. The waveform amplitude also decreases from ASO2 to ASIV and ASHV, consistent with an LPT source located close to the 1st quarter. In the subsequent CWT automatic detection scheme, we use the waveforms in the vertical component of stations ASO2, ASHV, and ASIV as the initial LPT waveform template.

CWT automatic identification of LPT
Having constructed the LPT waveform template, an automatic detection scheme is carried out by CWT and cross-correlation method against the 6-year continuous data. The waveform template and a continuous data window of 2 min are first decimated to 10 Hz, with the instrument response removed and normalized before calculating the CWT spectrograms and the time-frequency-amplitude map (TFAM). We compute a 2D cross-correlation between template TFAM and data TFAM. If the average cross-correlation coefficient is over 0.8, we mark the identification of LPT as well as its timing (Fig. 4). Data windows are shifted 12 s and the same process is repeated. In our approach, it is important to ensure that the detection scheme is robust but flexible enough to identify LPTs with diverse waveform shapes, opposite polarity, or/and variable spectra peaks. Therefore, the phase information is not utilized at this stage. Furthermore, we use fifty frequency points in the band of 0.02-2 Hz to ensure the resolution of the spectrogram and the constraint for detection.
Several stringent criteria are used to remove possible false detections and ensure high-quality detection. First, if the average amplitude ratio between the unfiltered seismogram and the filtered seismogram (10-30 s) of the three different channels is over 13, the detection is removed (Rejection criterion I). While the threshold of amplitude ratio is admittedly somewhat arbitrary, it warrants strong, long-period energy and quality LPT detection. Secondly, to remove duplicate LPT detection within a 2-min data window, we keep the detection only if 95% or more of the waveform energy concentrates in the central 30 s (Rejection criterion II). Thirdly, since LPT is located about 200 m southwest of the 1st crater, we expect the amplitude decays quickly away from ASO2 in the near field. Therefore, we retain LPT detections only if the maximum amplitude at ASO2 is at least 1.5 times greater than the maximum amplitude at ASHV and ASIV (Rejection criterion III). This eliminates false identifications where the source of energy originates far away from the 1st crater or the known LPT source location (e.g., Kawakatsu et al., 2000;Legrand et al., 2000).
Out of about 87,000 initial detections, 90% of the initial detections fail to pass these three rejection criteria and are removed. Finally, we perform an additional visual inspection over the remaining detection and remove several false detections caused by the calibration pulses from the short-period channel ASO2. The final number of LPTs detected by the automatic CWT scheme is 7455 and the detection threshold at the station ASHV is about 0.1 μm/s.

Identification of diverse LPT families
Following the CWT detection of LPTs, we systematically examine the stability/variety of LPT waveform shapes in a 30-day moving window with an overlap of 15-days. In many instances, LPT waveforms within a 30-day time window are quite similar and there are up to hundreds of LPTs during the active period between late 2014 and late 2015. On the other hand, we can observe that the waveforms of LPT-A and LPT-B are nearly symmetric but anti-correlated, whereas the waveforms of LPT-C and LPT-D are asymmetric, suggesting an energy interference between the main spectral peak near 15 s and the secondary peak at about 7.5 s are not uniform among these LPT families (Fig. 5e-h). In the next section, three-component waveform stacks of these 4 LPT families from ASO2, ASHV, and ASIV will be fed into the MF scheme to detect their activities in the 2011-2016 eruption cycle.

Detection of diverse LPT families in 2011-2016 with the MF
As demonstrated in previous sections, our CWT-based automatic detection is effective in picking quality LPT with high SNR. However, the detection threshold is high, and we potentially miss many weak LPTs. A more sensitive detection scheme, the multi-channel matched filter (MF) is a powerful technique that can help extract weak but coherent signals embedded in the noise using pre-defined templates (Turin, 1960), such as low-frequency earthquakes in subduction zones (e.g., Shelly et al., 2007) and LPT in Aso volcano (Ikeda, 2005). Here we adapt the MF scheme to improve the detection capability against the aforementioned diverse LPT families.
Specifically, in a period band of 10-30 s, we cross-correlate the data against each LPT template in a 2-minute moving window and compute the averaged CC function across 9 channels (Fig. 6). We disregard CC functions with negative peaks and a detection/classification is made against the remaining CC functions with the highest peak amplitude. Concerning a typical LPT waveform of 60 s, if there are multiple identifications within a 30-second window, the detection with the highest CC peak is retained, effectively removing redundant picks. Finally, we apply the Rejection criterion III previously introduced in the CWT automatic detection scheme to remove detections that are not consistent with the source location of LPT.
Before discussing the result of MF-based detection of LPTs in 2011-2016, it is useful to examine the capability of the detection scheme in detecting and classifying diverse LPT families. In particular, the construction of an MF-based catalog depends on the CC threshold and the minimum time separation of neighboring events. In particular, if the amplitude cutoff of the CC is too low, numerous false detections may occur. If two neighboring events occur too close in time, we may potentially miss one of the two events. Therefore, we construct a synthetic catalog by using four families of LPT waveform templates to examine the false-pick rate against SNR and CC thresholds in these two scenarios. As shown in Fig. S4, the MF worked reasonably well even when SNR and CC are low. With CC ≥ 0.45 and SNR ≥ 1.8, the percentage of correct picks can reach 99% and the percentage of false picks or misidentifications is near 1%. With CC ≥ 0.63 and SNR ≥ 2.92, the correct picks reach above 99.9%, the misidentifications are about 0.1%. Our MF scheme detects 2,740,745 events during 2011-2016, and with the criteria of CC ≥ 0.63 and SNR ≥ 2.92 and a nominal false-pick rate of 0.1%, we retain 208,518 events in the final catalog. This ensures the robustness of LPT event classification, allowing us to elaborate and contrast the activities of different LPT families over the eruption cycle.

Location of LPT families and temporal change?
To understand how the activities of LPT families may vary against surface volcanic activities, we first produce monthly waveform stacks of these LPT families at the closest station ASO2 and examine the temporal variations of their particle motions. As shown in Fig. 7, not only are the polarizations of these LPT families consistent with each other, but they also display relatively minor differences during 2011-2016, indicating a limited variation in depth or location by~100 m. It appears that the depth of LPT during the quiescent period of 2011-2013 is slightly deeper than that during the unrest and active period of 2014-2016. On the other hand, the particle motion of these LPT families is well reproduced by the source location and mechanism reported by Kawakatsu et al. (2000) and Legrand et al. (2000). The consistency also indicates that the source locations of these LPT families are not only in very close proximity during the 2011-2016 eruption cycle, but they also remain stable and consistent with the LPT source properties defined during a relatively quiet episode of low surface volcanic activity in the 1990s (e.g., Kawakatsu et al., 2000). While Hendriyana and Tsuji (2019) reported that LPT appears to migrate significantly after the 2016 Kumamato earthquake, steady polarization observed across the entire 2011-2016 eruption cycle does not support such a finding.

Activities of diverse LPT families during 2011-2016 eruption cycle
To explore how the activities of LPT families may correspond to surface volcanic activities, we first examine their median amplitudes and weekly numbers (Fig. 8). The weekly amplitude of these LPT families is typically weak at~0.2 μm/s during the period of relative quiescence between 2011 and early 2014. By July 2014, their amplitude increases rapidly, peaking at~5 μm/s in early 2015, when intermittent Strombolian eruptions occurred. However, except for the short episodes of two notable phreatomagmatic eruptions in September 2015 and October 2016, the amplitude of these LPT families after the Strombolian eruption activity is generally steady at about 0.5 μm/s, about a factor of 2 of the pre-eruptive level before mid-2014. Most notably, the amplitudes of these LPT families are comparable and they generally vary consistently during the 2011-2016 eruption cycle (Fig. 8b).
In contrast, the weekly numbers of these LPT families do show a considerable difference, and these observations do not depend on the choice of false-pick rate threshold (Figs. 8, S5). In particular, during the period of relative quiescence between 2011 and early 2014, the level of activities of LPT-A and LPT-B is generally comparable, with minor activities between mid-2012 and early 2013 (Figs. 8, see also S5). However, the activities of LPT-C and LPT-D are very low during this period. Between mid-2014 and the end of 2016, we identify several important features distinguishing the activities of these LPT families. Typically, LPT-A consistently becomes very active 1-3 months before the magmatic eruptions (e.g., November 2014 Strombolian eruption, September 2015 phreatomagmatic eruption, and 2016 phreatomagmatic explosion). Secondly, days to weeks before magmatic eruptions and other minor eruptions, LPT-B is typically quite active, and these eruptions frequently occur as the activities of LPT-B are in decline. Thirdly, the activities of LPT-C and D mostly begin after mid-2014, and their activities also peak days to weeks before several minor eruptions in 2015 and 2016. Overall, the level of activity between LPT-A and LPT-B/C/D over the 2011-2016 eruption cycle is about 1:2 (Fig. S6).

Classify the activities of LPT families against surface volcanic activities in the 2011-2016 eruption cycle
Since the noise level and the detection threshold at any given time are the same among all LPT families, the systematic differences in their activities discussed in Section 4.2 probably manifest genuine differences in the state of the shallow crack-like conduit, dictating surface volcanic activities over the 2011-2016 eruption cycle. Following the discussions above, we elaborate on the activities of these diverse LPT families concerning surface volcanic activities and classify the activities of LPT families into three phases. To highlight the episodes of elevated LPT activities compared with the background, we define the time windows when the weekly number of LPTs exceeds the average weekly number over an 8-week long time window, indicating the episodes where LPT-A or LPT-B/C/D activities are particularly elevated.
As shown in Fig. 8, Phase 0 is marked by the generally weak and low LPT activities associated with low surface volcanic activities during the period of relative quiescence between 2011 and mid-2014. Phase 1 marks the episodes of elevated LPT-A activities weeks-to-months before the major eruptions involving juvenile magmas in November 2014 (Strombolian), September 2015 (phreatomagmatic), and October 2016 (phreatomagmatic). Phase 2a proceeds after phase 1 and marks the episodes of elevated LPT-B/C/D activities of high amplitude before the major eruptions. On the other hand, phase 2b marks the elevated LPT-B/C/D activities of low amplitude, typically followed by minor, non-magmatic (e.g., phreatic, mud eruptions) in 2016. Notably, phase 1, phase 2a, or phase 2b are not necessarily a single, continuous phase and they can include multiple time windows (or episodes) before each eruption.
While the definition of elevated LPT activities is not unique, our classification scheme provides a useful way to coherently relate the activities of diverse LPT families to surface volcanic activities. In the following discussions, we summarize our observational findings and provide plausible mechanisms that link the activities of diverse LPT families near sea level to surface volcanic activities in Aso volcano. In particular, we use the observed attributes of diverse LPT families to discuss the inferences on the evolution of overpressure and wall/plug rheology in the crack-like conduit during the 2011-2016 eruption cycle.

Diverse LPT families: magmatic heating vs. outgassing; pressurization vs. depressurization
Here we document new findings of the four LPT families which are closely located during the 2011-2016 eruption cycle. As discussed in previous studies (e.g., Kaneshima et al., 1996;Kawakatsu et al., 2000), moment tensor inversion of LPT waveforms highlights a source of predominantly volume (or pressure) change in the crack-like conduit. It not only buffers the upward heat transport from a deep-seated magma chamber but also leaks gas or liquid up toward the crater lake or/and pathways toward surface fumarolic activities (e.g., Mori et al., 2008). Intermittent supply of hot fluids from the deep conduit or magma chamber can result in vaporization and a pressure increase in the crack-like conduit, producing LPT-A with a positive initial polarity. A transient upward leakage of gas out of the hydrothermal reservoir and a sudden pressure drop may drive the crack-like conduit to resonate, resulting in LPT B/C/D families with a negative initial polarity. Therefore, we regard the LPT-A family as a proxy of the transient pressurization episode and LPT-B/C/D families as a proxy of depressurization episode and outgassing in the crack-like conduit. Fig. 6. An example of waveform cross-correlation and LPT detection using the matched filter (MF). (a) Waveform cross-correlations between the 10-min data waveforms (black) and LPT waveform templates from two broadband stations (N.ASIV.BU and N. ASHV.BU) and a short-period station (V.ASO2.U). Cross-correlation functions between the data and the four LPT waveform templates are presented in the lower panel. The waveform templates of LPT-A (red) show the highest cross-correlation against data waveforms at~340 s (cross-correlation coefficient, CC, of~0.99 and signal-to-noise ratio, SNR, of~44). (b) is the same as (a), but with an LPT detection of lower SNR (~6.20). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Amplitude of diverse LPT family and changes in conduit overpressure
As pointed out earlier, diverse LPT families display a notable and systematic change concerning surface volcanic activities. However, the amplitudes of these LPT families generally fluctuate rather consistently with a coherent trend (Fig. 8b), independent of the type of LPT family. Considering the geometry of the crack-like conduit remains largely constant, the amplitude of the tremors can be related to the overpressure (e.g., Aki et al., 1977;Chouet, 1985). As LPT source is predominantly associated with a volume and pressure change in the crack-like conduit , we suggest that LPT amplitude is dictated by a characteristic volume of hydrothermal fluid readily vaporized or outgassed inside the crack-like conduit, proportional to the change in overpressure.

Activities of diverse LPT family and changes in conduit permeability/ strength
One important aspect yet to address is the permeability of the conduit wall/plug, which presumably plays a major role in regulating the likelihood of pressurization and depressurization episodes. As demonstrated in cyclic loading/unloading experiments (e.g., Eberhardt et al., 1999;Vinciguerra et al., 2005;Heap et al., 2010;Benson et al., 2012), repeated pressurization/heating and depressurization/cooling associated with these LPT families likely extend pre-existing cracks or/ and promote the nucleation of new cracks, potentially consistent with more frequent LPT activities. Not only can these extended crack networks and new fractures increase the porosity and permeability of the wall rock and the conduit plug (e.g., Farquharson et al., 2016;Heap and Kennedy, 2016;Heap et al., 2017), but they also weaken their overall strength. Furthermore, repeated heating and cooling of rising gassed magma and cooling from the inflow of groundwater are also likely to occur in the well-documented hydrothermal system beneath the Naka-dake crater, introducing a higher degree of cracking (e.g., Burlini et al., 2007;Browning et al., 2016) and potentially leading to a reduction in the conduit wall/plug strength and stiffness Mordensky et al., 2019).
To assess whether the conduit is prone to pressurization or depressurization and to infer the changes in conduit wall/plug permeability and their strength, we propose a simple attribute by contrasting the energetics of pressurization and depressurization events. We define the weekly moment as the product of the weekly number and median velocity to compute the weekly moment ratio (WMR) between LPT-A Fig. 7. LPT waveform polarization over calendar time at ASO2. To measure LPT polarization of the four LPT families, the waveforms of each LPT family are normalized and stacked in a 30day moving window, with an overlap of 15 days. Synthetic waveform polarization (black) is calculated with the F-K method (Zhu and Rivera, 2002) using the source parameters and the half-space velocity model of Legrand et al. (2000). The observed waveform polarization of each LPT family is steady over 2011-2016 and they are consistent with the synthetics. Note that the polarization of waveform stacks with b50 events is considered less reliable and not included in the plot. and LPT-B/C/D. Assuming that magma injection rate, conduit temperature, and overpressure are constant at a given week, we use the WMR as a proxy to infer changes in the permeability of the conduit wall/ plug. Since the WMR does not depend on LPT amplitude and the level of background noise, it is not biased by detection capability, providing a robust and uniform measure against episodes of LPT activities of different level (see also Fig. S7).
Typically, the WMR varies over 3 orders of magnitude (~0.005-20) (Fig. 8c), with a high (low) WMR corresponding to the episodes prone to pressurization (depressurization). During phase 0 and phase 1 between 2011 and late 2014 (Fig. 8c), it stays mostly within 1-5, except the period of mid-2012 to mid-2013 where the WMR is relatively low at 0.1-1. The WMR quickly decreases from~2 to about 0.1 during phase 2a in late 2014 and typically varies between 0.1 and 0.5 during intermittent strombolian eruptions in early 2015 (Fig. 8c). The end of the Strombolian eruption in early 2015 marks the beginning of a sharp increase in WMR and it remains at a high level of~2 during phase 1 in mid-2015. The WMR again decreases sharply from 2 to about 0.1 or lower in phase 2a in mid-late 2015, continuously fluctuating at a low level of~0.1-0.3 during phase 2b until mid-2016. In late 2016, the WMR again increases to a relatively high level of about 1 in phase 1, followed by a lower WMR in phase 2b before the 2016 phreatomagmatic explosion.
We suggest that, depending on the strength of the surrounding rock and the permeability of the crack-like conduit wall/plug, pressurization (e.g., LPT-A family) due to magmatic heat and vaporization is more likely to occur when a strong or/and less permeable conduit plug/wall can effectively keep the gas inside the crack-like conduit (e.g., Heap  (Fig. 26, http://www.data.jma.go.jp/svd/vois/data/tokyo/STOCK/monthly_v-act_doc/fukuoka/2016y/503_16y.pdf). (e) compares the outgassing potential (Amplitude × 1/WMR) and campaign SO 2 emission flux from JMA (grey) (https://www.data.jma.go.jp/svd/vois/data/fukuoka/rovdm/Asosan_ rovdm/gas/gas.html). Black triangles and hatched regions mark the onsets and intervals of the November 2014 Strombolian eruption, the September 2015, and the October 2016 phreatomagmatic eruptions, respectively. Grey arrows show minor eruptions (phreatic, ash or mud eruptions) reported by JMA (http://www.data.jma.go.jp/svd/vois/data/tokyo/ STOCK/monthly_v-act_doc/monthly_vact_vol.php?id=503) and the Global Volcanism Program (2011,2012,2015,2017). and Wadsworth, 2016). On the other hand, depressurization (e.g., LPT-B/C/D families) is prone to occur if the conduit wall/plug permeability is sufficiently high to allow gas to escape from the conduit. In the meantime, if the overpressure, at least locally, overcomes the strength of the caprock or/and the conduit wall, additional cracks can form, and the gas can escape toward the crater bottom of the edifice. Nevertheless, by contrasting the level of activity between LPT-A family and LPT-B/C/D families, or equivalently the tendency of pressurization events against the depressurization events, we may infer the changes in the balance between the overpressure and the permeability/strength of the conduit wall/plug.

A new observational framework and inferences of overpressure and conduit wall permeability/strength
While experimental data and theoretical considerations have shown that pore pressure can effectively modulate permeability (Zoback and Byerlee, 1975;Walsh, 1981;Berryman, 1992;Paterson and Wong, 2005;Nara et al., 2011), nonlinear feedback or two-way coupling between permeability and overpressure is often not included in the numerical simulation (e.g., Diller et al., 2006;Girona et al., 2019). Therefore, it is the motivation of this paper to present relevant observational attributes that can be used to directly infer overpressure (LPT amplitude), permeability (WMR), and how they may evolve at different stages of volcanic activities. As discussed in the following, the activities of diverse LPT and a simple attribute of WMR offer a new framework to evaluate the interplay and feedback among conduit overpressure and conduit plug/wall permeability and rheology (Fig. 9).
With a low LPT amplitude and a weak activity (e.g., phase 0), we suggest the strength of the conduit wall/plug is generally strong and less permeable when the WMR is high, making the gas less likely to escape at this stage. The conduit overpressure is relatively low, subjected to steady but episodic magmatic heating. On the other hand, with a high LPT amplitude, a strong activity, and a high WMR (e.g., phase 1), the conduit wall/plug remains sufficiently strong, potentially subjected to elevated magmatic heating and overpressure. With a high LPT amplitude, a strong activity, but a low WMR (e.g., phase 2a), the strength of the conduit wall/plug is significantly reduced due to more extended cracks, allowing the gas to escape more effectively and frequently. Presumably, magmatic injection or heating from the deep conduit or/and chamber remains at a high level and the conduit is still under a relatively high overpressure. We consider that this is the stage where the overpressure is high and the mechanical strength of the conduit plug/wall is low due to enhanced permeability (stress-dependent), leading to the failure and an imminent eruption. Finally, with a relatively low LPT amplitude, an elevated activity, and a low WMR (e.g., phase 2b), the conduit wall/plug with extended fractures remains weak and more permeable as the conduit overpressure decreases, possibly subjected to relatively weak magmatic heating.
As summarized in Fig. 9, we suggest that the 2011-2016 Aso eruption cycle likely begins with a strong or/and less permeable conduit plug/wall under a relatively low overpressure. The Strombolian eruption in late 2014 and the subsequent phreatomagmatic eruptions in late 2015 and late 2016 are preceded by an early episode of a comparatively strong or/and less permeable conduit wall/plug and an increasing overpressure (prone to pressurization). After an increase of overpressure over days/months, the conduit wall/plug subjected to a high overpressure becomes weaker or/and more permeable, allowing efficient gas to escape (prone to depressurization) and leading to the upcoming major eruptions. Several minor eruptions in the late stage of the eruption cycle in 2016 are also preceded by the episodes of weak or/and more permeable conduit wall/plug, responding to the episodic depressurization and pressurization under relatively low overpressure (Figs. 8, 9). 5.5. Diverse LPT activities vs. surface volcanic activities, eruption style, and SO 2 emission?
We find that the WMR also corresponds to surface volcanic activities systematically, indicative of whether the upcoming eruption is of the magmatic or hydrothermal origin (dominated by gas). As discussed earlier in Section 4.3 and 5, the extended period of phase 1 or high WMR (~1) are typically the prologues of magmatic eruptions in 2014, 2015, and 2016, which involve juvenile magma. It is consistent with magma rising from the crack-like conduit toward the crater bottom, evidenced by the episode of prominent temperature increase (Fig. 8d). On the other hand, phase 2a typically corresponds to the stages of rapid decline in the WMR (~0.01-0.1), suggesting a more substantial gas escape from the crack-like conduit toward the crater bottom, resulting in a rapid buildup of gas near the crater bottom before the imminent eruption.
As discussed in the above section, the amplitude of VLP or tremor in other volcanoes has been correlated with SO 2 emission flux (e.g., Kazahaya et al., 2011;Nadeau et al., 2011;Waite et al., 2013;Zuccarello et al., 2013). However, among other geometric factors (e.g., cross-section area, characteristic length scale), SO 2 emission is presumably modulated by pressure drop as well as permeability (e.g., Edmonds et al., 2003). In most instances, syn-eruptive VLPs or tremors only allow such a comparison during eruptions, whereas repetitive LPT in Aso volcano provides an excellent opportunity to systematically compare SO 2 emission flux against LPT activity throughout the entire eruption cycle in 2011-2016. Here we define the outgassing potential as the product of LPT amplitude (a proxy of pressure drop) and WMR (a proxy of permeability) and compare it against campaign SO 2 emission data (JMA) in the 2011-2016 eruption cycle (Figs. 8e, S5e). We find that the outgassing potential correlates satisfactorily well with the SO 2 emission data, especially during the unrest and active period between late 2013 and 2016. This result also implies that geometric factors remain either approximately constant or they are inherently included in WMR or permeability, which is geometry dependent. It is conceivable that joint modeling of outgassing potential, high-frequency tremors (e.g., Ichimura et al., 2018;Takagi et al., 2006Takagi et al., , 2009 and LP (e.g., Mori et al., 2008) near the upper end of the crack-like conduit, as well as continuous SO 2 emission data can shed light on the evolution of the shallow conduit system and the eruption potential in a persistent degassing volcano such as Aso.

Outlook and future work
In this report, we have detected diverse LPT families during the 2011-2016 eruption cycle at Aso volcano in Japan. Systematics of their activities, as discussed above, offer a new perspective and insight into the interplay and feedback among key components regulating the style and timing of upcoming eruptions. The hydrothermal system beneath the Aso crater has been fruitfully discussed in the literature (e.g., Tanaka, 1993;Hase et al., 2005;Terada et al., 2012;Shinohara, 2013;Hata et al., 2016Hata et al., , 2018Minami et al., 2018;Kanda et al., 2019). If the polarities of these LPT families can be associated with pressurization and depressurization (e.g., Kaneshima et al., 1996;Kawakatsu et al., 2000;Legrand et al., 2000), measurements of the amplitude, frequency of LPT as well as the WMR in real-time provide very useful attributes to monitor the evolution of conduit permeability/strength and pressure in situ. In particular, as the physical properties of the fluid (e.g., gas, vapor, magma, etc.) or/and the elasticity of the conduit plug/wall will probably change the attenuation, excitation period and modal interference of LPT (e.g., Chouet, 2000, 2001;Chouet et al., 2003;Namiki et al., 2018). Systematically analyzing the spectral characteristics of these LPT families over time potentially provides a further constraint on the temporal evolution of fluid properties and conduit plug/wall rheology over the eruption cycle (Niu and Song, in preparation).
One of the outstanding questions we did not address is the triggering mechanism of these diverse LPT families. While we have associated diverse LPT families with outgassing, the supply of gassed magma, and vaporization, it is not entirely clear if the outgassing process or upward transport of gassed magma and vaporization directly trigger LPT. It is conceivable that these processes may simply respond to an unknown source(s). As discussed earlier, the amplitude of LPT does not depend on the type of LPT family, which could be consistent with the latter scenario. Furthermore, it is not clear if the triggering source is necessarily collocated or spatially separated from the LPT source. Is the LPT activity related to the occurrence of deep low-frequency earthquakes at advocated in some volcanoes (e.g., Shapiro et al., 2017), or does the triggering manifest as a seismic source that is currently unknown? These outstanding issues will be explored in future work. Nevertheless, detailed waveform analysis may be performed against historical eruptions in Aso and similar efforts can potentially be extended to other volcanoes where VLP (or LP) is identified. Temporal seismic velocity changes may be measured to infer variations in the strength of conduit wall rock (e.g., Nara et al., 2011). When near-field seismic observations are not available, far-field observations of Rayleigh waves can potentially help provide continuous monitoring (e.g., Hashida, 1990;Kawakatsu et al., 1994;Sandanbata et al., 2015).

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Fig. 9. A conceptual model of diverse LPT activity against conduit overpressure and rheology. (a) sketches the shallow plumbing system and how it evolves during the eruption cycle. The shallow hydrothermal reservoir including cap rocks, fluid path, and altered rocks near the crater bottom (e.g., Kanda et al., 2008;Minami et al., 2018). The crack-like conduit (LPT source region, Yamamoto et al., 1999), surrounded by a deep aquifer, connects the crater bottom toward the deeper conduit or/and the magma chamber Kawakatsu and Yamamoto, 2007). Red and blue arrows indicate pressurization LPT and depressurization LPT events, respectively, and the size of the arrows indicates the level of LPT activities or energetics. (b) The conceptual model links the observed attributes (e.g., LPT activity, amplitude, and WMR) to the variations in conduit wall/plug strength, permeability, and conduit overpressure. Strong and less permeable conduit plug/wall and low overpressure are manifested by low LPT activity, low amplitude, and high WMR (e.g., phase 0). Increases in micro-cracks in the conduit plug/wall, elevated LPT activity, amplitude, and high WMR indicate the conduit plug/wall remains comparably strong and less permeable under an elevated overpressure (e.g., phase 1). Low WMR, high LPT activity and amplitude point to a weak and more permeable conduit plug/wall of more extended cracks under a high overpressure (e.g., phase 2a), whereas low WMR, low amplitude, and high LPT activity reflect a weak conduit plug/wall under a low overpressure (e.g., phase 2b).

Acknowledgments
J. Niu acknowledges the support by Prof. Caijun Xu and the China Scholarship Council to make the research visit at UCL possible. We gratefully thank NIED and JMA providing us the high-quality waveform data of V-net. We are also grateful to Prof. Mare Yamamoto and Prof. Takuhiro Ohkura for hosting J. Niu in Aso Volcanic Observatory and Tohoku University during the early stage of this work, which help improve our understanding of LPT and Aso volcanic system. The uses of the UCL Legion High-Performance Computing Facility (Legion@UCL), the UCL Grace High-Performance Computing Facility (Grace@UCL), the UCL Emerald High-Performance Computing Facility (Emerald@UCL) and associated support services are acknowledged. Data processing and production of figures are implemented in Python with relevant modules such as ObsPy and PyCUDA. J. Niu and T.-R. A. Song are supported by the Natural Environment Research Council, UK (NE/ P001378/1 & NE/T001372/1). We thank the editor Diana Roman for handling the manuscript and two anonymous reviewers for their comments, which help clarify and improve this paper.