Hydrothermal signature on episodic deflation/inflation ground tilt at Aso Volcano

Ground deformation in volcanic areas induced by geothermal fluid circulation can reveal useful information about the dynamical processes occurring in the subsurface hydrothermal system. In the present work, we investigate tiltmeter time-series recorded at Aso Volcano during 2011–2016, a time interval during which different phases of volcanic activity occurred. We performed polarization analysis of the data and identified peculiar long-lasting (hours) transients, defined as Very-Long-period Tilt Pulses. The transients were further characterized in terms of waveform cross-correlation, particle tilt pattern, energy, and time distributions. The analyses indicate that such signals, which appear like deflation–inflation (DI) events, are associated with a Poissonian process whose underlying dynamics evolves over time always driven by a Poissonian mechanism. The obtained results have been interpreted in light of the available geophysical, geochemical and volcanological information. In this framework, the Very-Long-period Tilt Pulses may be ascribed to the depressurization/pressurization of the shallow hydrothermal system according to a fault-valve mechanism, which was active with different efficiency throughout eruptive and inter-eruptive phases.


Introduction
Aso Volcano (southwest Japan) comprises the Aso caldera which formed between 270 and 90 ka, and post-caldera central cones (Fig. 1) (see, e.g., Ohkura et al. 2019 and references therein).Nakadake volcano, which is the only active central cone in the last 10,000 years, is a composite of seven craterlets aligned North-South (NS).Since the 1933 eruption, only the northernmost crater (the 1st crater) has been active, following cycles from quiescence to eruption (Sudo et al. 2006;Miyabuchi et al. 2006;Ikebe et al. 2008).
The quiescence is often indicated by the existence of a hyperacid crater lake (pH in the range [− 1:1]), characterized by hot temperatures (40 to 80 °C) (see, e.g., Shinohara et al. 2018a).At the beginning of the unrest phase, minor phreatic or/and mud eruptions can occur.The crater lake drainage is accompanied by ash emissions and finally by Strombolian or phreatomagmatic eruptions.Phreatomagmatic and explosive activity can also occur during the refill of the crater lake (Miyabuchi et al. 2006;Ikebe et al. 2008).The cycle may take 5-10 years, although the actual volcanic activities sometimes deviate from this cycle.The crater lake alternatively dries and refills: most of the water in the lake is groundwater/meteoric water, but part of the water level comes from magmatic or, more in general, hot fluids, rested at depths of 6-7 km and even deeper up to 10 km (see, e.g., Sudo and Kong 2001;Sudo et al. 2006;Abe et al. 2010).The cyclic activity in the crater lake, located at about 1200 m asl, is controlled by the dynamics of the shallow hydrothermal system, which is a region characterized by low resistivity below the crater, whose top is at an elevation of 800 m (Kanda et al. 2019;Matsushima et al. 2020) extending up to 2 km bsl.This conductive region is highly fractured and filled with acidic hydrothermal fluids (Kanda et al. 2019).
The volcano displays different styles of eruptions together with persistent degassing (Ono et al. 1995).Seismic data have revealed the existence of different types of signals: from Long-Period (LP), with characteristic periods of 0.5-1 s, to Very-Long-Period events (VLP) historically termed Long-Period Tremor (LPT) (see, e.g., Kaneshima et al. 1996;Yamamoto et al. 1999;Kawakatsu et al. 2000) with characteristic periods of seconds to tens of seconds.The repetitive waveforms of VLPs suggest a non-destructive source, whose location is about 200 m southwest of the Nakadake 1st crater, about 500 m bsl (Yamamoto et al. 1999;Legrand et al. 2000).They have been interpreted as the effect of a pressurization/depressurization mechanism in a shallow crack-like conduit due to the hydrothermal activity (Niu and Song 2020;Kawakatsu et al. 2000).The analysis of the geophysical signals associated with the volcano activity has been extensively carried out on timescales of seconds/minutes or timescales of years, looking at seismic (Niu and Song 2021a), tilt (see, Niu and Song 2021b) and SAR data (e.g., Nobile et al. 2017).
In the present work, we aim to fill the gap of knowledge of Aso activity, going to investigate the ground deformation over an hourly time scale, acquired by the tiltmeter network.To date, some studies have been conducted on this timescale (Genco and Ripepe 2010;Anderson et al. 2010Anderson et al. , 2015;;Zlotnicki et al. 2018) in other volcanic areas, revealing the presence of deflation/inflation (DI) transients characterized by different waveforms whose duration varies from minutes to hours.Such signals are generally connected to processes occurring at shallow depths in the magmatic plumbing system (Anderson et al. 2015).For instance, at Mt.St. Helens thousands of cyclic, near-event ground tilt transients lasting from minutes to hours were detected during the 2004-2008 eruption and were characterized by an asymmetric pattern of rapid tilt away from the vent followed by a gradual reversal (Anderson et al. 2010).Tilt vectors pointed southward to the 1980s lava dome.The absence of detectable tilt outside the crater suggests a shallow source (< 1 km), likely related to cyclic conduit pressurization/depressurization associated with dome slip and gas loss (Anderson et al. 2010).Tilt variations corresponding to DI episodes on hourly timescales have been observed at Taal Volcano (Philippines) and associated with gas release, fluid transfers and other processes in the hydrothermal system (Zlotnicki et al. 2018).Moreover, tilt signals have been found at Kilauea volcano (Anderson et al. 2015), where DI events of different shapes (DID-type, V-type, U-type) have periods of hours to days.They have been associated with perturbations in magma flow and pressure conditions in the shallow magmatic plumbing system and with the cyclic level variations of the summit lava lake.At a smaller timescale, Genco and Ripepe (2010), at Stromboli volcano (Italy), observed saw-toothed tilt deformations that appeared as repetitive waveforms, characterized by periods of hundreds of seconds.Such events were interpreted as a cyclic mechanism of inflation/deflation associated with a charge/discharge magma pressure controlling the explosions.A peculiarity of such deformation is that it is very local, almost disappearing at distances of 1 km.At Onikobe geyser, Nishimura et al. (2006) found similar inflation/deflation events, whose repetitive occurrence allowed the prediction of the duration time of an eruption.
The previous outcomes testify that the analysis of ground deformation over timescales of hundreds of seconds to days and the extraction of peculiar transients provide useful insights for the characterization of the volcano dynamics worldwide.The present study follows this line of research at Aso Volcano and adds a piece of information on the general behaviour of the Aso volcanic system.In addition, this investigation helps characterize Aso activity in the period 2011-2016.Our findings could be useful for tracking the dynamics of the shallow hydrothermal system.Indeed, monitoring hydrothermal activity, as well as magmatic activity, is also important since it can potentially help to better evaluate if the volcanic system is primed for phreatic and/or phreatomagmatic eruptions.Finally, our analysis takes advantage of some approaches based on polarization analysis and hierarchical clustering applied to tiltmeter time series.

Outline of the 2013-2016 activity of Aso Volcano
The period between September 2013 and May 2014 was characterized by high fumarolic activity in the Nakadake 1st crater (Cigolini et al. 2018).The crater lake, since April 2013, progressively reduced: the surface ratio of the lake water to the crater bottom fell to < 10% after December 2013 (Miyabuchi and Hara 2019) and fumaroles appeared at the crater bottom (Shinohara et al. 2018a).Between mid-January and late February 2014, several ash-gas emissions and small phreatic and low-energy phreatomagmatic eruptions took place at the Nakadake 1st crater (Ichimura et al. 2018).Strong variations in the amplitude of the continuous seismic tremor corresponding to shallow volcanic activity occurred between December 2013 and January 2014.In particular, a significant increase followed by a sudden decrease was observed before the opening of a new vent inside the crater.A similar pattern was observed in the SO 2 flux, which was about 500 t/d before 2013 and increased to 500-2000 t/d in late 2013, after the drying of the crater lake water (Shinohara et al. 2018a).Ichimura et al. (2018) estimated the source locations of continuous tremor from December 2013 to January 2014.The tremor source locations were distributed in a roughly cylindrical region of 100-150 m in diameter ranging from the surface to a depth of 400 m beneath the crater.The authors interpret the tremor variations before the January 2014 events as related to an increase in gas-rich volcanic fluids causing an enlargement of the conduit zone, followed by the migration of further magmatic fluid through other pathways, which resulted in a subsequent ash-gas emission.Upward and downward migration of the estimated source depth was also identified and was associated with changes in volcanic activity.These phenomena were considered as the first precursory anomalies for the latest 2014 Strombolian eruption (Ichimura et al. 2018).
The whole period between September 2013 and May 2014 was also characterized by significant gravity decreases (up to − 100 μGal) in the active area (Sofyan et al. 2016), while the volcanic activity resurfaces.This mass loss in more than half a year was the largest one ever observed throughout the period 2010-2014 (Sofyan et al. 2016), indicating the occurrence of hydrothermal dynamic processes at the subsurface.It was interpreted as a mass movement to the deeper area or a large mass loss due to a massive evaporation because of the high temperature, while the water level in the crater lake was decreasing rapidly.After the first 2 months of 2014, the phreatic phase ceased.The crater lake completely dried up in July 2014 and remained in that state up to June 2015, leaving an incandescent vent near the crater centre.Between May and December 2014, small variations in gravity measurements were observed.
The Strombolian phase started with a minor eruption on August 30, 2014 (Miyabuchi and Hara 2019).On November 25 there was a large Strombolian eruption with continuous ash emissions accompanied by intermittent explosions.Active eruptive events continued until May 2015 (Miyabuchi and Hara 2019).Afterward, the interaction with a small amount of hot water, entering the pyroclastic cone in June 2015, was followed by small eruptions in August and September.From the beginning of the eruption to the end of January 2015, the cumulative erupted mass increased at a high discharge rate (2.2 × 10 4 tons/day), then after February 2015, it decreased to 0.6 × 10 4 tons/day (Shinohara et al. 2018a;Ohkura et al. 2019).By comparing electromagnetic data before and after the Strombolian eruption of November 2014, Minami et al. (2018) reported an increase in resistivity beneath and outside the crater, likely due to the reduction of conductive groundwater in the upper part of a shallow aquifer (< 800 m asl).Specifically, the eruption broke the clay cap below the crater, decreasing the amount of hydrothermal fluid by enhanced pressure and temperature in the vicinity of the magma conduit, and by the generation of hot gases.
After the end of May 2015, the volcanic activity consisted of minor eruptions, then moderate phreatomagmatic explosions occurred on 14 September and 23 October 2015 (Morita et al. 2019).Since March 2016, a small amount of hot water appeared in the crater and the surface area of the lake increased to 70% of the total area of the crater bottom (Morita et al. 2019).During March and September 2016, minor mud emissions from the crater bottom were observed.A large and destructive earthquake (Kumamoto earthquake, Mw = 7.3) occurred on 16 April 2016 at about 40 km southwest of Aso with many large aftershocks (Yagi et al. 2016).The co-seismic rupture and subsequent crustal deformations associated with these earthquakes concentrated on the southwest side of the volcano and interacted with the magmatic system, likely causing magma and high-pressure hydrothermal fluid migration in the pre-existing cracks (Yagi et al. 2016;Ozawa et al. 2016;Yamada et al. 2019).On 8 October 2016, a violent explosive phreatomagmatic eruption occurred: the plume height reached 13-14 km asl and the total mass of volcanic ash was estimated to be 6.0-6.5 × 10 8 kg (Ohkura et al. 2019).
The eruptive and inter-eruptive phases which alternated in 2013-2016 have been accompanied by LP and VLP events located at a depth between 400 and 1500 m below the crater.They are preceded by deformation events, typically with a characteristic time scale of 50-100 s, located at a depth of about 3 km (Niu and Song 2021a,b).This interconnected seismic and deformation activity has been interpreted as the result of magma ascent, storage, and discharge in an over-pressured magmatic-hydrothermal system.After the strombolian eruption of November 2014, these episodic deformation events appear also to involve the ascent of gas or/and volatiles.

Data and methods
We analysed tilt data from two borehole tiltmeters, N.ASIV (hereinafter ASIV) and N.ASHV (hereinafter ASHV), belonging to the volcano observation network V-net (Tanada et al. 2017); the instrument resolution is ~ 1 nanoradian (Fig. 1).The borehole tiltmeters are distant from the crater lake a few kilometres, ASIV about 4.4 km and ASHV about 3.7 km.
Data have been recorded in the period January 2011-December 2016 (Fig. 2a), during which the volcano has experienced quiescence (2011)(2012)(2013), phreatic (January 2014), Strombolian (November 2014), and two phreatomagmatic (September 2015, October 2016) eruptions.Thus, the observational time range covers different phases of volcanic activity and represents a unique opportunity to study the transition mechanisms among different styles.
We perform time-series analyses in order to calculate the polarization parameters of the background ground tilt.Then, we extract peculiar long-lasting transients by compiling a catalogue.Such transients are characterized in terms of waveform cross-correlation, particle tilt pattern, energy and time distributions.
The polarization of the ground tilt is estimated following the approach described in De Lauro et al. (2018).The azimuth of the polarization vector, defined as the angle between the projection on the horizontal plane of the polarization vector and North (measured clockwise), is calculated by means of diagonalization of the covariance matrix.An alternative way of estimating the parameters can be done following the approach of De Lauro et al. ( 2016).The azimuth values, which provide the direction of the tilt in the horizontal plane, span the interval 0-180° due to the 180° ambiguity that affects the estimate of the polarization vector.In order to investigate the spatial pattern of the ground tilt on hourly timescales, the time-series are detrended and filtered between 0.5 and 8 h by applying a third-order Butterworth zero-phase filter.Using the NS and EW components, the azimuth is evaluated over sliding time windows with a length equal to the maximum period of the filtering band.
The recognition of peculiar events in a continuous time-series exploits two sequential steps: (1) detection of the potential events and compilation of a catalogue; (2) classification of the events based on the cross-correlation.This task is carried out by applying an automatic picking algorithm on the basis of the peak amplitude: a pick occurs when the absolute amplitude of the tilt pulse exceeds a fixed threshold.The choice of the threshold value is of course a compromise between background noise and the amplitude of the signal; larger thresholds will enable the picking only of the most energetic signals, while smaller values will also pick noisy low-amplitude events.On the basis of the pickings, signal waveforms are extracted using a time window whose length depends on the event duration.
The individuation of families is achieved by using a hierarchical clustering procedure based on cross-correlation to highlight waveform similarities (Petrosino and Cusano 2020).The correlation matrix calculated on the extracted waveform samples is rearranged according to the cluster distance estimated by the hierarchical clustering algorithm.The measure of the distance is 1 − C ij (i, j = 1,…, number of waveforms), where C ij is the maximum correlation coefficient between the ith and jth waveform samples so that the event pairs with high cross-correlation are grouped together (Greenfield et al. 2019).
The inter-event interval (IE) distribution will be obtained by evaluating the variability index I v , defined as , where σ �t is the standard deviation of all IEs and t is their mean value.We note that the value of I v can be a marker of the dynamic process.Indeed, I v = 1 indicates a Poissonian process, I v > 1 denotes a clustered process, and I v = 0 underlines a periodic one (Massey 1951).In addition, we applied the Kolmogorov-Smirnov (KS) test (Massey 1951) to check the hypothesis of Poissonianity in the IE distribution function.It is possible to compare the empirical (Φ E ) and the theoretical exponential cumulative distribution function � T (�t) = 1 − e − �t (Massey 1951;Rice 1995), where is a known parameter of the Poisson process.The null hypothesis is that the functions are equal; a measure of the discrepancy can be given by a coefficient I v2 = 1 − I v , complementary to I v , whose value is as close to zero as possible for a Poisson process.

Results
In each subsection, we briefly describe the primary outcomes derived by the polarization, hierarchical clustering, and statistical analyses of the ground tilt.Subsequently, we provide our interpretations of the deformation pattern in the framework of the geological setting and volcano dynamics.

Polarization analysis of the ground tilt
The results of the polarization analysis of the ground tilt time-series recorded in 2011-2016 show that in the time interval which roughly goes from September 2013 to June 2015, the particle tilt orientation at ASIV has a clear variation (Fig. 2b, c): the azimuth values concentrate in the 30°-100° range.
In this time interval, variations of the tilt amplitude are particularly evident on ASIV.Such amplitude variations correspond to a sequence of tilt pulses detected only at ASIV (see green ellipses in Fig. 2c).These signals have repetitive waveforms with a large negative pulse of about 2-3 h of duration on both components.We refer to these Very-Long-period Tilt Pulses as VLTPs.The most energetic VLTPs show some oscillations lasting up to 6-8 h; sometimes they are very close in time, merging in a nearly continuous "puffing".The comparison of the filtered and unfiltered waveforms makes us confident that the filtering operation preserves the overall shape of the original pulses (Additional file 1: Figure S1).
The azimuth of the VLTPs is about 45° (see the green ellipses in the bottom frame of Fig. 2c) which corresponds to ground tilt along the NE-SW direction.

Picking of long-period tilt events and catalogue compilation
We applied the automatic picking algorithm to the filtered tiltmeter time-series, fixing an amplitude threshold equal to two times the minimum value (among the four traces) of the Root-Mean-Square (RMS) calculated over the entire length.In addition, to avoid the overlap between the tilt pulses, we enabled the picking only if the time separation between two adjacent events was at least 8 h, considering their long duration.We used the automatic picking separately on the NS and EW components of the two tiltmeters, because there was no a priori information about the pulse energy partitioning between the two components.
The final catalogues for both the tiltmeters were obtained by merging the different picked times (not necessarily coincident) on both components.At ASIV the catalogue consists of 161 events, mostly concentrated in September 2013-June 2015.In this time interval, few signals are picked on ASHV tiltmeter, for which the final catalogue contains 79 events.We verified that almost all the pickings for ASHV are not coincident in time with those of ASIV, in particular in the time range September 2013-June 2015.This suggests that ASHV is detecting different kinds of signals with respect to ASIV.On the basis of the catalogues, we extracted signal waveforms from ASIV and ASHV filtered time-series by using an extraction time window with a length fixed to 480 min, starting 180 min before the timing of the maximum peak amplitude.

Hierarchical cluster analysis
The recovered waveforms were further analysed to identify events with VLTP features and eventually to discriminate them from other (spurious) signals.In order to highlight waveform similarities, we performed a hierarchical cluster analysis based on cross-correlation (Petrosino and Cusano 2020).The results are shown in Fig. 3.
Considering a correlation threshold C t = 0.9, on the ASIV NS component the algorithm separates one large cluster (orange branches of the dendrogram in Fig. 3, left panel) of 82 events (almost 50% of the entire dataset) from a few minor clusters, typically formed by less than 20 events (see, e.g., red dendrogram branches in Fig. 3, left panel) and isolated events.Similarly, on the ASIV EW component, the major cluster contains 71 signals (red dendrogram branches in Fig. 3, right panel), whereas the others (see, e.g., green dendrogram branches in Fig. 3, right panel) are composed of a few events (< 20).Looking at the correlation matrix, the major cluster could actually be merged, at the higher node (C t = 0.88) of the hierarchical cluster tree with the adjacent minor one composed of 19 events, thus forming an ensemble of 90 signals.In fact, using a lower correlation threshold C t = 0.88, the two clusters are no longer separated; this depends on the fact that the signal-to-noise ratio on the EW component is lower than the NS.From the comparison of the arrival times, we verified that all the 82 clustered events on the NS component are included in the ensemble of 90 signals on the EW component.
On the other hand, at ASHV tiltmeter, the cross-correlation analysis on both the components highlights mostly isolated events and a few small clusters (Additional file 1: Figure S2); the most numerous one, formed by 22 events, appears on the EW component (cyan branches of the dendrogram in Additional file 1: Figure S2 right panel).
In order to ensure statistical robustness, we carried out careful waveform and particle motion analyses on the cluster with the largest population at ASIV, thus on the 82 events present in both components at the same time.Such clustered signals have the polarization characteristics of the VLTPs previously defined.The timecoincident waveforms on the NS and EW components as well as their normalized stacking (with the 95% confidence intervals estimated by a bootstrapping procedure) are shown in Fig. 4a, b.The waveforms appear like deflation-inflation (DI) events (see, e.g., Anderson et al. 2015;Genco and Ripepe 2010): the signal associated with a tilt decrease indicates an inward motion of the tiltmeter towards the crater area (deflation), whereas the tilt increase indicates an outward motion away from the crater area (inflation).The amplitudes of the negative pulse range from 0.01 to 0.2 microrad on both components.A small positive pulse seems to precede it but looking at both raw and filtered waveforms, it is really difficult to distinguish if it is due to the low signal-to-noise ratio or a negligible artefact introduced by the acausal filter (see Additional file 1: Figure S3).The stacked signals from ASHV in the same time windows as those of ASIV (Fig. 4e) appear as fluctuations and do not reveal any characteristic waveform corresponding to the VLTP detected at ASIV.
Looking at the waveform, the shape of VLTPs that we observe at Aso Volcano is similar to that of the Kilauea V-type events identified by Anderson et al. (2015), although they have a different deflation rate (0.4-1 day).An interesting feature of Aso VLTPs is that, despite the different peak amplitudes, the deformation rate is similar: the steep tilt decrease (deflation phase) occurs in about 30 min, then the original deformation is recovered in about 120 min (inflation phase).In other words, the duration of such signals is independent of their amplitude, which implies that, whatever the energy involved in the process, the dynamical model underlying the phenomenon is the same.Recently, Ripepe et al. (2021) at Stromboli volcano (Italy) observed that both the inflation amplitude and the duration, in the DI tilt signals, scale with the magnitude of the explosions, whereas DI tilt events with variable duration were detected at Mt.St. Helens (see, Anderson et al. 2010).
By composing the NS and EW-filtered components of ASIV, we obtained the particle tilt in the horizontal plane for the clustered 82 VLTPs (Fig. 4c).The distribution of the 82 azimuths calculated by the polarization analysis provides a mean value of 52° ± 18°.The large tilt decrease corresponds to a deflation towards SW (i.e. the Nakadake 1st crater area), common to most of the events.Such deflation could be preceded by a small positive pulse (Fig. 4c and Additional file 1: Figure S3), as in the case of the V-type events (Anderson et al. 2015), which are sometimes preceded by small inflationary bumps just prior to the onset of deflation.
The same analysis performed on the 22 clustered events at ASHV provides waveforms with a similar rate of deformation (Fig. 5a, b).They also look like DI events, starting with a deflation towards NE (i.e. the Nakadake 1st crater area) and followed by an inflation towards SW (Fig. 5c, d).The stacked signals from ASIV in the same time windows as those of ASHV (Fig. 5e) do not show any characteristic waveform corresponding to the DI events detected at ASHV.

Statistical analysis of the clustered VLTPs
We calculated the time and frequency distributions of the amplitudes of the clustered 82 VLTPs recorded at ASIV.Indeed, both the inter-event interval distribution function (Cox and Lewis 1966) and the energy distribution can give information on the dynamic process of generating VLTPs.Such analyses were already successfully applied in the seismological domain (see, Bottiglieri et al. 2005;De Lauro et al. 2008, 2009;Sandanbata et al. 2015).The results show that most of the VLTPs occur between September 2013 and February 2014, then they become less frequent (up to June 2015).Their amplitude increased throughout June 2013-June 2015, with a steeper increase after November 2014 (Fig. 6a).Over this time period, the amplitude distribution indicates that the generation process involves tilt variations in the range of 0.04-0.18microrad; during September 2013-February 2014 tilt amplitudes are almost steady (and mostly concentrated) between/around 0.04-0.05microrad (Fig. 6a), suggesting the existence of a preferred energy scale (Matoza and Chouet 2010;De Lauro et al. 2012).
The IE analysis shows that the IEs are typically less than 4 days, especially in September 2013-February 2014, after which the IEs increase, indicating less frequent tilt episodes.Taking into account all the VLTPs, the variability coefficient I v is equal to 1.6, revealing an overall process that appears to be clusterized on a long time scale (see, Fig. 6b).We repeated the same analysis focusing on the cluster, inside of it IEs do not exceed 20 days.In this case, I v lowers attaining the value of 0.98, K-test = 0.097 passed with a significance level of 0.1, and I v2 = 0.14.
Figure 7 shows the empirical cumulative distribution function (CDF) compared with the theoretical one for an exponential distribution, together with the residuals.Combining these indicators together, they reveal an almost Poissonian nature of the cluster-generating process as discovered for other volcanoes worldwide, such as Stromboli (Bottiglieri et al. 2005;De Lauro et al. 2008), Erebus (De Lauro et al. 2009), and Campi Flegrei (De Lauro et al. 2012).The mean value of the Poissonianity, i.e. the mean value of the IEs (�t) , is 4.6 d, and the rate is such that 1/ = 4.93 d .In summary, the VLTPs belong- ing to the cluster occur randomly in time, but on average each VLTP happens approximately every 4.6 d.
On the basis of the statistical properties, we could identify three periods (see light coloured boxes in Fig. 6 Considering the evolution of Aso eruptive dynamics (described in "Outline of the 2013-2016 activity of Aso Volcano" section), P1, P2, and P3 roughly correspond to the phreatic phase, the inter-eruptive phase, and the Strombolian eruption phase, respectively.
The time and frequency distributions of the amplitude of the 22 events belonging to the small cluster of ASHV show that most part of them occurred between April and July 2016 (Additional file 1: Figure S4), suggesting that they could be related to the strong phreatomagmatic activity which characterized that period (Ohkura et al. 2019), or to crustal deformation and rupture processes induced by the Kumamoto earthquake sequence (Yagi et al. 2016;Ozawa et al. 2016;Yamada et al. 2019).Unfortunately, the small size of the sample does not allow us

Discussion
In the present work, we have tackled the problem of characterizing the activity of Aso Volcano by looking at the scale of hours.In particular, we have studied the ground deformation acquired by two stations of the V-net tiltmeter network, analysing the continuous recordings from January 2011 to December 2016.The period chosen is characterized by a variety of eruptive styles.The characterization of ground deformation in terms of polarization analysis together with a clustering methodology allowed the identification and extraction of VLTPs.VLTPs mainly belong to the same cluster observed from September 2013 to June 2015 and well recorded at ASIV.The stacked repetitive waveforms show a large negative pulse of about 2-3 h, indicating a clear deflation/inflation behaviour.The inter-event statistics of such a cluster provides evidence that the generation process of VLTPs is basically a Poisson process with a rate of about every 5 days, and a preferred energy scale in the period September 2013-February 2014 where the tilt amplitudes are almost steady coming from the estimated log-normality.
The detection of VLTPs at the sole ASIV site indicates that we are tracking a local phenomenon related to shallow dynamics.In fact, the size of the deforming region generally depends on the size and/or depth of the source: a deeper source causes deformation over a larger area of the surface (Pritchard et al. 2019), whereas surface mechanisms can generate DI events locally circumscribed, as also observed in other volcanoes (Anderson et al. 2010;Genco and Ripepe 2010).This also seems to rule out the hypothesis of a deep magmatic process as the case of the episodic deformations detected by all the tiltmeter network, which have been related to the transport of discrete magma batches from a deep overpressured chamber to a pre-eruptive storage zone at 3 km of depth (Niu and Song 2021b).
In the case of our VLPTs, the most probable origin appears to be in the hydrothermal system of the volcano and in the shallow fluid circulation; the structural features also account for the local nature of the phenomenon.In fact, the central area of the Aso caldera is characterized by high velocities (Vs > 3 km/s, H1 anomaly in Huang et al. 2018) at surface depths of 1 km surrounded by low velocities (Vs < 2 km/s, L1 anomaly in Huang et al. 2018) from the caldera surfaces to a depth of approximately 2.5 km (Huang et al. 2018).Velocity anomalies match the resistivity anomalies retrieved from magnetotelluric data (Hata et al. 2016;Matsushima et al. 2020).The high-velocity body H1 has been ascribed to the presence of a consolidated igneous rock with high density and resistivity.Conversely, the L1 low-velocity anomaly, which corresponds to an area of low density and resistivity, has been interpreted as a shallow hydrothermal reservoir (or pathways containing volcanic fluids) correlated to the surface geothermal activity of the Nakadake 1st crater.The L1 anomaly extends north-eastward of Nakadake 1st crater, as far as ASIV site, becoming progressively shallower (Huang et al. 2018).ASIV tiltmeter is located between the shallow low-velocity and high-velocity zones L1 and H1, in proximity of the lateral velocity contrast between the reservoir of geothermal fluids and the igneous rock, thus a favourable position (rather than Fig. 7 On the left, the comparison between the experimental (blue line) and theoretical (red line) CDFs according to the KS test; on the right, the residuals between the curves are within 10% in agreement with the hypothesis of Poissonian process ASHV site) for detecting the VLTPs induced by the activity in the hydrothermal system.All the primary structural information related to Aso volcano is sketched into a map shown in Fig. 8, where the location of the hydrothermal system, VLP and deep inflation sources, velocity and conductivity structures are reported.
We interpret the observed VLTPs occurred between September 2013 and June 2015 in light of the activity of Aso Volcano.Indeed, previous geophysical observations (described in Section Outline of the 2013-2016 activity of Aso Volcano) testify intense hydrothermal activity and shallow fluid migration before and during the minor phreatic phase in September 2013-February 2014 and a progressive decrease after the Strombolian eruption of November 2014.The VLTPs throughout September 2013-June 2015 accompany such a variety of states of volcanic activity with different rates.In particular, frequent and moderate rates are observed in P1 and P2, then the VLTP activity starts to decrease until some months after the Strombolian eruption (P3).In this context, we contemplate that the DI tilt episodes are possibly associated with the depressurization/pressurization of the shallow hydrothermal system occurring throughout the observational period.The VLTPs can be interpreted as the result of the stress release related to a decreasing fluid supply and/or lowering pressure.The independence between the amplitude and duration of VLTPs suggests that the time of discharge and recharge of a reservoir is always the same, either the fluid supply or the pressure involved.A pressure variation can be considered as initial condition in the volcano dynamics, not affecting the time evolution, thus suggesting a quite linear behaviour.
Taking into account these considerations, the VLTPs could be generated by a simple fault-valve mechanism related to the intermittent discharge of fluid from the hydrothermal system that is in a nearly critical state.The link between VLTPs and fluids accumulating in faults and fractures is not necessarily a peculiarity of Aso Volcano.Indeed, a similar phenomenon has been observed at Campi Flegrei (Italy), where episodic ground tilts were Fig. 8 SSW-NNE cross-section displaying the main features of Aso including the crater lake, the source location of VLPs and deep inflation/ deflation events, the hydrothermal system, the aquifer, some relevant velocity and conductivity anomalies, as well as the location of the tiltmeters on the horizontal plane ascribed to geological structures which become hydrologically active as a result of fluid increase (Petrosino et al. 2020).On the other hand, fault-valve mechanisms related to the periodic pressurization of the hydrothermal reservoir have been suggested in other volcanoes (Cucci et al. 2017;Calderoni et al. 2019;Cusano et al. 2020;Falanga et al. 2021;Niu and Song 2021b).In such conditions, a slight pressure fluctuation can trigger fluid outflow with a consequent net volume decrease and depressurization of the reservoir.This corresponds to the deflation part of the VLTPs.Subsequently, the system is gradually refilled by the circulating fluids (inflation phase), recovering its starting condition/initial state.The hydrothermal outflow can occur through the shallow fractured rocks, whereas the fluids can accumulate in discontinuities such as faults and/or fractures, with the effect of changing the local stress and producing fluid-induced tilt.
As evidenced by Shinohara et al. (2018b), hydrothermal fluids can flow from the hydrothermal system towards the volcanic conduit.Such mechanism also persists when the crater lake is dried (as well as during the Strombolian activity), as remarked by the study of salt shells and lumps that were observed also during that period.This suggests the persistence of the hydrothermal fluid circulation beneath the crater floor, even during the course of magmatic eruptions (Shinohara et al. 2018b;Kanda et al. 2019).
The repetitive pattern of the VLTPs indicates that their generation mechanism has remained almost stationary throughout September 2013-June 2015, always driven by a Poisson process.The rate and energy variations observed during the three eruptive and intereruptive phases are likely related to the efficiency of the source process.After the climax of the phreatic phase in January 2014, the progressive drainage of the crater lake indicates a lesser supply of hydrothermal fluids, which could explain a lower VLTP rate.The decrease of fluids in the hydrothermal system was probably even more conspicuous during the Strombolian phase, when the fault-valve mechanism was still active but with an even lower occurrence.This is supported by the energy distribution and IE analyses which provide evidence of a dynamical process (i.e. the fault-valve mechanism), which remains the same over different time scales.In other words, the valve opening/closure occurs more frequently involving lower energy, or less frequently but with higher energy.Overall, these observations testify a progressive depressurization of the hydrothermal system throughout September 2013-June 2015, first more continuously and involving constant energy release, then in a more discontinuous and energetic style.At the end of the Strombolian phase, the fault-valve mechanism stops, likely due to the drastic decrease in the amount of conductive groundwater in the upper part of the reservoir (Minami et al. 2018), which probably was not able to reach an overpressured condition for generating VLTPs.

Conclusions
In the present work, we have studied the ground deformation at Aso Volcano using a 6-year long (January 2011-December 2016) tiltmeter time series.We identified typical events, VLTPs which are radially polarized and mainly clustered in the period September 2013 to June 2015.They appear like deflation-inflation events, with a large negative pulse of about every 2-3 h.The inter-event statistics indicate a Poissonian-like generation process with a rate of about 5 days and a preferred energy scale in the period September 2013-February 2014.Such events are local as they have been recorded only at ASIV site.
The VLTPs observed in September 2013-June 2015 at Aso volcano are mostly related to the hydrothermal activity, as inferred from the available geophysical, geochemical and volcanological information.These signals can be considered as a signature of the depressurization/pressurization of a shallow source, unraveling surface processes involving geothermal fluids.The outcomes of the VLTP analyses shown in the present paper can be useful for tracking the dynamics and the evolution of the shallow hydrothermal system of Aso volcano.
distribution of the IEs; on the right IE temporal pattern and event cumulative number (red line) as a function of time.

Fig. 1
Fig. 1 Location of Aso Volcano and of the two tiltmeters belonging to the volcano observation network V-net (map data ©2022 Google)

Fig. 2 a
Fig. 2 a Ground tilt recorded along 6 years (2011-2016) at ASHV and ASIV instruments.The sign convention is positive for tilt towards the north and the east.b The first two panels show the filtered tiltmeter time series at ASHV (red = EW, black = NS) and ASIV (blue = EW, black = NS); the azimuth values obtained from polarization analysis at ASHV (red) and ASIV (blue) are plotted versus time in the third panel.c A zoom on 10 days (in September-October 2013) is shown

Fig. 3
Fig. 3 Correlation matrix ordered according to the hierarchical clustering for the NS (left) and EW (right) components of ASIV tiltmeter.The red colour corresponds to the highest values of correlation.The coloured branches of the hierarchical cluster tree (dendrogram) indicate correlation values greater than 0.9; their height is a measure of the correlation between two events or clusters.The orange dendrogram branches in the left panel and the red ones in the right panel correspond to the major clusters of 82 and 71 events, respectively

Fig. 4 a
Fig. 4 a Waveforms of the 82 clustered VLTPs (ASIV NS filtered component) and stacking of the signals (normalized amplitude) with bootstrap uncertainties at the 95% confidence level.b The same for ASIV EW-filtered component.c Particle tilt on the horizontal plane of the 82 clustered VLTPs.d Azimuth distribution obtained from the polarization analysis of the 82 clustered VLTPs.e Stacked waveforms (normalized amplitude) from ASHV (NS and EW components) in the same time windows as those of ASIV

Fig. 5 a
Fig. 5 a Waveforms of the 22 clustered events (ASHV NS filtered component) and stacking of the signals (normalized amplitude) with bootstrap uncertainties at the 95% confidence level.b The same for ASHV EW filtered component.c Particle tilt on the horizontal plane of the 22 clustered events.d Azimuth distribution obtained from the polarization analysis of the 22 clustered events.e Stacked waveforms (normalized amplitude) from ASIV (NS and EW components) in the same time windows as those of ASHV ): P1.September 2013-February 2014 characterized by VLTPs with high-rate and almost constant lowamplitude; P2.March 2014-October 2014 lower rate and smooth amplitude increase; P3.November 2014-June 2015 low rate and steep amplitude increase.

Fig. 6 a
Fig. 6 a On the left, frequency distribution of the amplitudes of 82 clustered VLTPs (ASIV NS component); on the right, VLTP amplitude and cumulative number (red line) as a function of time.The yellow line marks the beginning of the phreatic phase in January 2014, the purple line corresponds to the Strombolian eruption of 25 November 2014.The three periods P1 (September 2013-February 2014), P2 (March 2014-October 2014) and P3 (November 2014-June 2015) are highlighted by light blue, green and yellow boxes, respectively.b On the left, frequency distribution of the IEs; on the right, IE temporal pattern and VLTP cumulative number (red line) as a function of time