Solidified magma reservoir derived from active source seismic experiments in the Aira caldera, southern Kyushu, Japan

The Aira caldera, located in southern Kyushu, Japan, originally formed 100 ka, and its current shape reflects the more recent 30 ka caldera-forming eruptions (hereafter, called the AT eruptions). This study aimed to delineate the detailed two-dimensional (2D) seismic velocity structure of the Aira caldera down to approximately 15 km, by means of the travel-time tomography analysis of the seismic profile across the caldera acquired in 2017 and 2018. A substantial structural difference in thickness in the subsurface low-velocity areas in the Aira caldera between the eastern and western sides, suggest that the Aira caldera comprises at least two calderas, identified as the AT and Wakamiko calderas. The most interesting feature of the caldera structure is the existence of a substantial high-velocity zone (HVZ) with a velocity of more than 6.8 km/s at depths of about 6–11 km beneath the central area of the AT caldera. Because no high ratio of P- to S-wave velocity zones in the depth range were detected from the previous three-dimensional velocity model beneath the AT caldera region, we infer that the HVZ is not an active magma reservoir but comprises a solidified and cool remnant. In addition, a poorly resolved low-velocity zone around 15 km in depth suggests the existence of a deep active magma reservoir. By superimposing the distribution of the known pressure sources derived from the observed ground inflation and the volcanic earthquake distribution onto the 2D velocity model, the magma transportation path in the crust was imaged. This image suggested that the HVZ plays an important role in magma transportation in the upper crust. Moreover, we estimated that the AT magma reservoir in the 30 ka Aira caldera-forming eruptions has the total volume of 490 km3 DRE and is distributed in a depth range of 4–11 km.


Introduction
Kyushu island is located at the junction of the Southwest Japan arc and Ryukyu arc, where the Philippine Sea plate subducts northwestward beneath the Amurian plate at the rate of 74 mm/year along the Nankai Trough and Ryukyu Trench (Wallace et al. 2009) (Fig. 1a).As mentioned in Zhao et al (2018), a unique tectonic feature of Kyushu is subject to the influence of a northsouth extensional stress in the central part of the island.Under this extensional stress field, there exists a major tectonic structure, known as the Beppu-Shimabara graben, along the northeastern continuation of the backarc basin Okinawa Trough (Matsumoto 1979), and the large and active Aso caldera is also situated in the central southern rim of the graben.In the southern part of the island, there are four large calderas along a nearly straight chain almost parallel to the Nankai Trough: the Kakuto, the Aira, the Ata and the Kikai calderas (Aramaki 1984).Among these calderas, the most active caldera is the Aira distributed in the Kagoshima graben.The rifting activity in the Kagoshima graben is considered to have begun at 3 Ma (Hayasaka 1987).As shown in Fig. 1b, the approximate dimensions of the Aira caldera are 14 × 14 km extended across ~ 200 km 2 , and its present topographic outline is characterized by steep caldera walls (Geshi et al. 2020).Furthermore, the small Wakamiko caldera, with approximate dimensions of 8 × 6 km, is situated in the northeastern part of the Aira caldera (Moriwaki et al. 2017;Geshi et al. 2020).
An area in and around the Aira caldera is widely covered by Quaternary pyroclastic flows and volcanic rocks (Fig. 1b).Pleistocene sediments, tuff, and tuff breccia can be found under these surface pyroclastic flows (Nagaoka et al. 2001;Uchimura et al. 2014).Along the eastern inland area of the caldera, Upper Cretaceous Shimanto Group strata, consisting of sandstone and mudstone alterations, are widely distributed and regarded as the geological basement in southern Kyushu (Hayasaka 1987;Uchimura et al. 2014).Moreover, the Shimanto Group is recognized to be distributed at depths of > 1 km across this region (Aramaki 1984;Hayasaka 1987).Comparatively, the Kekura Formation, a marine sedimentary layer consisting of sand, silt, and tuff, lies above the Shimanto Group (Uchimura et al. 2014).
The Aira caldera-forming eruptions occurred approximately 30 ka (Okuno 2019;Geshi et al. 2020).However, the original caldera is considered to have existed before 30 ka (Aramaki 1984).In the eruptive history of the Aira caldera from 100 to 30 ka (Nagaoka et al. 2001), many small or medium explosive and effusive eruptions had occurred.These eruptions were concentrated in the eastern, north-eastern and southern areas of the Aira caldera, and the magma of these eruptions might be supplied from the center area in the Aira and Wakamiko calderas.Nagaoka et al. (2001) suggested that the ascending magma and formation of the initial magma reservoir began at least 100 ka.The Aira caldera-forming eruptions at ~ 30 ka consisted of five events (Tsukui and Aramaki 1990).According to Geshi et al. (2020), these are as follows: the Osumi pumice-fall deposit (with a volume of ~ 27-31 km 3 dense rock equivalent (DRE)), Tarumizu pyroclastic flow deposit (~ 10 km 3 DRE), Tsumaya pyroclastic flow deposit (~ 6 km 3 DRE), Kamewarizaka breccia, and the Ito ignimbrite and its coignimbrite ash ("Aira-Tanzawa ash").Aramaki (1984) reported that the caldera collapsed at the beginning of the Ito ignimbrite eruption.The last two pyroclastic eruptions (Ito ignimbrite and its coignimbrite ash) were likely to produce Fig. 1 a Tectonic map of Kyushu Island.Blue dashed lines are the Beppu-Shimabara graben "BSG" (Matsumoto 1979) and the Kagoshima graben "KG" (Hayasaka 1987).Red solid lines show five major calderas, namely, Aso, Kakuto, Aira, Ata and Kikai, in Kyushu.A red arrow shows Philippine Sea plate-Amurian plate relative motion (72 mm/years) (Wallace et al. 2009).b Enlarged geological map in and around the Aira caldera.An outline of the geological distribution of major rocks in and around the calderas is described (modified from Uto et al. 1997 andUchimura et al. 2014).Black broken line indicates the projection line of the present seismic profile pyroclastic deposits of ~ 350 km 3 DRE and collectively referred to as the "AT eruption".The post-caldera activity began approximately 26 ka and resulted in the formation of the Sakurajima volcano located near the southern rim of the Aira caldera, which remains the most active parasitic volcano to date.
The magma reservoir associated with the Aira caldera-forming eruptions has been investigated by many geological and petrological studies (Tsukui and Aramaki 1990;Yasuda et al. 2015;Geshi et al. 2020).In these studies, the petrographical characteristics, whole-rock chemical composition and water concentrations in glass inclusions in pumice-fall deposits and pyroclastic flow deposits were examined.Based on these features, pressure conditions and the depth range of the magma storage were estimated.Tsukui and Aramaki (1990) reported the top of the magma reservoir was distributed at depths of 8-10 km, whereas Yasuda et al. (2015) revised a thickness of the magma reservoir to 2-2.5 km, with a top depth of 4-5 km.Geshi et al. (2021) also presented a new magma reservoir model with a thickness of 4.4 km and a top depth of 5.3 km.In geodetic studies, Hotta et al. (2016) estimated a distribution of spherical pressure sources by applying a multiple-source model to the geodetic data and interpreted the present magma reservoir located at a depth of 9.6 km in the center of the Aira caldera.Employing finite element analyses to geodetic GNSS data, Hickey et al. (2016) proposed that an oblate pressure source was located at a depth of 13 km in the northeastern Aira caldera.Seismological studies have been conducted to estimate the velocity structure beneath the Aira caldera.For example, Alanis et al. (2012) determined the regional velocity structure via 3D travel-time tomography using local earthquakes in southern Kyushu.They detected a remarkably low S-wave velocity zone at a depth of 20 km in the center of the Aira caldera.From their model, however, it was difficult to investigate the detailed velocity model due to the low spatial resolution resulting from the coarse grid configuration.More recently, Tameguri et al. (2022) presented a local 3D velocity structure in and around the Aira caldera using travel-time tomography with a fine grid configuration, revealing a velocity zone > 20% lower than the initial S-wave velocity at a 15 km depth.Iguchi et al. (2009) conducted an active source seismic experiment focused on the Sakurajima volcano and Miyamachi et al. (2013) analyzed the observed travel-time data and imaged the 2D shallow velocity model to depths of ~ 4 km in the Sakurajima volcano and Aira caldera.However, the structure in the deeper part has been remained unclarified.
It is most important to understand the present structure in the Aira caldera; thus as a first step, this study aimed to clarify a 2D detailed P-wave velocity structure down to 20 km depths, focusing particularly on the structure beneath the Aira caldera.To address this aim, extensive active source seismic experiments were undertaken in 2017 and 2018 along the 195 km E-W transect crossing the Aira caldera region.Furthermore, we attempted to reveal the relation between the derived P-wave velocity model and the volcanic caldera system in the Aira caldera.

Travel-time data
The active source seismic experiments along the 195 km E-W transect crossing southern Kyushu was completed in 2017 and 2018, as shown in Fig. 2. The entire profile comprised 6 sub-profiles: Profile-A, -B, -C, -D, -E, and -F.
Onshore (Profile-A, -C, and -E), 829 vertical component seismometers (4.5 or 5 Hz) were deployed at 100-m intervals.Across the offshore areas, 18 ocean bottom seismographs with 4.5 Hz three-component seismometers were deployed at 1 km intervals in Profile-B; whereas 24 ocean bottom seismographs were deployed at 2 km intervals in Profile-D.There were no seismic stations along Profile-F.At all stations, seismic waves were digitally recorded using a 24 bit analog-to-digital conversion, with a 200 or 250 Hz sampling frequency.
In 2017, nine explosive sources (SP1-SP9) with a 200 kg dynamite or water gel charge were detonated.In 2018, a total of 206 airgun sources using an airgun array of 6000 in 3 (2000 psi), towing a 600 m-long streamer cable (48 channels) were fired at 50 m intervals in the western half of Profile-B.Along Profile-D and -F, 495 and 272 airgun sources towing a 2100 m-long streamer cable (168 channels) were fired at 100 m intervals, respectively.In addition, two vibroseis sources, consisting of four vibrator-vehicles, were conducted at two fixed points, one at the western end of Profile-A (referred to as the VA point) and another at the eastern end of Profile-C (referred to as the VC point).
Figures 3, 4 show examples of the observed seismic waves derived from the onshore sources (SP1-SP9) and offshore airgun sources in Profile-B and Profile-D, respectively.The first P-wave arrivals with an onset uncertainty < ± 0.10 s were picked from all observed seismic waves.A total of 7224 travel-time data from the onshore sources SP1-SP9, VA and VC and 118,306 data from the offshore airgun sources were obtained.These picked travel-time datasets are used to derive the detailed structure beneath the surveyed area.

Method of travel-time tomography
A 2D P-wave fine velocity distribution along the whole profile was determined using the travel-time tomography demonstrated in Shiraishi et al. (2010).The entire study area was divided into fine square cells, and an inversion algorithm was created based on the simultaneous iterative reconstruction technique (Trampert and Leveque 1990).Theoretical travel-times were calculated using the linear travel-time interpolation (Asakawa and Kawanaka 1993), and the velocity value at any point in the model was evaluated by averaging the velocity values of the surrounding cells.Additionally, in accordance with Shiraishi et al. (2010), tomographic inversions for 100 sets of randomly generated initial models were executed to reduce the dependence of the initial model on the tomographic solution.The mean velocity value (i.e., the final solution) and its standard deviation were then calculated for each cell from the 100 sets of tomographic solutions.The small or large standard deviation in a cell implies that the cell is well or badly resolved, respectively.
In this study, the locations of all stations, shots, and travel-time data were projected on a straight line that connected the easternmost station in Profile-A and the western end of the airgun sources in Profile-F (Fig. 2).First, the entire profile region (ranging in distance from 0 to 195 km and in depth from − 1 to 60 km) was divided into a 200 m square section (hereafter, called the 200 m cell model), and tomography was applied to all traveltime data.Second, to obtain a fine structural image for the Aira caldera, tomographic analyses were conducted only for the central profile (Fig. 2) at distance between 20 and 118 km (Profiles-A, -B, and -C) with a finer 100-m square (hereafter, called the 100-m cell model), where denser ray coverage was expected.The target region and the travel time data used in the 100-m cell model are different from those in the 200-m cell model.

Results of the tomography
Figure 5 plots the root mean square (RMS) of travel-time residuals in the tomography for the initial 100 models adopted in the 200 m and 100 m cell analyses.The mean values of the RMS residuals for the 200 m and 100 m cell models decreased from 5.36 to 0.14 s and 2.39 to 0.10 s in the 9th iteration, respectively.Although the RMS residuals showed a large variation at the early stage of the iteration process depending on the initial models used, the tomographic solutions predominantly converged within 8 iterations.The final RMS residuals in both cell models were also comparable to the P-wave onset uncertainty (± 0.10 s) from the travel-time picking.The convergence in the RMS residual does not imply that all 100 solutions obtained by the 100-m tomographic model converge to a similar solution.
Figures 6, 7 show the final velocity structures, their standard deviations, and the ray densities in the 200 and 100 m cell models, respectively.Solution reliability remained high for both cell models to depths of ~ 10 km because of the high ray density distribution and small standard deviations (< 0.1 km/s).Over a depth range of 10-15 km in the central profile, only the solutions in the 200 m cell model appeared stable owing to the medium standard deviations (< 0.2 km/s).In regions deeper than 15 km, however, the solutions in both cell models were less reliable because of poor ray density and larger standard deviations (> 0.3 km/s).Figure 8 shows an example of histograms of traveltime residuals (O-C = the observed travel time-the calculated one) for onshore sources (SP3-SP8) in the 100-m cell models.In this study, we assumed that an allowable misfit was twice the P-wave onset uncertainty, namely ± 0.20 s.Notably, 73% and 64% of all residual data in the 200 m and 100 m cell models, respectively, ranged within the allowable misfit.These results suggest that the velocity model in each model obtained by tomography roughly fulfills the travel-time data.Moreover, 5% of all data in both cell models have large misfit values over ± 0.50 s and the systematically large misfit is found especially in the onshore section on Profile-C for SP5 (Fig. 3e).There is no definite explanation for the large misfits; however, there are several possible causes: (a) insufficient ray coverage (in the sedimentary layers) over the offshore areas, (b) local velocity variations in the complicated shallow region, and (c) smoothing processes for tomographic modelling.

Checkerboard and restoring resolution tests
The checkerboard resolution test (CRT) and restoring resolution test (RRT), based on Zhao et al. (1992) and Nakajima et al. (2010), respectively, were adopted to evaluate the reliability of the spatial resolution of the obtained velocity model.
In the CRT, an initial one-dimensional (1D) velocity cell model was assumed.The true model was defined; positive and negative velocity perturbations of 5% were assigned to alternating cells in the 100 and 200 m cell models at 10 km horizontal and vertical intervals.Synthetic travel-times were calculated for all stations and shot pairs used in the actual tomographic analyses.The calculated synthetic travel-time data were inverted using the initial 1D model without any velocity anomalies.All CRT results are shown in Fig. 9.In a horizontal distance range of 50-80 km in the both cell models, the upper region (depth range < ~ 10 km) appeared to be resolved.These resolved regions in Fig. 9 agree well with the ray high density and small standard deviation regions shown in Figs. 6, 7. Comparatively, there was limited resolution in the deeper regions of both models.As noted by Zhao et al. (1992), the simple CRT may be unsuitable for evaluating reliability of spatial In the RRT, synthetic travel-times were calculated using the synthetic model (upper panels in Fig. 10a,  b) obtained by inverting the observed data.Synthetic travel-times were then inverted to obtain a restored model by tomography (lower panels in Fig. 10a, b).A comparison of the synthetic and restored models demonstrated that the main velocity trends and characteristic structures in both cell models were well-restored.
From the results of Figs. 6, 7, 9 and 10, we concluded that the tomographic solutions to depths of ~ 10 km in the central profile were reliable in the 100-m cell model, whereas those of the 200 m cell model were deemed reliable to depths of ~ 15 km.

Subsurface region characteristics (0-5 km depth)
Figure 11 shows the marine reflection section image along the western half of Profile-B created by applying common midpoint reflection stacking, post-stack time migration and depth conversion (Yilmaz 2001).For the depth conversion, the result of the tomographic analysis was used as the velocity model.Because of the relatively short length of the 500 m streamer cable (40 channels), the reliable depth range was limited to 2-3 km.
In Fig. 11, a horizontally stratified sedimentary structure (Holocene) was observed at a depths of 0.7 km.Marked reflection events occurred at a depth range of 0.5-0.7 km, depicting the top of the Pleistocene Kekura Layer (Hayasaka 1987).The strong reflection events at a depth range of 1.0-1.4km represent a geological boundary between the Kekura Layer and the Shimanto Group basement, distributed widely around the Aira caldera.Some small vertical faults < 1 km in length (yellow broken lines in Fig. 11) were also recognized at depth ranges of 1-2 km.
The velocity and its perturbation distributions from the 100-m cell model are shown in Fig. 12.The velocity perturbation (%) was defined as the deviation of the velocity of a given cell from the mean velocity at each depth and averaged from 20 to 118 km (Profile-A, -B, and -C).This figure (Fig. 12b) clearly indicates that the shallow structure inside the Aira caldera area is dominated by the wide and thick low-velocity region.Most notably, the thick low-velocity material abruptly appears east of SP6.

Upper region characteristics (5-11 km depth)
The upper crustal structure beneath the entirety of the Aira caldera is shown for the 100 m cell model in Fig. 12, displaying the complicated structural variation observed across a depth range of 5-11 km.
Firstly, an extremely high-velocity zone (HVZ) is distributed over a depth range of 6.5-11 km beneath the center of the Aira caldera.This HVZ has a P-wave velocity > 6.8 km/s (i.e., 6% faster velocity perturbation), with a horizontal extent of 9.5 km and a thickness ranging from 1.0 to 2.5 km, inclining towards the west at ~ 20°.Higher velocity is concentrated in the eastern part of the HVZ.The HVZ is also identified in the 200 m cell model (Fig. 13), indicating that the imaged HVZ is independent of the cell size in the tomography analysis.Additionally, the main part of the HVZ is located just beneath the AT caldera.
Second, a low-velocity zone (LVZ), ~ 5 km in width, was recognized at a depth of ~ 8 km depths beneath the Wakamiko caldera (hereafter referred to as the upper LVZ); however, this LVZ could not be distinguished in the 200 m cell model (Fig. 13), indicating that the magnitude and extent of the upper LVZ are not well resolved from our travel-time data coverage.

Lower region characteristics (> 11 km depths)
To examine the deeper velocity structure, the velocity section from the 200 m cell model is shown in Fig. 13.The velocity contour lines of 6.0 and 6.2 km/s were convex beneath the center of the Aira caldera, while the HVZ mentioned in previous section was located at the top of this convex structure (Fig. 13a).Comparatively, an equivelocity line of 6.4 km/s abruptly deepened from 12 to 17 km towards the east and showed a valley-like geometry at the bottom, suggesting the existence of another weak low-velocity zone around a depth of ~ 15 km (Fig. 13b), which is hereafter referred to as the deep LVZ.

Comparison with previous velocity models
Using a subsurface sonic survey in the sea area of the Aira caldera, Hayasaka (1987) reported that the western half area of the Aira caldera, which corresponds to the current seismic reflection line (Fig. 11), did not originate from the caldera depression structure, because the sedimentary layers with obvious bedding surface are notably recognized in the area.Hayasaka (1987) also considered that the eastern half of the caldera originated from the volcanic activity and the caldera depression structure owing to the existence of a fracture zone with prominent normal faults a few kilometers wide between the Aira and Wakamiko calderas.Unfortunately, the reflection profile in this study does not cross this fracture zone.This fracture zone is almost identical to the volcanically active zone of the small and medium eruptions that occurred during the 80-30 k.y.period estimated by Nagaoka et al. (2001).
From the 2008 seismic experiment (Iguchi et al. 2009), Miyamachi et al. (2013) estimated the velocity structure along a NW-SE profile across the sea area of the Aira caldera and revealed a surface sedimentary layer with a P-wave velocity of 2.3-2.8km/s, ~ 1 km in thickness, as well as a low velocity zone (LVZ) with a P-wave velocity of 4.2-4.4km/s over a depth range of 1.4-3.2km in the center of the Aira caldera.This surface sedimentary layer and LVZ correspond well with the horizontally stratified sedimentary layer shown in Fig. 11, in addition to the shallow, thick, low-velocity region in the AT caldera indicated in Fig. 12b.Assuming the P-wave velocity at the bottom of the low-velocity region is 4.4 km/s estimated by Miyamachi et al. (2013), the extent of the  2013), respectively.Volcano-tectonic (VT) earthquakes (white circles) are also plotted.Other symbols are same as those in Fig. 7 low-velocity region is evaluated in this study.As a result, in the western part of the Aira caldera, the 4.4 km/s contour line extends to a depth of > 4.3 km.Comparatively, this contour line gradually becomes shallower to the east and reaches 2 km in depth at the eastern end of the caldera.The western boundary of the Aira caldera rim is sharply defined just west of SP6, while its eastern boundary remains less well defined.The velocity perturbation in Fig. 12b clearly depicts this remarkable structural difference: the bottom of the low-velocity region in the western part is located at a depth of ~ 4.3 km, whereas in the eastern part, it occurs at a depth of ~ 1.7 km.This noticeable difference in thickness implies that the lowvelocity region in the eastern part is distinct from that in the west.The western low-velocity region (with a 12 km horizontal extension) corresponds to the AT caldera associated with the AT eruption 30 ka, while the eastern low-velocity region may be an expression of the younger Wakamiko caldera (13 ka) with a horizontal dimension of 6 km (Geshi et al. 2020).These results suggest that the  2013), respectively.Volcano-tectonic (VT) (white circles) and low-frequency (LF) (yellow circles) earthquakes are also plotted.Other symbols are same as those in Fig. 6 greater "Aira caldera" comprises at least two sub calderas: the AT and Wakamiko calderas.
Alanis et al. ( 2012) investigated a 3D regional velocity structure in southern Kyushu via tomography using local earthquake travel-time data and revealed a region with a high ratio of Vp/Vs region at a depth of 20 km beneath the Aira caldera.More recently, Tameguri et al. (2022) revised the 3D local velocity structure in and around the Aira caldera via tomography by adopting a finer 5-km grid configuration in both the horizontal and vertical directions.They found that the high Vp/Vs region was distributed at a depth of 15 km beneath the center of the Aira caldera and was interpreted as the active magma reservoir, which likely corresponds to the deep LVZ in Fig. 13 (although the horizontal location of their high Vp/ Vs region is ~ 5 km west of the current deep LVZ).Alternatively, the upper LVZ over depth ranges < 5 km (Fig. 12) and the HVZ ranging in depth between 6 and 11 km (Fig. 13) in the current model were not observed in the 3D velocity model of Tameguri et al. (2022).Kiser et al. (2016) clarified the existence of a relatively high P-velocity zone with a high Vp/Vs corresponding to the active magma reservoir in the depth range of 4-12 km beneath Mount St Helens.In the case of the Aira caldera, the 3D velocity model by Tameguri et al. (2022) showed a high Vp/Vs zone with an evident low S-wave velocity at ~ 15 km depth.These two studies suggest that active and probably partially melted magma reservoirs are characterized as high Vp/Vs zones.In and around our mapped HVZ, in contrast, the Vp/Vs values imaged by Tameguri et al. (2022) were not obviously high.The characteristics of the HVZ suggest a volcanic origin but is composed of the solidified magma in contrast with the deep LVZ at a depth of 15 km, although the hypothesis that the HVZ was emplaced much earlier than the AT eruption is not discarded.

Seismic activity in Aira caldera
Figure 14 shows the Japanese-unified hypocenter distribution in and around the Aira caldera provided by the Japan Meteorological Agency from 2000 to 2019.These hypocenters are classified into two categories: volcanotectonic (VT) and low-frequency (LF) earthquakes.
From the VT epicenters (Fig. 14a), high seismic activity was recognized in the Sakurajima inland, as well as in the fracture zone (Hayasaka 1987) along the western edge of the Wakamiko caldera.There were no VT earthquakes inside the HVZ region and the Wakamiko caldera.Recently, however, Miyamachi et al. (2022) conducted a hypocenter determination using arrival-time data obtained by ocean bottom seismograph observations in and around the Wakamiko caldera.They noted that many of the Japanese-unified VT earthquakes located in the fracture zone had relocated to inside the Wakamiko caldera.
Most VT earthquakes occurred at depth ranges of < 11 km, and the LF earthquakes were primarily distributed over depths of 17-30 km.From the distribution of the VF and LF earthquakes, an apparent aseismic region was distributed over a depth range of 12-18 km (Fig. 14b), corresponding to the deep LVZ observed in this study (Fig. 13).
Figures 12 (100 m cell model) and 13 (200 m cell model) show that the VT and LF hypocenters within widths of 6 km along the profile were projected onto the velocity sections.In Fig. 12, the VT hypocenters over a depth range of 5-10 km were vertically concentrated in a region between the HVZ and upper LVZ.Contrastingly, VT earthquakes rarely occurred in or around the HVZ.This fact is important in understanding the role of the HVZ in the magma transportation path.The vertical distribution of LF earthquakes deeper than 18 km may be interpreted as magma ascent in the lower crust (Fig. 13).

Pressure sources and velocity model
Pressure sources beneath the Aira caldera deduced from GNSS data have been previously estimated by a variety of geodetic studies, shown in Figs. 12, 13.For example, Iguchi et al. (2013) analyzed the horizontal displacements measured using GNSS from September 2009 to October 2010 by assuming two sources based on the alleged "Mogi model" (Mogi 1958).They estimated a major source ("pressure source-B" in both figures) located at a depth of 12 km beneath the center of Aira caldera and an additional source at a depth of 5 km in the Sakurajima volcano.By applying a three-source model to continuous-combined geodetic data from a GNSS, as well as tilt, and strain meters over a period from October 2011 to March 2012, Hotta et al. (2016) presented the spherical pressure source-A, located at a 9.6 km depth beneath the center of the Aira caldera.They also revealed the spherical pressure sources-K and -M in the Sakurajima volcano (Fig. 14).Non-spherical pressure sources were proposed by Hickey et al. (2016), who identified an oblate pressure source at a depth of 13 km in the northeastern Aira caldera by applying finite element analysis to geodetic GNSS data from 1996 to 2007.A portion of the deep LVZ in the current model was included within their oblate pressure source region.We interpret these pressure sources represent magmatic activity corresponding to the individual phases of observed crustal deformation.The pressure source-A, the latest activity with respect to our seismic experiment, is located just beneath the HVZ in the current velocity model (Fig. 12).This spatial relationship led to the interpretation that magma ascending from the deeper region was impeded owing to the roof-like HVZ shape.As a result, magma may be accumulating beneath the HVZ behaving as a pressure source.
In general, a pressure source originates from a change in the magma supply volume.From a volume of pyroclastic ejecta in the Sakurajima volcano, a mean magma supply rate is estimated as ~ 1 × 10 -2 km 3 years −1 (Ishihara 1981).If this supply rate continued for 1000 years, the total volume of magma accumulated in the pressure source would reach ~ 10 km 3 , corresponding to a spherical volume with a diameter of ~ 2.7 km.However, the supplied magma is probably released by the eruptive activity of the Sakurajima volcano.Although the active magma reservoir is estimated as a pressure source from the geodetic data (e.g.Hotta et al. 2016), its qualitative size remains unclear.We found that the dimensions of the pressure source area may be laterally small (< ~ 2 km) and are not resolved by our active source seismic data.Iguchi et al. (2013) and Hotta et al (2016) presented magma transportation models beneath the Aira caldera and Sakurajima volcano based on their pressure sources determined from geodetic data.Combining their interpretations with results of the present study, we show a schematic model of magma transportation beneath the Aira caldera in Fig. 15.We speculate that magma in a lower crust ascends through the focal region of LF earthquakes and accumulates in the deep LVZ at ~ 15 km depth.Afterwards, the magma moves toward the HVZ to reach the pressure source-A (Hotta et al. 2016) and is likely a smaller temporarily magma reservoir.The low seismicity surrounding the HVZ (see Figs. 12,14) suggests that the ascending magma does not penetrate the HVZ, but likely moves upward along the lower boundary of the inclined HVZ.

Magma transportation
According to Hotta et al. (2016), the magma in pressure source-A moves to the shallow pressure source-K located at a depth of 3.3 km in the northern part of the Sakurajima volcano (Fig. 14).However, its transportation path from the pressure source-A to the pressure source-K remains unclear.The present study, on the other hand, revealed the detailed 2D velocity model in the Aira caldera along the E-W direction but not in the N-S direction.However, a gravity anomaly in the Aira caldera is distributed concentrically around the center of the caldera (Yokoyama and Ohkawa 1986;Yamamoto and Shichi 2004).Assuming that the crustal structure has a roughly concentric geometry about the center and the HVZ also extends in the N-S direction (as shown in Fig. 14a), the pressure source-K is found to be situated near its southern end.From this spatial relationship, the magma is considered to move along the lower boundary of the HVZ and ascend vertically beneath the pressure source-K.Consequently, the inclined HVZ was considered to play an important role in the magma transportation path.

Solidified and deep magma reservoirs
Our velocity model shows a deep LVZ at ~ 15 km depth (the area enclosed by a black solid curve in Fig. 13), although its spatial resolution is not high.This deep LVZ is nearly in the same location as the aseismic region in Fig. 14.Furthermore, this location is near the high Vp/ Vs region detected by Tameguri et al. (2022) and in the vicinity of the oblate pressure source revealed by Hickey et al. (2016).These correspondences suggest that the deep LVZ represent the active magma reservoir, hereafter referred to as the deep magma reservoir.
Next, we discuss the inclined HVZ, which is defined as the 2 km thick region with a P-wave velocity > 6.8 km/s and is localized in the western part of the Aira caldera (Fig. 12).
First, we consider the origin of the HVZ from geological point of view.The Shimanto Group (Fig. 1b), which is the Cretaceous-Paleogene accretionary complex composed of alternating sandstone and mudstone layers, is widely distributed as a basement rock in and around the Aira caldera (Hayasaka 1987).According to seismic refraction experiments (Ando et al. 2002;Miyamachi et al. 2013), the Shimanto Group is situated in a depth range of 1-4 km with a P-wave velocity of 4.6-4.9km/s.The igneous intrusive body, which penetrated the Shimanto Group at 13-14 Ma (Uchimura et al. 2014), is located at shallow depths of 1-3 km approximately 6 km southeast of the Aira caldera rim, with a velocity of 5.2-5.4km/s (Miyamachi et al. 2013).Considering a higher velocity of the HVZ, however, it is difficult to iconclude that the HVZ originated from the Shimanto Group or as an igneous intrusive body.
Meanwhile, the early volcanic activity in and around the Aira caldera started at 3.1-2.9Ma (Kaneoka et al. 1984), and the activity has continued up to the present.The caldera-forming AT eruptions originated from a rhyolitic magma at 30 ka (Tsukui and Aramaki 1990) and the dacitic-andesitic magmas eruptions in the Sakurajima volcano after the AT eruption are well mapped.However, the high P-wave velocities (> 6.8 km/s) suggest a mafic composition of the magma.Thus, it is difficult to reconcile the interpretation of HVZ being rhyolitic or dacitic-andesitic magmas.During the caldera-forming process, basaltic magma eruptions were identified at 550 ka and 350 ka along the caldera rim (Sudo et al. 2001), at 100-95 ka near the northwestern area of the caldera rim (Nagaoka et al. 2001), and at 8.2-8.1 ka near Fig. 15 Schematic model of the magma transportation path in the Aira caldera.Gray and orange ellipses depict the solidified and deep magma reservoirs, respectively.Red and pink arrows express the paths of magma and heat transportation, respectively.Two orange cross symbols indicate the pressure sources estimated in previous studies (Iguchi et al. 2013;Hotta et al. 2016).Volcano-tectonic (VT) (white circles) and low-frequency (LF) (yellow circles) earthquakes are also plotted.Other symbols are same as those in Fig. 6 the northwestern area of the caldera rim (Moriwaki et al. 1986;Okuno 2019).Based on petrological evidence, Takahashi et al. (2013) reported mantle-derived rocks (basaltic) from 61 ka eruption.From these facts, the HVZ is presumed to represent a consolidated part of basaltic magma in the Aira caldera.Moreover, pooled basaltic melt will differentiate into its silicic and mafic components, and the HVZ alternatively may represent a mafic remnant or cumulate material from fractional crystallization prior to the AT eruption.Regardless, the origin of the HVZ in relation to the Aira caldera eruption history requires further investigation.
Second, we focused on the eruption history of the Aira caldera.According to Nagaoka et al. (2001), the formation of the original magma reservoir of the Aira caldera started in its central area at least 100 ka, followed by small and medium eruptions from 100 to 30 ka and the final caldera-forming eruptions at 30 ka.The magma supply to the original reservoir had been possibly continued during the period of 100-30 ka.The cumulative discharge mass of magma in this period was estimated to be less than 63 × 10 12 kg (Nagaoka et al. 2001), whereas that of the Aira caldera-forming eruptions at 30 ka amounts to more than 414 × 10 12 kg (Aramaki 1984).These esti- mates suggest that the magma supply during a period of 100-30 ka made only a small contribution to the growth of the original magma reservoir.
Third, we consider the possibility of the magma cooling.This cooling process is complicated and affected by various factors, such as the physical properties of magma, heat conduction, convection, and crystallization (Matsumoto et al. 2006).Ehara et al. (2001) performed numerical simulations based on a simple 2D hot-water convective model with a conductive cooling heat source (i.e., intrusive trapezoidal magma), where an initial temperature of 1000 °C was assumed over a depth range of 3-5 km.They showed that a growth of ~ 1 km-thick solidified outer shell of the magma requires 32 k.Hence, it is reasonable to think that the magma possibly has been cooled and solidified to be HVZ since the AT eruption at around 30 ka.
The total magma volume of the Aira caldera-forming eruptions at 30 ka, which consisted of five events, was estimated as 400 km 3 DRE (Ueno 2007).The AT eruption at 30 ka had the volume of 350 km 3 DRE while the total volume of other events was small (Nagaoka et al. 2001;Ueno 2007;Geshi et al. 2020).In the pre-AT eruptions (100-30 ka), the total magma volume was estimated to be less than 22 km 3 DRE (Geshi et al. 2020).In the post-AT eruption (30 ka to the present), 17 Plinian and sub-Plinian eruptions occurred around the AT caldera rim (Kobayashi and Tameike 2002; Okuno 2019).However, the largest P14 eruption at13 ka ejected ~ 11 km 3 tephra, and the total volume of ejected tephra from other eruptions was less than ~ 1.3 km 3 (Kobayashi et al. 2013).Hence, the eruptions in the pre-and post-AT eruption could not make significant contribution to the magma supply in the volcanic activity in the Aira caldera.Therefore, we hypothesize that the magma remaining in the AT magma reservoir has been continuously cooled since 30 ka, particularly in the dormant period between 24 and 13 ka, and was detected as the HVZ in this study.

Magma reservoir associated with Aira caldera formation
Caldera-forming eruptions, the AT eruption at ~ 30 ka, including the associated magma reservoir, have been studied using geological (Aramaki 1984;Nagaoka et al. 2001) and petrological (Tsukui and Aramaki 1990;Yasuda et al. 2015;Geshi et al. 2021) approaches.In contrast, there were no geophysical estimates of this magma reservoir.In this section, a volume of the AT magma reservoir at 30 ka is estimated from the present velocity model.For this estimate, the following four points were presupposed: (1) the AT magma reservoir is cylindrical; (2) the horizontal extent of the magma reservoir is not well known, and we assume it comparable to the HVZ (9.5 km); (3) at the beginning of the AT eruption, most of the associated magma accumulated in the AT magma reservoir, while supply to the reservoir was limited; and (4) the maximum volume of the AT magma reservoir is the sum of the magma ejected during the AT eruption and the remaining in the AT magma reservoir, the latter of which is considered to contribute to the solidified magma reservoir.
The concentric distribution of the gravity anomaly in the Aira caldera (Yokoyama and Ohkawa 1986;Yamamoto and Shichi 2004) suggests that the solidified magma reservoir in the E-W section generally extended to the N-S.Taking the diameter of a cylindrical magma reservoir to be the horizontal dimension of the HVZ (~ 9.5 km), the area of the cylindrical magma reservoir is ~ 70 km 2 at its bottom.As the mean thickness of the HVZ, namely the solidified magma reservoir, is 2 km, the resultant volume of the cylindrical reservoir of the solidified magma reservoir is calculated as 140 km 3 (Fig. 16).Meanwhile, the pyroclastic ejecta from the AT eruption is estimated as ~ 350 km 3 DRE (Ueno 2007).This corresponds to the 5 km-thick cylindrical magma reservoir.Therefore, the total thickness of the cylindrical magma reservoir becomes 7 km and the total volume of the AT magma reservoir is estimated as 490 km 3 DRE.This means that 70% of the total volume was ejected during the AT eruption with the remaining 30% left in the reservoir.Because the lower depth of the solidified magma reservoir is ~ 11 km, the cylindrical magma reservoir ranges from 4 to 11 km in depth.As shown in Fig. 16, this upper depth of 4 km is consistent with the bottom of the subsurface low-velocity area in the present velocity model.
Based on petrological analysis, the chemical composition of the phenocrysts and their melt inclusions in pyroclastic deposits of the Aira caldera-forming eruption of 30 ka, Yasuda et al. (2015) estimated a upper depth of the magma reservoir to be 4-5 km with a thickness of 2-2.5 km.Geshi et al. (2021) also analyzed water content in quartz glass embayments and inclusions from pyroclastic deposits throughout the AT eruption sequence.Through estimation for change of decompression inside the magma reservoir, they proposed an evolution of decompression of magma reservoir and deformation of the host rock toward the caldera collapse and reported that the magma reservoir at the initial stage in the evolution was distributed in a depth range of 5.3-9.4km.Because our present magma reservoir model shows an upper depth of 4 km and thickness of 7 km (= 5 km for pyroclastic deposit (ejecta)) + 2 km for the solid magma reservoir), our estimate is consistent with the petrological model by Yasuda et al. (2015) or Geshi et al. (2021).
Lastly, as discussed above, the Aira caldera-forming eruption requires a significantly large magma reservoir.The present velocity model, however, shows that there is no such LVZ or HVZ along the seismic profile down to a depth of 10 km.This suggests that there may be no active large upper crustal magma reservoir in the present Aira caldera region.

Conclusions
In this study we developed a detailed 2D P-wave velocity model across a 195 km-E-W transect crossing southern Kyushu, Japan, and beneath the Aira caldera.
Variation in the thickness of the low velocity layers in the subsurface depth of 0-5 km suggested that the composite Aira caldera comprises at least two sub calderas, the AT and Wakamiko calderas.The boundary region between two calderas is a fracture zone, which seems to be volcanically active with small and medium eruptions during a 80-30 k.y.period.
The high-velocity zone (HVZ), with a P-wave velocity of more than 6.8 km/s at depths of 6-11 km, was observed at the center of the Aira caldera.This HVZ may constitute a part of the main magma reservoir in the Aira caldera, having existed for ~ 100 k.y.Considering the magma volume ejected in the pre-AT, AT and post-AT eruptions, the AT caldera eruption at 30 ka is considered to make a significant contribution to the growth of HVZ.According to the result of tomography analysis by Tameguri et al. (2022), the HVZ is not a high Vp/Vs zone.Hence, the HVZ may represent the inactive and cooled remnant of a shallow magma reservoir associated with the 30 ka AT eruption.Alternatively, the HVZ may represent mafic cumulate material from fractional crystallization prior to the AT eruption.Regardless, the origin of the HVZ in relation to the Aira caldera eruption history requires further investigation.
A deep LVZ at ~ 15 km, beneath the Aira caldera, coincides with low seismicity and is located in the vicinity of the high Vp/Vs body derived by tomography analysis by Tameguri et al. 2022.Accordingly, we interpreted the deep LVZ as the current active magma reservoir.
We develop a new magma transportation model based on the distribution of the VT and LF earthquakes, the mapped geometry of the modelled pressure sources derived from the GNSS observed ground inflation (Hotta et al 2016), and our seismic transect velocity model.We found that the inclined HVZ plays an important role in the magma transportation path, especially from the center of the Aira caldera to the Sakurajima volcano.In addition, from the velocity model, we estimated that the Other symbols are same as those in Fig. 7 cylindrical magma reservoir associated with the Aira caldera formation at 30 ka was distributed in a depth range of 4-11 km at the center of the Aira caldera.
To our knowledge, this study is the first attempt to construct a quantitative caldera model of the Aira caldera, and further studies are expected to increase our understanding.

Fig. 2
Fig. 2 Map of the 195 km E-W transect constructed during the 2017/2018 seismic experiments.a Observation sites in Profile-A, -B, -C, -D, and -E are denoted by black dots.Onshore explosive sources (SP1-SP9) and vibroseis sources (VA and VC) are indicated by inverse triangles.Offshore airgun sources in Profile-B, -D, and -F are indicated by blue lines.The red broken straight line depicts the projection line used in travel-time tomography.b the profile of topographic heights at all stations and all seismic sources

Fig. 3
Fig. 3 Examples of the recorded sections for the onshore explosive sources: (a) SP1, (b) SP2, (c) SP3, (d) SP4, (e) SP5, (f) SP6, (g) SP7, (h) SP8, and (i) SP9.Each trace was digitally bandpass filtered (4-20 Hz) to suppress ground noise and normalized by its maximum amplitude.Reduction velocity was taken to be 6.0 km/s.Red and blue circles indicate the observed and calculated travel-times, respectively.The locations of the controlled explosive sources are shown in Fig. 2

Fig. 4
Fig. 4 Examples of the recorded sections for the offshore airgun sources: (a) AB1188, (b) AB1314, (c) AD2215, (d) AD2537, and (e) a map showing the airgun locations (red crosses) of AB1188, AB1314, AD2215 and AD2537.Each trace was digitally bandpass filtered (4-20 Hz) to suppress ground noise and normalized by its maximum amplitude.Reduction velocity was taken to be 6.0 km/s.Red and blue circles indicate the observed and calculated travel-times, respectively

Fig. 12
Fig. 12 Velocity cross sections in the 100-m cell model for the upper crust beneath the Aira caldera: (a) velocity distribution and (b) velocity perturbation distribution.The region surrounded by red contour lines in (a) indicates a velocity of 6.8 km/s for the high-velocity zone (HVZ).The upper LVZ indicates a slightly lower-velocity zone (5.6 km/s).Cross and plus symbols indicate the spherical pressure sources-A and -B estimated by Hotta et al. (2016) and Iguchi et al. (2013), respectively.Volcano-tectonic (VT) earthquakes (white circles) are also plotted.Other symbols are same as those in Fig. 7

Fig. 13
Fig. 13 Velocity cross sections for the 200-m cell model in the crust beneath the Aira caldera: (a) velocity distribution and (b) velocity perturbation distribution.Red contour lines in (a) indicate a velocity of 6.4 km/s; whereas the region surrounded by the white thick solid line denotes the deep low-velocity zone (LVZ).In (b), the region surrounded by the black thick solid line expresses the deep LVZ (same as that in (a)).In (a) and (b), the ellipse drawn by a broken line indicates the pressure source reported by Hickey et al. (2016).Cross and plus symbols indicate the spherical pressure sources-A and -B estimated by Hotta et al. (2016) and Iguchi et al. (2013), respectively.Volcano-tectonic (VT) (white circles) and low-frequency (LF) (yellow circles) earthquakes are also plotted.Other symbols are same as those in Fig. 6

Fig. 14
Fig. 14 Map depicting the distribution of earthquakes in and around the Aira caldera.Blue and red circles indicate volcano-tectonic (VT)and low-frequency (LF) earthquakes detected between 2000 and 2019 (data acquired from the Japanese-unified hypocenter catalog by the Japan Meteorological Agency).In (a), large and small circles drawn with a black broken line indicate the locations of the HVZ (high velocity zone) and the deep LVZ (low velocity zone) estimated in this study.Cross symbols denoted by "A", "K" and "M" are the pressure source-A, -K, and -M, respectively, as estimated byHotta et al. (2016).Red straight line indicates the projection profile in the tomography; whereas the red broken arrow line describes a theorized magma transport path.In (b), the green hatched zone delineates an aseismic region, while the deep LVZ is also presented.All other symbols are the same as those in (a)

Fig. 16
Fig. 16 Image of the magma reservoir associated with the AT eruption 30 ka. White and red rectangles indicate the HVZ (= the solidified magma reservoir with a volume of ~ 140 km 3 (this study)) and a pyroclastic ejecta of ~ 350 km 3 DRE from the AT eruption (Geshi et al. 2020), respectively.Other symbols are same as those in Fig.7