Tilt and strain change during the explosion at Minami-dake, Sakurajima, on November 13, 2017

We herein propose an alternative model for deformation caused by an eruption at Sakurajima, which has been previously interpreted as being due to a Mogi-type spherical point source beneath Minami-dake. On November 13, 2017, a large explosion with a plume height of 4200 m occurred at Minami-dake. During the 3 min following the onset of the explosion (November 13, 2017, 22:07–22:10 (Japan standard time (UTC + 9); the same hereinafter), phase 1, a large strain with changes up to 120 nstrain was detected at the Arimura observation tunnel (AVOT) located approximately 2.1 km southeast from the Minami-dake crater. After the peak of the explosion (November 13, 2017, 22:10–24:00), phase 2, a large deflation was detected at every monitoring station due to the continuous Strombolian eruption. Subsidence toward Minami-dake was detected at five out of six stations, whereas subsidence toward the north of Sakurajima was detected at the newly installed Komen observation tunnel (KMT), located approximately 4.0 km northeast from the Minami-dake crater. The large strain change at AVOT as well as small tilt changes at all stations and small strain changes at the Harutayama observation tunnel (HVOT) and KMT during phase 1 can be explained by a very shallow deflation source beneath Minami-dake at 0.1 km below sea level (bsl). For phase 2, a deeper deflation source beneath Minami-dake at a depth of 3.3 km bsl was found in addition to the shallow source beneath Minami-dake, which turned inflation after the deflation that occurred during phase 1. However, this model cannot explain the tilt change of KMT. Adding a spherical deflation source beneath Kita-dake at a depth of 3.2 km bsl can explain the tilt and strain change at KMT and the other stations. The Kita-dake source was also found in a previous study of long-term ground deformation. Not only the deeper Minami-dake source MD, but also the Kita-dake source deflated due to the Minami-dake explosion.


Introduction
The tilt and strain changes accompanying Vulcanian eruptions of Sakurajima have been interpreted as being due to the deflation of a Mogi-type source (Mogi 1958) beneath Minami-dake, because downward tilts associated with eruptions are oriented towards Minami-dake in the underground tunnels Harutayama (2.7 km northwest of the Minami-dake crater) and Arimura (2.1 km south-southeast of the crater) ( Fig. 1). For the eruptions at Minami-dake starting in 1955, a deflation source was found beneath Minami-dake at a depth of 2-4 km below sea level (bsl), with a deflation volume of 10 3 -10 5 m 3 . This source was also found during an episode of frequent Vulcanian eruptions at Minami-dake that peaked from 1974 to 1992 by precise leveling surveys (Eto 1989). A shallower deflation source was also found at the onset of the eruptions (e.g., Ishihara 1990;Iguchi et al. 2013). For the eruptions at Showa crater starting in 2006, only a shallower deflation source was found beneath Minami-dake at a depth of 0-1.5 km bsl, with a deflation volume of 10 2 -10 4 m 3 , one-order smaller than that Open Access *Correspondence: hotta@sus.u-toyama.ac.jp 1 Faculty of Sustainable Design, University of Toyama, Toyama 930-8555, Japan Full list of author information is available at the end of the article of the Minami-dake eruptions . During an episode of eruptions at the Showa crater, a different source was found beneath Kita-dake at a depth of 3.3 km bsl by continuous GPS observation data. This source was considered to exhibit repeating inflation and deflation with a period of several years (Hotta et al. 2016a). In August 2016, additional tilt and strain observations were undertaken at a new observation tunnel at Komen, 4.0 km northeast of Minami-dake, and changes were detected, associated with a large explosion with a plume height of 4200 m at Minami-dake on November 13, 2017. More than 1 year had passed since installation of tiltmeters and strainmeters in Komen observation tunnel when the eruption occurred. Ground inflation started at 8:00 (Japan standard time (UTC+9); the same hereinafter) on November 7, 2017, which became a gradual deflation on November 11, 2017. When an explosion occurred at 22:07 of November 13, 2017, the gradual deflation became rapid and continued for two hours. The total duration of tilt/strain change of nearly a week is longer than the typical eruptions at Minami-dake in the 1980s of at most 7 h (Ishihara 1990). The downward tilt direction at the Komen tunnel deviated by 90° clockwise from the direction of the Minami-dake crater, in a different manner than at the Arimura and Harutayama tunnels, which showed a downward tilt in the direction of the crater. This tilt direction is substantially different from the deformation expected from the inflation of the Mogi-type spherical point source beneath Minami-dake found in previous studies, and thus the deformation at the Komen tunnel is considered to be affected by a different source.
In the present study, we propose a composite source model for tilt and strain changes at ground deformation stations, taking into account the tilt change at the Komen tunnel. We then discuss the magma plumbing system of Sakurajima for a large-scale Vulcanian eruption.

Observations and data
At Sakurajima, tilt and strain observations have been performed by the Sakurajima Volcano Research Center (SVRC) of Kyoto University; the Osumi Office of River and National Development Bureau, Ministry of Land, Infrastructure, Transport and Tourism of Japan (MLIT); and the Japan Meteorological Agency (JMA). The tilt and strain stations for the present study are distributed as shown in Fig. 1. We used data obtained by watertube tiltmeters and linear strainmeters at the Arimura, Harutayama and Komen underground tunnels (AVOT, HVOT and KMT, respectively) and borehole tiltmeters at Yokoyama, Amida-gawa and Komen (JMA-A, JMA-F and KOM, respectively). The tiltmeters and strainmeters at AVOT, HVOT, KMT and KOM were installed in directions radial and tangential to Minami-dake and the tiltmeters at JMA-A and F were in the north-south and east-west directions. The sensitivities were 0.056 V/μrad for the water-tube tiltmeters and 0.056 V/μstrain for the linear strainmeters. The sensitivity of the borehole tiltmeters was 0.02 V/μrad. The analog-to-digital resolution was 24 bits and the resolution of tilt and strain was 10 −6 nrad.
The tilt and strain changes at AVOT, HVOT and KMT during the period of November 5-14, 2017, are shown in Fig. 2. Upward tilts in the direction of Minami-dake and extension strains started to appear at about 6:00 of November 7 at AVOT and HVOT, whereas a tangential tilt was dominant at KMT. The tilts turned to gradual downward tilts on November 11, 2017. Then, the downward tilt suddenly became rapid at the onset of the explosion on November 13. A deflation tilt change of 283 nrad occurred over two hours at HVOT, which is larger than the tilt changes associated with explosions at Minamidake crater in the previous reports (Ishihara 1990;Tateo and Iguchi 2009).
The tilt and strain changes associated with the explosion are magnified in Fig. 3. In phase 1 (3 min from the onset of the explosion at 22:07), rapid extensions were detected in the radial strains as strain steps. The step at AVOT, the nearest station to the Minami-dake crater, is the largest at 120 nstrain, and the steps are smaller at

Deformation source model
We applied a finite spherical source (McTigue 1987), a dislocation source (Okada 1992) and a combination of these in a half-infinite homogenous elastic solid. We assumed a Poisson's coefficient of ν = 0.25 . Battaglia et al. (2013) obtained the volume change of a spherical source V by an internal pressure change of P , for a modulus of rigidity µ and depth and radius of the source d and a , respectively, to be From Eq. (1), ΔP/μ in the equations of deformation by McTigue (1987) can be obtained using ΔV. First-order topographical effects were corrected by adding the elevation of the station to the source depth (Williams and Wadge 1998). The effect of the topography using the finite element method (FEM) is discussed later.

Fig. 2
Tilt and strain records at the AVOT, HVOT and KMT stations during the period from November 5 through 14, 2017. δR and δT represent radial and tangential tilt changes, respectively. εR, εT and εO are strain changes in the radial, tangential and oblique components, respectively. The increase in the tilt represents an upward change in a direction toward the Minami-dake crater in the radial component and in a direction rotated clockwise 90° from the crater in the tangential component. The increase in the strain represents an extension of the radial component in a direction toward the crater and a direction rotated clockwise 90° and 45° from the crater in the tangential and oblique components, respectively. The two vertical gray lines represent the start of inflation on November 7, 2017 and the onset of the explosion at Minami-dake on November 13, 2017, respectively Hotta and Iguchi Earth, Planets and Space (2021) 73:70 To expand the original elastic source models described above to a multiple-source model, we assumed that the sources do not interact with each other and that the calculated tilt and strain are linear sums of the tilt and strains caused by the sources (Hotta et al. 2016a).

Residual evaluation
We searched for the optimal parameters for each source (Table 1) using a grid search method. For models with more than seven parameters, i.e., composite model of a tensile fault and dual spherical source, we used a genetic algorithm (e.g., Carroll 1999;Tiampo et al. 2000;Murase et al. 2006) instead of a grid search method. In this case, we set the number of individuals to be 50; the probabilities of crossovers, jump mutations and creep mutation to be 70, 2 and 18%, respectively; and the number of iterations for the genetic algorithm to be 200,000 (Hotta et al. 2016a;b). Here, individuals have model parameter information as binary code; a crossover is a random exchange of part of the binary code among an arbitrary pair of individuals; jump mutation is a random change in a bit in the binary code of an individual; and creep mutation is a ± 1 step change in the model parameter information in an individual. To avoid a local minimum, we executed the iteration process several times (Hotta et al. 2016a, b). To minimize the difference between the observed and theoretical values, we used the weighted residual sum of squares (WRSS): where 2 and E 2 are weighted residual sums of squares of tilt and strain, respectively: and where the subscripts obs and cal indicate observed and calculated values; δR i and δT i indicate the radial and tangential tilt change at the i-th tilt station, respectively; εR j , εT j and εO j indicate the radial, tangential and oblique strain changes at the j-th observation tunnel, respectively; and σ indicates the observation error. The search spaces and step sizes for the parameter estimates are given in Table 1.

Phase 1
A single spherical source model was determined for phase 1. Through the grid search method, we obtained a very shallow deflation source beneath Minami-dake (source M S ) at 0.1 km bsl (Fig. 4a, Table 1). Its radius was the minimum value of the search range, 0.1 km. The large strain change detected at AVOT, small tilt change including AVOT and strain change at the other stations can be explained by this source (Fig. 4a and b).

Phase 2 Deformation sources beneath Minami-dake
During phase 2, subsidence in a direction toward Minami-dake was detected except at the KMT observation tunnel where a subsidence in the direction toward the north flank was detected. Excluding the data of KMT, we applied dual spherical sources: a shallow source (source M S ) fixed at the location obtained for phase 1 and a deep source, M D . M D is a deeper (3.3 km bsl beneath Minami-dake) deflation source, whereas M S causes inflation as shown in Fig. 5a and Table 1. We note, however, that these sources beneath Minami-dake cannot explain the oblique strain change at HVOT. These sources also cannot explain the tilt and strain changes at the KMT observation tunnel, where subsidence toward the north area of Sakurajima was detected (Fig. 5). A significant cause of the difference in tilt change between KMT and KOM, which are spatially located close together, is an additional source. Other causes of the difference may be topography or heterogeneity.

Improvement by the addition of sources for KMT data
Although the combination of the shallow and deep Minami-dake sources can explain most of the tilt/strain data, there are still misfits in the oblique strain change at HVOT and tilt/strain data at KMT. In the period of 2011-2012, a source was identified beneath Kita-dake (Hotta et al. 2016a), as mentioned in the introduction, and this source may have re-activated during the explosion of Minami-dake. Thus, we add a third source and test whether it was spherical or dike/sill-like.
Triple spherical source model Hotta et al (2016a) found a Mogi-type deformation source beneath Kitadake in addition to that beneath Minami-dake for the  Table 1 Search range and step for the model parameters and their optimal values and uncertainties (square brackets) for each period The origin of the horizontal coordinate is the summit crater of Minami-dake (130.66° E, 31.58° N). The uncertainties are the 95% confidence intervals estimated from the F-test (Árnadóttir and Segall 1994). The position of the dislocation model (dike) is at the center of the top edge of the fault. The dip is clockwise from the horizontal, and the strike is clockwise from the north. We also checked smaller radius of sources M S and K, but WRSS value did not improve. Therefore, we did not extend the search range of radius any smaller  long-term deformation event of October 2011-March 2012. We applied a triple spherical deformation source model adding a spherical source, K, beneath Kita-dake. The two original sources, M S and M D , are fixed in their previously determined locations (Fig. 5), and we solved for the location of the source K and volume changes of the three sources. The new source is a deflation source beneath the northeast flank of Kita-dake where an inflation source was found by Hotta et al. (2016a) during the long-term deformation period from October 2011 to March 2012 ( Fig. 6a; Table 1). The fit of the calculation to the observed oblique strain change at HVOT and tilt and strain changes at KMT are improved by this addition (Fig. 6), and the WRSS value reduced from 55,494 to 13,408. Hotta et al. (2016b) applied a tensile fault model to the deformation of Sakurajima. Here, we add a tensile fault to M S and M D . In addition to the inflation of M S and deflation of M D , a sill with a length and width of 4.7 and 0.4 km, respectively, and striking in the northsouth direction is located beneath northeast off Sakurajima where submarine eruptions occurred around this area during the An'ei eruption in 1779 at a depth of 0.7 km bsl (Fig. 7a, Table 1). The fit of the calculation to the observed oblique strain change at HVOT and tilt and strain changes at KMT are improved by this addition (Fig. 7), and the WRSS value reduced from 55,494 to 9,140.0.

Selection of the optimal model for phase 2
To determine the optimal solution, for each model we evaluated the correction of Akaike's Information Criterion (c-AIC) values (Sugiura 1978), which is a finite correction of Akaike's Information Criterion (AIC; Akaike The WRSS values of the three spherical sources and the dual spherical sources with tensile fault models are 13,408 and 9,140.0, respectively, and the number of free parameters is seven and ten, respectively. These yield c-AIC values of 222.19 and 233.53, respectively. Although the residual is larger than the dual spherical sources with a tensile fault model, the three spherical sources model is better, based on the obtained c-AIC values. In another respect, the volume increase of the obtained tensile fault is nearly 10 7 m 3 (/1.83 h), which is too large based on the typical magma intrusion rate of Sakurajima, 10 7 m 3 /year (Ishihara 1981). The intrusion rate during the inflation period of October 2011-March 2012 was 1.4 × 10 7 m 3 / year (Hotta et al. 2016a). From GPS observations, a rate as high as 10 7 m 3 in ~ 2 h was not detected around the November 13, 2017 Minami-dake explosion, as shown in Fig. 8, indicating that the intrusion rate at this time was less than on the order of 10 7 m 3 /year. Therefore, we concluded that the three spherical sources model is the optimal model for phase 2.
As for a finite spherical source model, Pascal et al. (2014) noted that the effects of interactions between sources are negligible when the sources are separated by four or more radii. The distance between sources M S and M D is 3.6 km, that between M D and K is 4.8 km, and that between K and M S is 5.7 km. Meanwhile, the optimal radii of M S , M D and K are 0.1 km, 0.5 km and 0.1 km, respectively. From these distances and radii, these three spherical sources are separated enough to ignore the effects of interactions between them.

Interpretations of the sources
During the 3 min from the onset of the explosion at 22:07, phase 1, a very shallow deflation source was obtained beneath Minami-dake, similar to previous eruptions at the Minami-dake (Ishihara 1990) and Showa craters . This means that the deflation starts from a shallow part of the magma plumbing system of Sakurajima. This shallow source might represent the uppermost part of the magma conduit from a deeper Minami-dake source to the summit crater, which has been interpreted as a high-pressure gas chamber (e.g., Iguchi 1994). After the peak of the Vulcanian eruption, phase 2, the deeper source beneath Minamidake deflated while the shallower source inflated. The deep deflation source beneath Minami-dake had been identified based on tilt and strain changes during previous eruptions at the Minami-dake and Showa craters. The obtained depth of 3.3 km bsl was comparable with previous studies (2-4 km), and the horizontal location of the source was similar to the fixed location beneath Minami-dake of previous studies (e.g., Ishihara 1990;Iguchi et al. 2013). In addition to the deeper source, inflation of the shallower source was detected. The shallow inflation source beneath Minami-dake found in the present study was not found in the previous studies. During the Minami-dake eruption including this event, the shallower source is considered to deflate for the first phase of each eruption. In this event, the shallower source inflated after the peak of the explosion probably because the amount of accumulated magma from the deeper source exceeded that of the emission from the Minami-dake crater due to continuous Strombolian eruptions (Fig. 9). In addition to the two Minami-dake sources, we found a deflating spherical source at a depth of 3.2 km bsl at the northeastern flank of Kita-dake. The Kita-dake source might correspond to the source K found beneath Kita-dake at a depth of 3.3 km bsl from a combination analysis of GPS, tilt and strain data during the period from October 2011 to March 2012 (Hotta et al. 2016a). The Kita-dake source was also found based on vertical displacement during  Hotta et al. 2016a) and the assumption that the intrusion rate and the deformation rate are linearly related. Fig. 9 Schematic diagrams of magma migration during the Minami-dake explosion revealed in the present study. The left diagram represents phase 1 (22:07-22:10) and the right represents phase 2 (22:10-24:00). The open and solid circles represent inflation and deflation of the spherical sources, respectively. The thickness of the arrows between each source or the Minami-dake crater represents the amount of magma migration based on the assumption that uncompressible magma migrates from sources K to M D , M D to M S , and M S to Minami-dake crater, and magma is not injected towards source K the period from November 2008 to November 2009, obtained from a precise leveling survey . Not only the deeper Minami-dake source M D but also the Kita-dake source deflated due to the Minamidake explosion (Fig. 9). The total deflation volume of sources M D and K was 88,900 m 3 , which is comparable with the previous Minami-dake eruptions of 10 3 -10 5 m 3 (Ishihara 1990).
A clear reflector at a depth of 5.8 km bsl beneath Kitadake, identified as an overcritical PP reflection, was found by seismic experiments during the period from December 2009 to December 2014. The reflection became clear in December 2014, and the reflector was interpreted as a sill (Tsutsui et al. 2016). The source K obtained in the present study is shallower than the reflector, but the horizontal location is similar. The 99% confidence interval of the depth of source K is 2.6-5.8 km bsl, as shown in Table 1, and the depths of source K and the reflector may overlap considering the uncertainty in the depth of the reflector of ± 90 m (Tsutsui, personal communication on December 17, 2020). If source K and the reflector are the same, the difference in the optimal depth may be due to the source shape, such as a penny shaped source (Fialko et al. 2001), or heterogeneity, which we will discuss for volume changes later. These should be considered closely in future studies.
Recent studies on magma plumbing systems (e.g., Cashman et al. 2017) indicate that magma recharge of upper crustal reservoirs, where silicic melts segregate from mafic magma, grow large magma bodies and become a trigger for eruptions. The mafic magma is intermittently resupplied from deeper levels in the system. As for Sakurajima, mafic magma has been frequently injected since the twentieth century (Nakagawa et al. 2011). The intrusion rate from the Aira caldera to the Kita-dake source is up to 500 m 3 /day during inflation periods (Hotta et al. 2016a) when the SiO 2 content of volcanic glass decreases, which indicates the injection of mafic magma (Matsumoto et al. 2013). This causes inflation of the Kita-dake source, which acts as a buffer for increases in eruptive activity (Hotta et al. 2016a), becoming a trigger of magma migration to deeper and shallower Minami-dake sources. As a result, magma is finally emitted from the Minami-dake or Showa crater, causing eruptions.

Effect of topography and heterogeneity
There are many limitations such as topography and heterogeneity in the homogeneous analytical model we used in the present study (e.g., Del Negro et al. 2009;Gottsmann and Odbert 2014;Hickey et al. 2016;Hotta and Iguchi, 2017). Although realistic modeling using FEM is time consuming, we checked for the effect of volume change for phase 2. We applied the FEM using Flex PDE Professional software version 6.50 provided by PDE Solutions Inc. We set an FEM domain size of 100 × 100 × 50 km 3 (130.13325-131.18675° E, 31.12906-32.03094° N, 0-50 km bsl). The maximum grid size was set to be 30 per 100 km and the grid size in the domain was set automatically by the software. We set the top of the domain to be a free surface, and the sides and bottom of the domain to be fixed boundaries. Topography was introduced in the area within Sakurajima using digital elevation model (DEM) data with 100 m resolution, and the elevation of the outside area was assumed to be sea level. Heterogeneity of the P-wave velocity was introduced referring to Miyamachi et al. (2013), as 1, 2.5, 3.6, 4.8 and 6 km/s for the depth ranges of 0-0.2, 0.2-1, 1-2, 2-3 and 3 + km from the surface, respectively, as shown in Fig. 10. Assuming a Poisson's ratio of 0.25 and V P /V S velocity ratio of 1.73 (e.g., Hotta et al. 2016b), the modulus of rigidity μ can be written using the density of ground ρ as (5) µ = ρV 2 P 3 Fig. 10 Schematic description of P-wave velocity structure introduced in the finite element method in "Effect of topography and heterogeneity" section Currenti et al. (2007) found the following empirical formula between the density of the ground and P-wave velocity: From Eqs. (5) and (6), the modulus of rigidity of the P-wave velocity structure can be calculated to be 1.02, 3.99, 6.84, 10.5 and 15.1 GPa for the depth ranges of 0-0.2, 0.2-1, 1-2, 2-3 and 3 + km from the surface, respectively. Using the FEM domain described above, we first calculated the tilt and strain change due to a pressure change of + 1 MPa inside each source (M S , M D and K) while fixing the horizontal location, depth and radius that were obtained using a homogeneous analytical model. We then searched for the optimal internal pressure change by a grid search method based on the fact that tilt/strain changes are proportionate to internal pressure changes in FEM calculations. The search range and step were set to be − 100-+ 100 MPa and 0.01 MPa, respectively. The optimal internal pressure changes of M S , M D and K were + 0.71 MPa, − 1.41 MPa and − 82.98 MPa, respectively. These yielded volume changes of + 450 m 3 , − 36,200 m 3 and − 16,700 m 3 , respectively. The yielded WRSS value of 15,997 was slightly larger than that of our homogeneous analytical model, but the data can be well explained (Fig. 11). All of the yielded volume change amounts were smaller than those obtained using a homogeneous analytical model, especially that for the shallowest source M S , the result for which (+ 450 m 3 ) was outside of the estimated 99% confidence interval of the homogeneous analytical model (600-5000 m 3 ). This indicates that modeling by a homogeneous analytical model overestimates of the volume change, particularly for the shallow source. The weight of ashfall in November 2017 is estimated to be 1.75 × 10 8 kg. Although six eruptions occurred in November 2017, five eruptions excepting that on November 13, 2017, were extremely small-scale (plume heights were not exceeding 2000 m). Therefore, the weight of ashfall by the November 13, 2017 Minami-dake explosion can be assumed to be nearly 1.75 × 10 8 kg. Assuming a dense rock equivalent (DRE) density of 2500 kg/m 3 , the volume of ejected magma is nearly 70,000 m 3 . Since the total deflation volume of sources M D and K from FEM analysis was 52,900 m 3 and the inflation of the source M S is very small (450 m 3 ), magma was compressed by up to 25% in these deep sources. The pressure change value of the source K is very large. In addition, the 95% confidence interval of the radius of source K is 0.1-10 km. This result indicates that the radius of the source K can be larger than 0.1 km and therefore the pressure change value would be smaller.
In this study, we checked only for volume change. The location, depth and radius of the source should be investigated in a future study.

Conclusions
We modeled the tilt and strain change during the eruption at Minami-dake on November 13, 2017 by an elastic model. During the 3 min from the onset of the explosion at 22:07, a shallow source beneath Minami-dake was deflated, which was interpreted as the uppermost part of a magma conduit from the deeper Minami-dake source to the summit crater. After the peak of the Vulcanian eruption, the deeper source beneath Minami-dake deflated while the shallower source inflated. The deep deflation source beneath Minami-dake had been identified based on tilt and strain changes during previous eruptions at the Minami-dake and Showa craters. Inflation of the shallower source may be caused by accumulated magma from the deeper source exceeding that of the emission from the Minami-dake crater due to continuous Strombolian eruptions. In addition, a deflating spherical source beneath Kita-dake was also found, which  Fig. 6, but for the analysis using the finite element method in "Effect of topography and heterogeneity" section was identified during the long-term deformation event in 2011-2012. We found that the Kita-dake source, which was considered to be in a cycle of inflation and deflation with a period of several years, may also cause deformation during each eruption.