Temporal variation in the depth of the magma surface at Aso volcano in 2014–2015

Monitoring the depth of the magma surface at open-vent volcanoes can be a practical tool to infer temporal variations in the magma supply during an eruption. We focus on the magmatic eruption of Aso volcano in 2014–2015 to estimate the temporal change in the depth of the magma surface, and show that this needs to be coupled with an understanding of the shallow conduit geometry if it is to be done in a representative manner. The eruption lasted 5 months from November 2014 and ending with a crater floor collapse in May 2015. During the eruption, we recorded seismo-acoustic waveforms related to frequent Strombolian explosions. The infrasound signals show several distinct peak frequencies derived from acoustic resonance inside the vent. We estimate the depth of the magma surface using the time delay of seismo-acoustic signals and the peak frequency of infrasound signals. In addition, the temporal variation in the shape of the conduit is constrained by the overtone frequency of the acoustic resonance. From the beginning of the eruption to early-January 2015, the magma surface was located at a depth of ∼\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\sim $\end{document}200 m, and the conduit was a cylindrical pipe. Later, between January and February, the magma surface rose to ∼\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\sim $\end{document}120 m, and the shape of the conduit changed to a conical frustum flaring inside. This finding indicates that the magma was injected into the shallow conduit and that it heated and weakened the conduit wall near the magma surface. Before the cessation of the magmatic eruption, the magma surface dropped by approximately 70 m. This magma drainage and, primarily, the instability of the conduit shape caused the crater floor to collapse. We show the possibility of tracking and assessing the depth of the magma surface and the shallow conduit geometry even with limited seismo-acoustic observations.


Introduction
The depth of the magma surface in the conduit fluctuates to reflect the behavior of plumbing systems, resulting in a change in the eruptive style (cf. Ripepe et al. 2002;Harris and Ripepe 2007;Patrick et al. 2019). Temporal changes in depth have been observed at several open-vent volcanoes such as Kilauea (e.g., Tilling 1987;Patrick et al. 2015), Villarrica (e.g., Palma et al. 2008;Richardson et al. 2014;Johnson et al. 2018a), Stromboli (e.g., Ripepe et al. 2005Ripepe et al. , 2015Marchetti et al. 2008), and Etna (e.g., Sciotto et al. 2013;Cannavo' et al. 2019;Watson et al. 2020). Such changes in depth can be induced by a change in the magma supply rate, a change in the density of the ascending mixture, the formation of a new vent at the flank, and/or conduit constriction (e.g., Patrick et al. 2016;Hamling et al. 2019). Periodic depth change can also suggest cycles of fresh magma intrusion (Richardson et al. 2014), and sudden drops in the magma surface can trigger crater collapse and phreatic eruptions (Dvorak 1992;Anderson et al. 2019). In addition, changes in shape of the upper-most portion of the conduit can influence the style of explosive activity and/or the dynamics of ejecta (e.g., James et al. 2004;Capponi et al. 2016;Spina et al. 2019;Salvatore et al. 2020;Schmid et al. 2021). Therefore, monitoring the depth of the magma surface and the shape of the shallow conduit can elucidate the mechanism of the transition in activity style.
Aso volcano began its first magmatic eruption in approximately 20 years in November 2014 (Yokoo and Miyabuchi 2015). This activity lasted for 5 months, and was characterized by Strombolian explosions and ash emissions, which are the ordinary eruptive style of Aso volcano (Ono et al. 1995). A collapse of the crater floor occurred at 22:04 JST on May 3, 2015 (all dates and times in this paper are in Japan Standard Time, UT+9 h), which consequently brought an end to the magmatic eruption. Although the magma surface inside the upper conduit could not be seen from the crater rim during the eruption, the depth of the magma surface was estimated using several methods. Yokoo et al. (2019) estimated the depth as 118 m using the peak frequency of infrasound signals in late-March 2015, assuming that the peak frequency was derived from the acoustic resonance of the space above the magma surface in the conduit. Ishii et al. (2019) estimated the depth as <400 m using the time delay between seismic and infrasound signals observed for Strombolian explosions during late-April 2015. Ishii and Yokoo (2021) estimated the depth using both the peak frequency of infrasound signals and the time delay between seismic and infrasound signals. In addition, Ishii and Yokoo (2021) constrained the shape of the shallow conduit from the frequency ratio of the fundamental mode and a higher mode of acoustic resonance. As a result, it was found that the near-surface conduit pipe was a conical frustum flaring inside and that the magma surface was located at a depth of 40-200 m in late-April 2015. However, the change in the depth during the 5 months of activity has not yet been considered. Estimating the temporal variation in depth should provide answers to the following questions: (1) How variable was the rate of magma supply during the 5 months of eruptive activity?
(2) Why did the eruption end with collapse of the crater floor?
This study thus aims to estimate the depth change in the magma surface during the 2014-2015 eruption of Aso volcano. We apply the method of Ishii and Yokoo (2021) to obtain the depth from the seismo-acoustic data recorded during the 5-month activity. In applying this method, we consider the temporal change in the conduit shape for the depth estimation. Temporal variation in conduit shape has rarely been considered in depth estimations (e.g., Richardson et al. 2014;Spina et al. 2015;Johnson et al. 2018a), even though it significantly affects the peak frequency of the infrasound signals (Watson et al. 2019).
Based on the estimated temporal change in the depth of the magma surface and the shape of the conduit, we elucidate the relationship between magma supply and volcanic activity during Aso's 2014-2015 eruption.

Activity at Aso volcano and the 2014-2015 eruption
Aso volcano is located in southwest Japan. In the last 80 years, its eruptive activity has occurred only at Nakadake 1st crater, which is located at the most northern part of seven craters of Aso volcano (Fig. 1a). During the recent non-eruptive periods, the Nakadake 1st crater contains a hot water lake ( Fig. 1e; Terada et al. 2012). Although a crater lake had existed in the Nakadake 1st crater since February 1993, lake water had evaporated since 2013, and it was confirmed that the crater floor was entirely dry by the survey on July 17, 2014 (Japan Meteorological Agency 2014). The magmatic eruption started on November 25, 2014, from a vent opening in the crater (cross in Fig. 1a). The principal activity during this period was Strombolian explosions and continuous ash emissions (Fig. 1b). Ash emissions continued for approximately 2 months, accompanied by the frequent ejection of bombs. The tephra discharge rate gradually decreased starting in February (Miyabuchi and Hara 2019). Although gas emissions were dominant from around March 21, ash emissions and ballistics ejections temporarily resumed from around April 16 (Fig. 1c). At this time, the vent had expanded to a radius of ∼25 m, and a thick layer of ash was deposited on the crater floor (Fig. 1f). This magmatic activity ended with a collapse of the crater floor on May 3, 2015 (Fig. 1d). The southern part of the crater floor subsided by ∼20 m (Fig. 1g).

Data acquisition and methods
We used seismo-acoustic data recorded in the vicinity of Nakadake 1st crater of Aso volcano from November 28, 2014, to May 31, 2015. A short-period seismometer (>1.5 Hz) and a low-frequency microphone (>0.1 Hz) were equipped in KAF and ACM stations, respectively (Fig. 1a). Both stations were located on the western rim of the active crater. The data were sampled at 100 Hz. A surveillance camera (a frame every 20 s) was also set at the CAM station during the eruption (Fig. 1b, c). Figure 2 shows the overview of the seismo-acoustic data. The RMS amplitude of the seismic signal increased from January to February (Fig. 2a) (Kazama et al. 2022) contained between 2 and 15 Hz (Fig. 2b). The infrasound signals were mainly composed of low frequencies (<1 Hz; Fig. 2d). Large amplitudes were observed after mid-March accompanied by high SO 2 fluxes (Fig. 2c). When a Strombolian explosion occurred, characteristic seismo-acoustic signals were observed (Fig. 2e-h;Ishii et al. 2019). In the most frequent cases, the seismic signal arrived at KAF station first, and then the infrasound signal arrived at ACM station ∼1 s later. The infrasound signal constantly oscillated at low frequencies (0.4-0.8 Hz), which was amplified during an explosion. This oscillation was observed even when no explosion occurred. Since similar characteristics were observed at multiple stations around the active vent, this frequency is not likely to be due to local propagation effects but is likely a result of the fundamental resonance mode in the conduit ). In addition, the frequency spectrum of the infrasound signal showed another peak around 2 Hz, regarded as an overtone of the acoustic resonance (Ishii and Yokoo 2021). Then, ∼0.3 s after the arrival of the leading compression phase, the high-frequency component of the infrasound signal (∼15 Hz) was superimposed on the low-frequency component . Waveforms of the infrasound signal changed during the 5 months of activity (described in the result section; Fig. 2e-h).
The explosion events were identified from the seismoacoustic signals using STA (short time average)/LTA (long time average) method (Allen 1978). The short and long time Fig. 2 Overview of seismo-acoustic signals. a RMS amplitude of seismic signals (>1.5 Hz). It is calculated using a 10-min window. The dashed line shows the date of crater collapse (May 3, 2015). b Amplitude spectra of seismic signals. It is also calculated using a 10min window. c RMS amplitude of infrasound signals (>0.1 Hz) and SO 2 discharge rates (provided by the Japan Meteorological Agency).
d Amplitude spectra of infrasound signals. e-h Highpass-filtered waveforms of seismo-acoustic signals (upper: >1.5-Hz seismic signals, lower: >0.1-Hz infrasound signals). Time 0 is the arrival time of the infrasound signal. These are observed at the time indicated by the white triangles in Figs. a and b (January 21, February 13, March 28, and April 16) windows of 1 s and 40 s were applied to the seismic data of KAF station and the infrasound data of ACM station. Seismic and infrasound events were detected separately. Seismic events were detected by the STA/LTA method for high-pass filtered waves (>1.5 Hz), where the threshold was 2. Events whose maximum positive amplitudes were <10 μm/s were excluded. Infrasound events were detected by the STA/LTA method using four bandpass filtered waves: (1) the 0.5-Hz component (0.4-0.8-Hz bandpass filtered wave), (2) the 2-Hz component (1.5-3-Hz bandpass filtered wave), (3) the 4-Hz component (3-6-Hz bandpass filtered wave), and (4) the 15-Hz component (10-20-Hz bandpass filtered wave). Unlike the detection of seismic events, these four filtered waveforms were used to make the detection more rigorous because the amplitude of the infrasound signal can oscillate even without explosions due to continuous gas and ash emissions. The thresholds of STA/LTA for these four infrasound components (0.5, 2, 4, and 15 Hz) were 1.5, 2, 2, and 2.5, respectively. When the 15-Hz component exceeds its threshold and any of the 0.5-, 2-, or 4-Hz components also exceed their threshold, then we detected the signal as an infrasound event. To reduce the effect of wind noise, the periods when the wind speed (10-min value) at the JMA weather station (Fig. 1a) exceeded 7.74 m/s were excluded from the detection . Events with a maximum positive pressure amplitude of <3 Pa in the raw waveform were excluded to prevent wrong detections. Most seismic events were detected from 4 s before to 1 s after the infrasound detection time. Therefore, we considered that when a seismic event is detected in this time range, the seismic and infrasound signals were related to the same explosion and counted these waveforms as applying to an explosion event.
We used the model of Ishii and Yokoo (2021) to estimate the depth of the magma surface from the crater floor. During the 2014-2015 eruption of Aso volcano, a distinct peak frequency of the infrasound signal f 0 was observed whenever any explosion occurred ( Fig. 2d; Yokoo et al. 2019). As many papers reported (e.g., Ripepe and Gordeev 1999;Hagerty et al. 2000;Ripepe et al. 2001), the seismic and infrasound signals accompanying explosions were observed with a time delay t (e.g., Fig. 2e-h). Now, following Ishii and Yokoo (2021), the depth of the magma surface can be estimated from these two measured parameters, f 0 and t. The two parameters can be expressed as functions of d, the depth of the magma surface from the crater floor (1140 m a.s.l.), and c, the sound velocity in the conduit (i.e., f 0 = F(d, c) and t = G(d, c)). These two functions can be combined to estimate the depth (d = H(f 0 , t)).
The peak frequency of the infrasound signal f 0 is thought to be a resonant frequency excited by explosions and/or steady gas emissions in the space above the magma surface ( Fig. 3a; Yokoo et al. 2019). The space acts as a sound resonator (Johnson et al. 2018a;Richardson et al. 2014). The resonant frequency is generally defined by the size of the resonator and the sound velocity in the conduit (Johnson et al. 2018b). Therefore, f 0 can be expressed by f 0 = F(d, c). Function F also depends on the shape of the resonator, such as a cylindrical pipe or horn (Watson et al. 2019). If the conduit pipe is a conical frustum, the resonant frequencies of the pipe, the fundamental frequency (f 0 ) and the overtone component (f 1 ), depend on the radius ratio between the open and closed ends of the frustum (a o /a c ; Fig. 3c). When a o /a c < 1, the resonant frequencies are lower than those of the cylinder. On the other hand, frequencies are higher than those of the cylinder when a o /a c > 1. The rate of change in the overtone component with a o /a c is smaller than the change in the fundamental frequency, especially in the range of a o /a c > 0.3. Therefore, the infrasound frequency ratio f 1 /f 0 can be used to determine a o /a c ( Fig. 3b; Ishii and Yokoo 2021). We constrained the pipe radius ratio using spectral peak frequencies of observed infrasound signals. Moreover, the relationship among f 0 , d, and c (f 0 = F(d, c)) was obtained by conducting a 3-D numerical simulation of infrasound propagation (Ishii and Yokoo 2021). In computing this simulation, we systematically changed the depth of the magma surface and the sound velocity in the conduit to obtain the peak frequency at the station. The temporal developments of the vent size and the conduit shape were reflected in the elevation data used for the numerical simulation. All other simulation settings were the same as those in Ishii and Yokoo (2021).
The time delay between the seismic and infrasound signals was interpreted as a propagation time difference assuming that a Strombolian explosion simultaneously excites these signals on the magma surface (e.g., Hagerty et al. 2000;Petersen and McNutt 2007;Sciotto et al. 2013). They propagate through the ground and air, respectively (Fig. 3a). Strictly speaking, the seismic and infrasound sources may not be the same. For example, the seismic source may be deeper than the magma surface (e.g., Ripepe et al. 2001;Harris and Ripepe 2007;Gurioli et al. 2014). In In the space above the magma surface, acoustic resonance is continually excited by explosions or gas emissions. The resonant frequency is observed as the peak frequency of the infrasound signal. When an explosion occurs at the magma surface at depth d from the crater floor (1140 m a.s.l.), the infrasound signal propagates in the conduit at sound velocity c and then in the air at c air . The seismic signal first propagates through the lava pond at velocity v p0 and in the ground at P-wave velocity v p . L i is the path length of the infrasound signal from the vent to the station (ACM). L s is the path length of the seismic signal from the margin of the lava pond to the station (KAF). b Relationships between the pipe radii of a resonator and the resonant frequencies.
The horizontal and vertical axes represent the ratio of the resonant frequencies (the fundamental mode f 0 and the first overtone f 1 ) and the pipe radii of the open and closed ends (a o /a c ), respectively. A solid line represents the relationship calculated after Ayers et al. (1985) assuming a o is fixed. When the resonant pipe is a cylinder, f 1 /f 0 is three (black circle). If the radius of the open end is smaller than that of the closed end, f 1 /f 0 becomes larger than three (red line). On the other hand, the radius of the open end is larger than that of the closed end, and f 1 /f 0 decreases to less than three (blue line). Dotted lines are the relationships obtained from the numerical simulation assuming the sound velocity inside the vent is 300 and 600 m/s. c Changes in the resonant frequencies with the pipe radius ratio, assuming the depth of the magma surface is 200 m and a o is fixed. The frequencies (f 0 and f 1 ) are normalized by f 0cylinder , which is the fundamental resonant frequency of the cylinder closed on one side and open on the other (a o /a c = 1) this study, however, a simple model was adopted to decrease the number of unknowns. Ishii and Yokoo (2021) assumed that the seismic wave first propagates at a very slow velocity, v p0 , through a lava pond and then travels through the ground at P-wave velocity, v p (Fig. 3a). The observed distribution of t can now be explained by considering the horizontal change in the explosion source position with the assumption of such a slow v p0 (Ishii and Yokoo 2021). Although we cannot discount the possibility that there may be other factors causing variation in t (Supplementary information 1), the model of Ishii and Yokoo (2021) is used here for the sake of simplicity. The infrasound signal propagates inside the vent at a sound velocity c modified by the temperature of the gas and the ash content, outside of the conduit the sound velocity is that in the air c air (Ripepe et al. 2001). The time delay between the arrivals of seismic and infrasound signals, t, can be now also expressed as a function of d and c: Here, r is the distance from the center of the magma surface to the explosion source (0 ≤ r ≤ a c ), L i is the infrasound path length from the vent to the station, and L s is the seismic path length from the margin of the magma surface to the station (Fig. 3a). P-wave velocity v p is assumed to be 2.7 km/s (Tsutsui et al. 2003), and c air is obtained using the temperature measured at the JMA weather station (Fig. 1a). The maximum value of v p0 can be constrained so that the variation in the observed distribution of t is explained as the variation in the horizontal location of the explosion source in the range of 0≤ r ≤ a c (see Ishii & Yokoo, 2021 for details). Therefore, we used the maximum value of v p0 for Eq. 1.

Classification of seismic and infrasound waveforms
We detected 107,988 events by the STA/LTA method, and then classified them into several waveform types in order to automatically read the arrival time of signals (Supplementary information 2). Seismic waveforms associated with explosions could be classified into two types based on their waveforms. All seismic signals had their primary energies in the 2-15-Hz range. However, some events were predominantly composed of high frequencies (>10 Hz) and others were predominantly composed of low frequencies (<10 Hz) ( Fig. 4a, b). The slope of the regression line of the frequency spectrum can thus be used to classify seismic waveforms into two types: the first type with a positive slope (SH; Highfrequency) and the second type with a negative slope (SL; Low-frequency).
Infrasound waveforms associated with explosions observed at different times during the 2014-2015 eruption were distinctly different ( Fig. 2e-h). The explosion-related infrasound waveforms could be classified into six types (Fig. 4c). The 2-, 4-, and 15-Hz components of the normalized spectrum by the 0.5-Hz component were used as features for the classification (Fig. 4d). When these three features were plotted in 3-D space, each type was distributed in a different coordinate position (Supplementary information 2). A representative position was defined for each type, and all events were classified into six types based on the smallest distance from the representative positions (Supplementary information 2). Figure 5 shows the result of the classification. Although the major type of seismic signal recorded before April 2015 was type SH, type SL was dominant in April and until the crater collapse. The type of infrasound signals transitioned as time passes. Multiple types coexisted in December 2014, and then the major type changed progressively from IN to IH, then to Ih, and finally to IL between January and April 2015. The ratio of type IN temporarily increased in late-April 2015. After the crater collapses, type IA emerged in May 2015. However, the number of events decreased significantly after the crater collapse on May 3 (Fig. 5a).
The transition of the infrasound type can be related to the change in eruptive style. Types IN and IM were mainly observed when explosions occurred during the period of continuous ash emission during December 2014 and late-April 2015. Types IH and IL were observed when levels of ash emission were relatively low during February and mid-April 2015. Type Ih was generated from gas-dominant . They are mainly composed of high-frequency waves (>10 Hz) and low-frequency waves (<10 Hz). b Amplitude spectra of the waveforms shown in a. They are calculated from 3-s waveforms (indicated by colored lines in a and normalized with averages of 10-20 Hz. Colored lines are regression lines of 3-15 Hz (the red line is SH and the blue line is SL). c Infrasound waveforms of the six types (IL, IM, IH, IN, Ih, and IA). d Average amplitude spectra of events classified into each type of infrasound signal. These spectra are calculated from 7-s waveforms and normalized with 0.5-Hz components. Clear peaks are recognized at 0.5, 2, 4, and 15 Hz (hatched areas). Type IL, IM, and IH have more distinct peaks at 2, 4, and 15 Hz, respectively. Type IN is the most common type. Type Ih is characterized by a relatively large 0.5 Hz and small high frequencies (especially at 15 Hz). Type IA appears to be a shock-like signal at the initial phase. This type has no distinct peaks in its spectrum, unlike the other types. The amplitudes of the 2-, 4-, and 15-Hz components are used as features for the classification (Supplementary information 2) explosions without ballistics during mid-March; and type IA was associated with gas leakage from the collapsed crater floor accompanied by mud explosions. On the other hand, change in seismic waveform type did not clearly relate to explosive style.

Temporal change in the crater geometry and the shape of the conduit pipe
The size of the active vent was constrained by field observations , where the vent diameter rapidly increased until early-January 2015 (Fig. 6a). In addition, a pyroclastic cone developed in the crater during the same period. Based on these diameter changes, the analysis period was divided into five periods (I-V; Fig. 6a), and elevation data for the 3-D numerical simulation of infrasound propagation (Ishii and Yokoo 2021) were created for each period. The elevation data obtained using drone photogrammetry in March was used for period V (Figs. 1f and 6a). We assumed that the elevation profiles of the pyroclastic cone in periods I-IV have similar geometries to that in period V (Fig. 6b, c). The profile of each period was determined from the cone's rim diameter to maintain the aspect ratio of the cone of period V. The elevation profiles of the inner side of the cone's rim were interpolated to meet the observed vent size, and a conical frustum pipe was set inside the vent.
The radius ratio of the conical frustum was determined by the ratio of resonant frequencies, the fundamental frequency (f 0 ), and the overtone component (f 1 ) of infrasound signals. Although f 0 changed during the 2014-2015 eruption, f 1 was stable at ∼2 Hz (Fig. 7a). The probability density The percentage of each type of waveform function (PDF) of f 1 /f 0 was thus calculated every day using the amplitude spectra for data recorded every 10 min ( Fig. 7b; Ishii and Yokoo 2021), and we chose three or four radius ratios so that the sum of the PDF values of these shapes exceeded 0.75. The depths of the magma surface were then estimated assuming these radius ratios, and then a weighted mean of these results was taken (Ishii and Yokoo 2021). The observed vent diameter and the radius ratio provided the temporal change in a o and a c (Fig. 7c), where a o is approximated by the radius of the vent (Fig. 6a). From the beginning of the eruption to January, a o /a c was close to one, indicating that the space above the magma surface was a cylindrical shape. This period coincided with the period when the diameter of the vent and the pyroclastic cone rapidly grew . After this point in time, a c became larger than a o , meaning that the vent broadened so that the upper conduit became a conical frustum whose radius widens toward the inside.

Depth estimation
We obtained the relationship for the peak frequency (f 0 = F(d, c)) from the 3-D numerical simulation of infrasound propagation using the time-varying elevation data, and then estimated the depth of the magma surface by combining it with the relationship for the time delay ( t = G(d, c); Eq. 1). The peak frequency of the infrasound signal (f 0 ) was calculated from a 40-s long waveform (-10-30 s from the arrival time) for each event. The arrival times of the infrasound and seismic signals were automatically determined using the STA/LTA method again, with shorter time windows (0.1 s and 1 s) and different thresholds for each waveform type (Supplementary information 2). Figure 8b and c show the temporal changes in the distributions of f 0 and t. The peak in the probability density function for f 0 increased from 0.5 to 0.7 Hz until January, but then decreased to 0.6 Hz in February. It was stable at 0.6 Hz until early-March, but dropped to ∼0.4 Hz in late-March, and increased again to ∼0.5 Hz in April. The distribution of t had a peak at ∼0.5 s with no significant fluctuation from the beginning of the eruption through to mid-March. After mid-March, the peak decreased to 0 s; but increased to 0.7 s in April. These trends were independent of the waveform types of the seismic and infrasound signals. The distributions of both f 0 and t for days when the daily number of explosion events was <210 are not included in the depth estimation analysis because the distributions of both f 0 and t had no clear peaks (Fig. 8a).
The depth of the magma surface in the conduit was estimated every day from the distributions of f 0 and t (Ishii and Yokoo 2021). Figure 9a shows examples of estimated depths from given pairs of f 0 and t, assuming that the explosion source is located at distance r from the center of the magma surface. The distributions of f 0 and t The representative diameter is determined for each period (horizontal lines) to create the elevation data. A pink line represents the date when the drone photogrammetric survey was conducted. b Cross-section of the crater topography (east-west direction). The horizontal axis shows the distance from the center of the vent (vertical dashed line). Although the data of periods I-IV are created based on the representative diameters (Fig. a), only the data of period V is measured by a drone in March 2015. The horizontal dotted line represents the crater floor before the eruptions (1140 m a.s.l.). The horizontal bars above the crater indicate the range of the rim of the pyroclastic cone and exit of the vent of periods I and V. The lower bar (exit of the vent) corresponds to twice a o (Fig. 3a). c Contour maps of the active crater for each period were indicated by gray contours in Fig. 9a, which are the distribution of April 25 as an example. The depth indicated by the color in the area that overlaps with the gray contours in Fig. 9a is thus the estimated value, and gives a range of 60-200 m (Fig. 9c). Following Ishii and Yokoo (2021), this method can simultaneously constrain the sound velocity in the conduit (Fig. 9b), for which we obtain 300-790 m/s (Fig. 9c). Magma surface depth evolved with time ( Fig. 9d). At the beginning of the eruption in December 2014, the depth was located at ∼200 m. The magma surface then rose to ∼120 m by late-January 2015, and remained at 120-140 m until late-April. The magma surface then dropped by ∼70 m before the crater collapsed in May (Fig. 1d). During this time, though, sound velocity did not significantly change (Fig. 9e).
In our model, the seismic signal propagates at a slow velocity (v p0 ) in the lava pond (Fig. 3a). Equation 1 shows that the value of t depends on v p0 . We can thus constrain the maximum v p0 that can account for the observed variation in t (Fig. 8c) when the horizontal position of the source changes within the range of 0 ≤ r ≤ a c at the magma surface following Ishii and Yokoo (2021). For example, v p0 on April 25 was determined to be 35 m/s. The constrained v p0 for every day during the analysis period was used for Eq. 1 (Fig. 8d). Although the maximum v p0 was <10 m/s at the beginning of the eruption, it increased from mid-January to exceed 30 m/s. Then, it decreased to 20 m/s in mid-March and increased again in April to 30 m/s. In general, the seismic velocity of the magma depends on the volume fraction of gas (Morrissey and Chouet 2001). The v p0 change shown in Fig. 8d may reflect a temporal change in the amount of volatiles contained in the upper part of the magma column. The low velocities in December and mid-March suggest that the magma contained more volatile components at these points in time.

Discussion
The magma surface rose between January and February 2015, and dropped before the crater collapsed (Fig. 9d). At the same time, the space above the magma surface changed from a cylindrical shape to a frustum whose radius widens toward the inside (i.e., flaring inside; Fig. 7). We first compare our estimation of magma surface depth with other observed records and discuss the relationship among them. Based on the temporal variation in activity and waveform, we can divide the eruption into six periods (P1-P6; Fig. 10). The waveform types and surface phenomena for each period are summarized in Table 1. The magmatic eruption comprises P1-P5, and the collapse of the crater floor that coincided with cessation of activity occurred at the end of P5, so that P6 is a post-eruptive phase.

Temporal change in the depth of the magma surface and volcanic activity
At the beginning of the eruption (P1), the magma surface was located at ∼200 m (Fig. 10a), which was the deepest level of any period. P1 was characterized by a variety of infrasound waveform types (Fig. 10a), and the tephra discharge rate fluctuated ( Fig. 10d; Miyabuchi and Hara 2019). This suggests that the magma supply had not yet stabilized the magma column in the conduit. In P2, the magma surface gradually rose from 162 to 129 m during January 10-27 (Fig. 10a). Prior to this change, strainmeters and tiltmeters at HDT station included by our observation network (Kyoto University) recorded deformation (Fig. 10b). For example, the tilt vector on January 5-8 pointed toward an area south of the active . The magma surface rose from ∼200 to ∼120 m in January and dropped ∼70 m in May (black arrows). e Temporal change in the representative sound velocity crater (Fig. 11a, d). Assuming a spherical inflation source (Mogi 1958) located at the intersection of the extension line of the tilt vector and a crack-like conduit beneath the craters ( Fig. 11d; Yamamoto et al. 1999), the depth of the source is estimated to be 1.5 km using the ratio of tilt change and the area strain (cf. Ishihara 1990). The volume change in the pressure source is 6.5 × 10 3 m 3 . The intrusion of new magma into the shallow part of the volcano is thus inferred at this time. We can interpret the rise of the magma surface in P2 to be related to this magma intrusion. In addition, Miyabuchi and Hara (2019) reported that the proportion of pale brown glass included in the ash samples increased in mid-January. The pale brown glass is generally created by rapid cooling of magma (Miyabuchi and Hara 2019). Therefore, this change in the ash component is again consistent with an injection of new magma at this time. The increased magma supply derived from magma injection may provide a high tephra discharge rate (Fig. 10d) and large seismic amplitudes (Fig. 2a). The major type of infrasound waveform becomes IN in P2, which can be related to the continuous ash emission (Fig. 10a).
The rise in the magma surface stopped in P3 (Fig. 10a). At the same time, the tephra discharge rate decreased (Fig. 10d). P3 thus likely marks a decline in magma supply rate to the shallow conduit. In this period, the baseline length between GNSS sites across the magma chamber, which is located at a depth of 4-6 km (Sudo and Kong 2001;Sudo et al. 2006) also stagnated and then began to contract  Fig. 5a). b Strain and tilt data recorded at HDT station (Fig. 1a). Radial and tangential components with respect to the active vent are represented. The positive direction is extension or up. c Baseline length between two GNSS stations (HND-960701; Fig. 11d) across the magma chamber (gray dots) and the running mean (31 days; black line). Red markers are SO 2 discharge rates (provided by the Japan Meteorological Agency). d Tephra discharge rate (Miyabuchi and Hara 2019) and precipitation measured at station JMA (Fig. 1a ( Fig. 10c). This suggests that the supply to the magma chamber began to decrease. The major type of infrasound waveforms changed from IN to IH (Fig. 10a). This change can be linked to the decrease in ash emission, which in turn can be related to changes in the dynamics of the explosions and/or the density of the propagation medium (Meier et al. 2016;Johnson and Aster 2005).
The surface phenomena changed drastically in P4 (Table 1). Tephra discharge rate became low (2 × 10 3 t/d; Miyabuchi and Hara 2019), and ejections of ballistics were also rare. On the other hand, SO 2 flux increased (Figs. 2c and 10c). At night, we observed continuous glowing inside the vent through P4. These observations imply that the degassing rate increased, probably because volatile-rich magma was injected into the shallow region. The low values of the seismic velocity in the lava pond (v p0 ; Fig. 8d) in P4 also suggest that the magma contained a large number of bubbles, consistent with an injection of volatile-rich magma. Miyabuchi and Hara (2019) pointed out that the magma also became less viscous due to changes in magma temperature, magma composition, and/or water content from late-March, which would have played a role in changing the emission style (cf. Gurioli et al. 2014;Leduc et al. 2015). Although the magma surface is estimated to remain at a shallow level (120-140 m) in P4 (Fig. 10a), a slight rise in the magma surface at the end of P3 (by ∼20 m during March 19-20) may reflect magma injection (Fig. 10a). Strain and tilt also changed before and at the beginning of P4 (Fig. 11b). The tilt vectors began to point to a more southern area than the series of craters (Fig. 11d). This change in direction may be due to rainfall events occurring just before P4 when 123 mm of rain fell during March 18-19, which is equivalent to 58% of the total rain in March 2015 (Fig. 10d). The ground deformation by increasing the amount of underground water around tiltmeters may overlap the volcanic deformation, resulting in an apparent change in the tilt direction. The infrasound signal continued to oscillate with ∼50 Pa during P4 (Fig. 11b) even if no explosion occurs. It was occasionally amplified up to ∼200 Pa (Fig. 2g), which was detected and classified as type Ih (Fig. 4c). This suggests that high gas fluxes promoted passive degassing (Burton et al. 2007), which continuously excited acoustic resonance in the conduit. Another possibility is that the flow pattern in the magma column changed due to the increase in the gas flux (Valade et al. 2018). The infrasound waveform of type Ih can be related to a ballistic-poor explosion (like Type 0; Leduc et al. 2015), which can strongly amplify the resonance. The volatile-rich nature of the magma surface at this time is consistent with the shallow conduit conditions when Type 0 occurs at Stromboli (Leduc et al. 2015;Gaudin et al. 2017). After P4, explosions with ballistics resumed, and the tephra discharge rate increased in mid-April ( Fig. 10d; Table 1). After April 26, the occurrence of explosions and infrasound amplitude decreased (Fig. 11c), and continuous ash emission became predominant. Such variation in surface phenomena likely reflects the process of declining volatile content in the magma. Infrasound signals of type IL, which have a distinct 2-Hz component, increased in P5, implying that overtones were strongly excited because the resonance-exciting sources (pressure perturbations inside the vent) changed to higher frequency components. The type of seismic signals also changed from SH to SL. However, the cause of this change in the seismic nature is not clear. Since low-frequency components (3-5 Hz) of background tremor are amplified in this period, it is possible that superposition with the amplified low-frequency tremor changes the spectral characteristics of the seismic events. Although the depth of the magma surface did not change through most of P5, the magma surface dropped by ∼70 m beginning on April 26 (Fig. 10a). This suggests an abrupt decrease in magma supply at the end of P5. This decline in magma surface level was followed by the crater floor collapse on May 3. The volume of the space above the magma surface at the end of P5 is estimated to be ∼7 × 10 5 m 3 . This volume is consistent with the volume of the southeastern subsidence area calculated from the elevation profile before and after the crater collapsed (4.5 × 10 5 m 3 ; Fig. 1f, g). Therefore, this collapse would have filled more than half of the space below the vent created by evacuation of the shallow system. The center of the collapsed area is south of the active vent, implying that the shallow conduit was inclined to the south (Fig. 1d, g). The source locations of the continuous tremor before the magmatic eruption (during December 2013-January 2014) are also estimated to be south of the center of the crater at <400-m depth (Ichimura et al. 2018). The plumbing system thus likely extended to the south, as also indicated by the tilt vector in Fig. 11d(1).
Although the depth of the magma surface could not be estimated in P6, changes in ground deformation and waveform of the infrasound signals were observed. Both are considered to be caused by blockage of the conduit due to the crater collapse. Immediately after the collapse, although tilt changes were very small, strain changes were observed (Fig. 11c). This included contraction in the radial direction and extension in the tangential direction. These changes are similar to the strain records observed before the explosions  Yamamoto et al. 1999). A gray area represents the approximate position of the magma chamber inferred from low P-velocity and leveling data (4-6-km depth; Sudo and Kong 2001;Sudo et al. 2006). The black star represents the epicenter of the earthquake on May 8. The black solid line represents the GNSS baseline across the magma chamber (Fig. 10c). HND is one of the GNSS stations of our observation network (Kyoto University) and the other is the 960701 GEONET site (130.9962 • E, 32.8707 • N) of Showa crater of Sakurajima (Iguchi et al. 2013), and can be explained by the spherical inflation of a source near the ground surface. Using the ratio of radial and tangential strains (cf. Iguchi et al. 2013), and assuming that the Mogi source is at the same horizontal position as before P2 (Fig. 11d), the depth is estimated to be ∼200 m. The volume change in the source is 2.0 ×10 2 m 3 , which is comparable to, for example, the shallow inflation before Vulcanian explosions of Suwanosejima and Sumeru volcanoes (1-5 × 10 2 m 3 ; Iguchi et al. 2008) and the local inflation before the 2018 phreatic eruption of Iwo-Yama volcano (∼10 2 m 3 ; Narita et al. 2020). This very small inflation can be interpreted as a temporary increase in pressure beneath the crater due to blockage of the conduit. After the crater collapsed, the waveform type of the infrasound signal changed to IA, which was characterized by a shocklike waveform (Figs. 10a and 4c). These shock-like signals were simultaneous with mud explosions inside the crater, as that impulsive gas leakage repeatedly occurred through the collapsed crater floor. The characteristic strain changes became ambiguous after a local earthquake that occurred on May 8 (Fig. 11c). The magnitude of this earthquake was 2.7, and its source was located ∼600 m south of the crater at ∼1km depth ( Fig. 11d; Japan Meteorological Agency 2015). The earthquake may have been caused by the pressure release within the plumbing system extending to the area south of the crater. Thereafter, the GNSS baseline length across the magma chamber turned to extension (Fig. 10c), possibly heralding a new phase of magma accumulation in the magma chamber.

Pipe shape of the conduit
One of differences between this study and other studies in terms of depth estimation is the consideration of temporal variation in the shape of the conduit (cf. Richardson et al. 2014;Spina et al. 2015;Johnson et al. 2018a). As the pipe shape changes, the relationship between the pipe length and resonant frequency also changes (Watson et al. 2019(Watson et al. , 2020. Figure 12a compares magma surface depths estimated when the pipe shape is fixed and varied. When the shape is fixed as a cylinder, the magma surface slightly drops during P2 because the fundamental resonant frequency (f 0 ) decreases. In contrast, the magma surface rises in the case that the shape of the conduit is assumed to change based on the temporal variation in f 1 /f 0 (Fig. 7). Our estimated depth during P2 is consequently more consistent with the ground deformation in early-January as we consider the evolution in conduit shape (Fig. 11a). Likewise, our estimated depths for P3-P5 are shallower than those assuming a fixed shape. Therefore, we emphasize that considering temporal changes in the shape of the conduit plays a vital role in the depth estimation using the resonant frequency. The recurrent monitoring of the local topography of an active vent is essential in constraining conduit shape.
We note, though, that our method to constrain the shape of the conduit involves a controversial point. That is, we used the analytical relationship of Ayers et al. (1985) to constrain the shape (Fig. 3b), a selection based on the low computational cost of this method. However, this relationship is not strictly appropriate for volcanic Fig. 12 Depths estimated by different conduit shapes. a Depths estimated by assuming a cylindrical conduit. Red circles are representative depths and error bars indicate the variation in the estimated value. Black circles are estimated depths considering the temporal change in the shape of the conduit pipe (same as Fig. 9d). b Temporal change in the radii of the conduit pipe (same as Fig. 7c). The depths indicated by black circles in Fig. a are estimated considering this radius change scenarios, especially in terms of boundary conditions. Figure 3b shows the relationships between f 1 /f 0 and a o /a c computed by the numerical simulation. Note that, the radius ratio based on Ayers's relationship is slightly smaller than the radius ratios based on the computed relationship. In addition, the computed relationship depends on the sound velocity inside the conduit. This means that the analytical relationship of Ayers et al. (1985) can result in an exaggeration of conduit shape (more flaring) and an underestimated depth for the magma surface. Although it does not affect the trend of temporal changes in shape and depth in this study, iterative calculations of a lowcost numerical simulation (Watson et al. 2019) are more appropriate to rigorously constrain the conduit shape to provide a more accurate value for the absolute depth.
The constraint of the temporal change of the conduit shape also aids in developing insights into the activity. The frequency ratio (f 1 /f 0 ) indicates that the pipe shape changed from a cylinder to a conical frustum flaring inside (Fig. 7). Although such a shape is thought to be physically unstable (cf. Orr et al. 2013, Eychenne et al. 2015, a flaring space inside the conduit has also been observed at Kilauea and Villarrica (Orr et al. 2013;Palma et al. 2008). Such a form is the result of rock erosion due to the heat of the magma column (Fee et al. 2010) or the adhesion of spatter to the inner wall of the vent (Palma et al. 2008). In our case, the shape change occurred during the same period as the rise in the magma surface (Fig. 12). This may suggest that heat from the rising magma surface made the inner wall of the conduit fragile. Consequently, the wall fell especially near the magma surface to widen the conduit. The period when the conduit is similar to a cylinder (until early-January) matches the period when the vent size rapidly grows ( Fig. 6a; Yokoo et al. 2019). During this period, the conduit may have been abraded by the ash and gas emissions while preserving a cylindrical shape because the conduit was relatively narrow.
The physically unstable shape of the flared conduit, coupled with magma withdrawal, probably led to crater floor collapse at the end of the eruption. At the end of P5, the depth of the magma surface dropped by ∼70 m, which is considered to be one of the triggers of crater collapse. However, the amount of this drop was not great and caused the magma level to return to that of P1. During P1, the shape was cylindrical and stable. This suggests that the physical instability of the conduit would have allowed the crater floor to collapse even with a small change in magma level. Therefore, the eruption may have been forcibly stopped by the collapse. The shallow inflating deformation after the crater collapse (May 4-8, Fig. 11c) implies that the magma supply continued so as to increase the pressure beneath the crater.
The change in the fundamental resonant frequency (f 0 ) suggests that the magma drained during April 25-27 or 28. Considering a 70-m drop during this period, the rate of drainage is 23-35 m/day or 1.8-2.7 × 10 5 m 3 /day (assuming a 50-m radius conduit). This rate is much lower than that of the lava lake withdrawal at Kilauea in 2018, 40-53 m/day (i.e., 1.3-1.7 × 10 6 m 3 /day assuming a 100m radius; Anderson et al. 2019). This difference in the rate may be due to the difference in the mechanism of magma drainage. The drainage at Aso volcano can be derived from a decrease in magma supply from the magma chamber to the shallow conduit, which can be interpreted as a termination process of eruptive activity. Instead, drainage at Kilauea was caused by the formation of a new magma pathway toward the volcano flank (the intrusion of a dike) in a phase where volcanic activity was getting higher ). This dike intrusion was likely caused by pressurization of the plumbing system, which is corroborated by the rise of magma levels at Pu'u O'o and the summit lava lakes prior to the dike intrusion (Neal et al. 2019). The drops of magma level observed at Stromboli (Ripepe et al. 2005(Ripepe et al. , 2015 and Etna (Spina et al. 2015;Cannavo' et al. 2019) were related to lava flow or dike intrusion at the flank when the activities were getting higher, thus were similar to the case of Kilauea. In these cases, the rate of magma drainage from the summit (and old) conduit is mainly determined by the magma flow rate into a new branch of the conduit and/or the magma discharge rate from a new vent. These magma flows are driven by gravity and/or the pressure of the plumbing system (Ripepe et al. 2015;Anderson et al. 2019). In Stromboli, the gravity-driven model (Ripepe et al. 2015) can explain the observed effusion rate (∼1.0 × 10 6 m 3 /day at the beginning of the effusive event in 2007), which is also a higher rate than the drainage rate of Aso as well as that of Kilauea. Given there is no evidence to support the formation of a new magma pathway or vent at Aso at that time, the drop in magma level at Aso should be different from these examples, resulting in a lower drainage rate. Therefore, it may be better to distinguish between the drainage rate in the process of terminating eruptive activity and that associated with flank eruptions.
The observed frequency ratio (f 1 /f 0 ) slightly decreases during the drop in the magma surface, which means that the shape of the resonator slightly reverts to like a cylinder (Figs. 7 and 12). It is thus possible that the drop in the magma surface expanded the space acting as the acoustic resonator and made the space approximate a shape close to a cylinder.

The 2014-2015 eruption of Aso volcano: a conceptual model
Based on our findings, we can summarize the activity of the 2014-2015 eruption and propose a conceptual model, as given in Fig. 13. This model progresses through six periods of activity (P1-P6) as follows.

P1:
The magma surface remained at relatively deep level (∼200 m) at the beginning of the eruption. The conduit pipe was narrow and cylindrical, but then expanded while preserving its shape during continuous ash emission. P2: The magma surface rose due to magma intrusion at depths of 1-2 km. The heat from the magma surface weakened the conduit wall, and the conduit above the magma surface widened to gradually take on the form of a conical frustum flaring inside. At the beginning of the eruptions, the magma surface was located in a deep region (P1). The conduit was a cylindrical shape. The diameter of the vent rapidly increased while maintaining the cylindrical shape. As the magma supply increased, the magma surface rose, and the conduit pipe became a frustum flaring inside (P2). In P3, the magma surface stops rising, and the tephra discharge rate decreases due to the decline in the magma supply. Although gas emissions occurred predominantly in P4, probably because of the additional injection of volatile-rich magma, the eruptive style changed into ash emissions in P5. At the end of P5, the magma surface dropped, and the crater floor collapsed. The space above the magma surface was filled with collapsed materials. Thus, the pressure increased beneath the crater floor (P6)

P3:
The magma supply decreased. The magma surface stopped rising, and the tephra discharge rate also decreased. The decline in magma supply (and ash emission) resulted in an evolution of the infrasound waveforms. P4: The gas flux increased due to the additional intrusion of volatile-rich magma. Continuous passive degassing effectively excited the acoustic resonance, and the amplitude of the infrasound signal became the largest of the 2014-2015 activity. Ballistic-poor explosions often occurred. P5: The flux of gas emissions decreased. The eruptive style transitioned into ash emissions. A slight decline in the magma supply, and hence magma level, coupled with the unstable shape of the conduit induced crater collapse. The collapse blocked the upper conduit and stopped the magmatic eruption. P6: The pressure beneath the crater floor increased due to blockage of the conduit. Signals of infrasound pulses were observed when gas escaped through the collapsed crater floor.

Conclusions
We tracked the temporal change in the depth of the magma surface and the evolution of shallow conduit shape using seismo-acoustic signals for the 2014-2015 eruption of Aso volcano. In terms of eruption mechanisms and controls, the primary conclusions of this study are as follows: (1) The depth of the magma surface rose to a level of ∼120 m below the crater floor between January and February 2015. Magma intruded at a 1-2-km depth in early-January induced the rise of the magma surface. The shape of the conduit changed from a cylinder to a conical frustum due to wall-heating during the rise in the magma surface. An additional magma intrusion in mid-March also occurred. A clear change in depth was not seen, but due to the volatile-rich nature of the magma, a change in P-wave velocity through the lava pond was suggested.
(2) A slight decline in the magma supply leading to a small decrease in magma level, and the unstable shape of the conduit (conical frustum flaring inside), caused the crater floor to collapse. This stopped the magmatic eruption by plugging the upper conduit. This finding implies that the termination process of the eruption can be affected not only by a decrease in magma supply but also by blockage of the conduit by collapse of a crater floor.
Monitoring the depth of the magma surface and the topography around the active vent can aid in understanding the evolution of, and termination mechanisms for, magmatic eruptions. We stress, though, that this can only be done through the simultaneous constraint of conduit shape and magma surface depth. However, we present here a method and blueprint as to how this can be done using one seismic station and one infrasound station, so as to allow near-real time assessment and tracking of conduit conditions.
Acknowledgements KI acknowledges constructive comments from M. Iguchi, H. Nakamichi, K. Kuge, and Y. Miyabuchi. We thank all staff of the Aso Volcanological Laboratory for their support. We also thank L. Watson and an anonymous reviewer whose comments helped improve our study. Editors, M. Ichihara and A. Harris, gave us helpful feedback. SO 2 flux and weather station data were provided by the Japan Meteorological Agency. GNSS data of 960701 GEONET site were provided by the Geospatial Information Authority of Japan.
Author contribution KI performed all analyses in this study. KI and AY contributed to the planning and interpretation of this study. TO acquired geodetic data of the underground tunnel and GNSS and provided insights into the ground deformations. TK and TO made the DEM data of the crater. All authors read and approved the final manuscript.
Funding The present study was financially supported by the Ministry of Education, Culture, Sports, Science and Technology of Japan under its Earthquake and Volcano Hazards Observation and Research Program, JSPS KAKENHI Grant Number JP20J11915 to KI, and ISHIZUE 2022 of Kyoto University to AY.

Data availability
The seismo-acoustic data in this study are available from the corresponding author upon request.

Competing interests
The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.