Bull Volcanol (2001) 63:326–344 DOI 10.1007/s004450100145

R E S E A R C H A RT I C L E

Y. Sudo · L. S. L. Kong

Three-dimensional seismic velocity structure beneath Aso Volcano, Kyushu, Japan

Received: 19 July 2000 / Accepted: 3 March 2001 / Published online: 22 June 2001 © Springer-Verlag 2001

Abstract Tomographic results for P- and S-wave velocity structure beneath the active Aso Volcano, Kyushu, Japan, using 800 well-recorded earthquakes and ten shots recorded by an eight-station seismic network, are presented. A 68% variance reduction was achieved upon simultaneous inversion for hypocenter and velocity structure. Well-resolved velocity anomalies associated with the active crater reveal heterogeneity up to 26% slower and 18% faster in P velocity, and up to 31% slower and 22% faster in S velocity, than the one-dimensional model. The largest anomaly is seen over the upper 11 km in the central and northern parts beneath the central cones. Two low-velocity regions are imaged. The first region, a 10×15-km region encompassing the upper 3 km centered near the caldera wall at Tateno Valley, is characterized by P velocities up to 19% slower (20% for S). The second low-velocity region is associated with the central cones and active magma conduit system at 6 km depth. Velocities as low as 4.3 km/s (up to 26%) in P and 2 km/s (31% slower) in S characterize the 7-km-wide volume. The magma chamber is roughly spherical in shape, centered at 6 km depth, flattens at 10 km depth, and is located between Mt. Kishima, Mt. Eboshi, and Mt. Naka, the present focus of magmatism. A sharp velocity contrast at the depth of 3 km, with high velocities to the southwest and lower velocities to the northeast, characterizes different abutting structures associated with the Oita-Kumamoto Tectonic Line.

Editorial responsibility: J.F. Lénat Y. Sudo (✉) Aso Volcanological Laboratory, Kyoto University, Choyo, Aso, Kumamoto, 869-1404, Japan e-mail: [email protected] L.S.L. Kong Hawaii Institute of Geophysics and Planetology, University of Hawaii, 2525 Correa Road, Honolulu, HI 96822, USA

Keywords Magma chamber · Aso Volcano · 3D-velocity structure · Low-velocity region · Magma conduit system · Tomography · Oita-Kumamoto · Tectonic Line

Introduction Aso Caldera, one of the largest in the world, is a 18×25-km elliptical caldera containing more than 17 volcanoes of ages up to 0.3 Ma. The current caldera was formed ca. 70,000–80,000 years ago during four large pyroclastic eruptions (Ono and Watanabe 1985), and consists of two caldera floors, Nango Dani and Aso Dani, averaging ca. 400–500 m in elevation, with the younger central volcanic cones ranging up to 1600 m in height. The caldera sits along the Oita-Kumamoto Tectonic Line (Fig. 1a), a northeast/southwest-trending geologic boundary separating Tertiary granitic rocks in the north from pre-Tertiary sedimentary and metamorphic rocks to the south, and Aso’s central cones lie along a NE–SW volcanic belt connecting the magmatically active Tsurumi, Kuju, and Unzen volcanoes of the central Kyushu volcanic field (Fig. 1a). Mount Naka, one of the youngest pyroclastic central cones, is active and characterized by Strombolian-type, phreatic eruptions of ash, mud, and scoriae ejecta occurring every few years; other cones comprise the central area include Mt. Kishima, Mt. Eboshi, and Mt. Taka (Fig. 1b). Volcanism is accompanied by continuous volcanic tremor activity located beneath the crater (Sudo 1978, 1984). Presently, steam emissions characterize crater activity during the period of volcanic quiescence. At Aso Volcano, numerous geothermal spring inns and resorts, with surface-water temperatures at some locations exceeding 100°C, are found throughout the active area. The highest-temperature hot spring spas are located at the Yunotani geothermal region west of Mt. Kishima and Mt. Eboshi. Temperatures below 1500 m depth are estimated to be in excess of 200°C in the Yunotani region (NEDO 1996). The local temperature gradient is estimated to be 34°C/100 m based on temperature–depth

327

Fig. 2 Aso region seismicity (RMS travel time residual ≤0.2 s) during the period from 1981 to 1989. Earthquakes (circles) were located by the eight-station AVL seismic network (triangles) using a 1D velocity model (Fig. 7) with station delays (Fig. 6a)

Fig. 1 Maps of a northern Kyushu Island and b Aso Caldera showing main tectonic structures. FF corresponds to Futagawa Fault which breaches the western caldera wall in Tateno Valley. Yunotani is the geothermal area. AVL is the center of Aso Volcanological Laboratory, Kyoto University. T.L. indicates Tectonic Line

data from drill holes, compared with an average of 4.5°C/100 m for Japan in general (Okubo and Shibuya 1993). Recently, numerous low-temperature springs have been found in the north caldera floor (Aso Dani) and in the south caldera (Nango Dani; Fig. 1b). These are indicative of sub-surface temperature anomalies and can be correlated with the observed seismic heterogeneity. In this paper, we investigate the P- and S-wave velocity structure in the Aso Caldera region using travel times from earthquakes during the 1981–1989 time period and shots to simultaneously relocate the earthquakes and determine the three-dimensional (3D) structure.

Background Seismicity in the Aso Caldera region trends southwest to northeast, shoaling over a 50-km distance from 15 to

20 km depth to less than 5 km depth beneath the caldera as shown in Fig. 2 (Sudo 1975b, 1981, 1993). In January 1975 many earthquakes with magnitudes up to 5.5 (maximum magnitude event recorded is 6.2) occurred in this region (Kubotera and Mitsunami 1980; Sudo 1975a). For a 13000-event data set observed during the period from 1964 to 1979, an average b-value of 1.08 was determined (Sudo 1981). Over the 1981–1989 time period, over 3500 earthquakes were recorded. During this period, seismicity was recorded around the western and northern caldera wall and outside the caldera (Sudo 1981, 1990). Little seismicity was observed beneath the active crater of Mt. Naka (Naka Dake crater), east of station ksm at Mt. Kishima (Fig. 2). Previous seismic structural results have shown that significant P- and S-waveform attenuation (Qp <100, Qs <50) and large P-wave travel-time residuals (27% >0.1 s) are characterized by a 10×15-km area beneath the central and northward extension of the caldera (Sudo 1988, 1991). Ray tracing suggests that the largest traveltime delay times are attributed to velocity anomalies located in a region at 6–9 km depth. Sudo (1991) studied a 24×24×18-km region and documented large travel-time delays beneath the central cones with 50% of the seismic rays passing through the 3-km blocks having delays of more than 0.3 s, 70% with delays 0.2–0.3 s, and 60% of the rays with greater than 0.1-s delays. The western wall of the caldera is breached by the Tateno Valley, a collapsed portion of the wall that is thought to have formed tectonically as active SWW–NEE (Futagawa Fault, FF in Fig. 1b) and NW–SE faults are found in this region (Watanabe et al. 1979; Watanabe 1984; Sudo 1993). Futagawa Fault is part of the southwest/northeast-trending Oita-Kumamoto

328

Tectonic Line, which is characterized by a steep gravity gradient and linear Bouguer anomaly that reflects the density differences between the granitic and metamorphic sedimentary rocks (Komazawa 1995). The gravity anomaly lineament is unidentifiable within Aso Caldera. The focal mechanism solutions from earthquakes occurring in this region are consistent with right lateral strikeslip motion in each of these regions (Sudo 1975b, 1993; Sudo et al. 1984). A positive elongate magnetic anomaly extending from the central cones to the west and coincident with the Futagawa Fault is consistent with the presence of a highly magnetized, 2D buried prism source dipping north from 10 km on the west to 3 km inside the caldera (Okubo and Shibuya 1993). The anomaly suggests that igneous rocks were intruded along the tectonic line. Existing temperature–depth data show a sudden change from low temperature outside the caldera to high temperatures inside, suggesting that the lower depth limits of the magnetic source material are temperaturedependent and related to the presence of volcanic activity within the caldera (NEDO 1996). Gravity modeling suggests that the presence of a 5-km-diameter, cylindrical block of low-density material (5% lower) extending 6–17 km below sea level beneath the northern slope of the central cones (Komazawa 1995), is consistent with heat flow measurements which suggest that partial melting of the basaltic rocks occurs at a depth of ca. 15 km (Ehara 1984). The gravity data also indicate the presence of five local low Bouguer gravity anomalies within the caldera, generally corresponding to the regions of Nango Dani, the northern slope of the central cones, Aso Dani, and towns of Uchinomaki, and Ichinomiya; the lows are bordered by a steep gradient inside the caldera rim (Komazawa 1995). The flat-bottomed local anomalies were attributed to lower-density, shallow sedimentary material or lake deposits, formed by the erosion, collapse, and recession of the central cones and caldera rim, and are consistent with seismic reflection and seismic microtremor observations that indicate the presence of several hundred meters (up to 400 m) of sediments beneath Aso Dani (Kubotera and Otsuka 1971; Watanabe et al. 1986) and Nango Dani (Tsutsui et al. 1997).

Data set In order to investigate the sizes of velocity anomalies associated with Aso Caldera, we carried out tomographic inversions using P- and S-wave travel-time residuals from well-recorded earthquakes occurring between 1981 and 1989. During this period, over 3500 earthquakes were recorded, with seismicity trending southwest to northeast. Earthquake P- and S-wave data were available from six to eight stations for over 1500 earthquakes during the 1981–1989 time period. These data were recorded by an eight-station (three-component, 1-Hz seismometers) network which encompassed the caldera and the active Naka Dake crater (Fig. 2). The network was oper-

Fig. 3 Tomography data set consisting of 800 earthquakes (RMS travel time residual ≤0.2 s and number of readings ≥10; open circles) located by AVL network (triangles), and 10 shots (four explosions, stippled diamonds, recorded by open-diamond stations, and 6 shots, stippled squares, recorded by open-square stations). Box delineates tomographic inversion region

ated by Aso Volcanological Laboratory (AVL) of Kyoto University as part of their geophysical monitoring of the active volcano. Data from the network are radio- and telephone-telemetered in real time to AVL where they are digitized at 200 Hz and stored on magnetic tapes. The tomography data set consisted of 800 well-recorded earthquakes covering a 40×47-km area (0–16.1 km depth range, 8.8 km average depth; Fig. 3). These events were located with more than ten readings (at least five P and S readings for each event) and each had final overall root mean square (RMS) travel-time residuals of less than 0.20 s (0.07 s average, 0.02 s standard deviation). Arrival times could be read to within 0.005 s accuracy from the digitized data. The average formal 1-sigma horizontal and vertical errors were 0.15±0.06 and 0.30±0.14 km, respectively. To help resolve the uppermost crust, ten quarry-

329

Fig. 5 Summary of Vp/Vs. An overall Vp/Vs ratio of 1.704 (standard deviation of 0.0078), taken as the average of the individual station Vp/Vs ratio, was used to locate all the earthquakes

Fig. 4 Examples of linear regression fits to P- and S-wave traveltime differences for a stations avl and b mak from 372 wellrecorded earthquakes having seven or eight P- or S-arrival times. The slope of the line provides an estimate of the Vp/Vs ratio

blast shots were included. Shot origin times were calculated from the travel times through the air as recorded by a tripartite array (distance span of 600 m) located at AVL. The earthquakes were located by first estimating a Vp/Vs ratio using the best-recorded events for the network. A modified Wadati diagram (Chatterjee et al. 1985) was constructed by plotting the P-wave travel time differences vs S-wave travel time differences for each station from 372 earthquakes having seven or eight P- or S-arrival times (nearly 2800 data points for each station) and linear regression fits were determined for each station (Fig. 4). The distribution of sources for the Vp/Vs ratio analysis is very similar to the distribution of sources for the tomography analysis. The Vp/Vs ratio corresponds to the slope of the best-fitting line of all the data, and the estimated error is the standard deviation in the slope fit (Francis 1976). The individual Vp/Vs ratios ranged from 1.693 for station nbr to 1.715 for station hdk. The Vp/Vs ratio for station ksm, near the active Naka Dake crater where volcanic activity is centered, was 1.701. Because the Vp/Vs ratio-range was very small, we used an overall Vp/Vs ratio of 1.704, taken as the average of all the station Vp/Vs ratios, to locate all

Fig. 6 a Station delays determined using 372 well-recorded earthquakes. The final delays are those which produce near-zero average travel-time residuals with the smallest standard deviation about the average. b P-wave travel-time residuals from these earthquakes. c S-wave travel-time residuals from these earthquakes

the earthquakes. The standard deviation about the average was 0.0078 (Fig. 5). The earthquakes were located using the Hypoinverse earthquake location program (Klein 1978, 1989), and a 1D P-wave velocity model previously determined for Aso Volcano (Sudo 1991) with station delays (Figs. 6a, 7; Tables 1, 2). Station delays were included to account for near-station velocity heterogeneity and determined using the 372 well-recorded events and requiring the average

330

hdk, mak, mkn, and tat. The delays ranged from –0.17 to 0.17 s (Fig. 6a), and were slightly larger (by <0.05 s) in magnitude than the delays assuming only station elevation differences. The two exceptions were station mkn, whose delay became more negative (by 0.10 s, implying actual velocities greater than the model) and station nbr, whose delay became positive (by 0.15 s, implying velocities slower than the model). The delay for station trm (0.22 s), located far outside the study area (90 km northeast of station avl), corresponds to the elevation delay and was not recalculated because no earthquake sources using this station were used. Excluding station trm, the largest positive delay was for station ksm (0.17 s), whereas the largest negative delay was beneath station tat, located 10 km northwest of Aso Caldera.

Fig. 7 One-dimensional velocity models. Solid line is plane-layered model used to locate the earthquakes. Dashed line is the starting model for the tomographic inversions

Table 1 One-dimensional velocity model(Vp/Vs=1.704)

Velocity (km/s)

Depth (km)

3.80 4.40 4.80 4.95 5.10 5.30 5.50 5.70 5.90 6.20 6.60 7.00 7.40 7.70

0.00 0.40 0.75 1.40 1.90 2.65 3.75 4.75 5.75 6.95 17.00 27.00 37.00 47.00

Table 2 One-dimensional station delays Station

Location

avl mkn nbr tat mak hdk mgr trm ksm

32° 52.93′N 33° 00.55′N 33° 06.45′N 33° 00.79′N 32° 55.04′N 32° 47.27′N 32° 47.78′N 33° 37.00′N 32° 53.55′N

131° 00.54′E 131° 11.92′E 131° 03.20′E 130° 53.94′E 130° 55.68′E 130° 54.52′E 131° 05.58′E 131° 25.91′E 131° 04.04′E

Elevation (m)

Delay (s)

550 670 450 300 380 475 700 1356 1105

0.039 –0.079 0.110 –0.175 –0.101 –0.070 –0.018 0.217 0.173

travel-time residuals at each station to be near zero. The iterative calculation of station delays reduced the scatter in the P- and S-wave travel-time residual about the average by 32 and 21%, respectively, over runs which used elevation delays only (Fig. 6b, c). The delays indicate that velocities less than the 1D model are located beneath stations ksm and nbr, whereas velocities greater than the model characterize the region below stations

Analysis method P- and S-wave travel-time residuals were inverted for velocity perturbations and hypocenters within a 3D volume using the damped-least-squares method developed by Thurber (1981, 1983, 1993), later modified to include S waves by Eberhart-Phillips (1993) and Thurber (1993). The method has been used successfully in numerous volcanic settings, including Hawaii (Thurber 1987), Iceland (Arnott and Foulger 1994; Foulger et al. 1995), and the Geysers geothermal area (Eberhart-Phillips 1986). Parameter separation (Pavlis and Booker 1980) of the hypocentral and velocity components decouples the problem, greatly reducing computational time by permitting the two components to be treated independently. Tests indicate that the simultaneous inversion for hypocenter and structure, rather than the use of hypocenters as fixed sources, significantly improves the determination of the true velocity model and hypocentral parameters (Thurber 1981; Kissling 1988). Each inversion involves the relocation of earthquakes using the better-known structural model followed by the re-determination of velocity perturbations. The iterative process is repeated as long as model improvement is judged significant at the 95% confidence level, as determined by an f-test on the reduction in variance values. Travel times were calculated using a computationally efficient and fast approximate ray tracing algorithm (Thurber 1981), which calculates accurate travel times over local distances (<50 km; Evans et al. 1994). Velocity adjustments for each iteration were multiplied by a damping factor to suppress large model changes for nodes with near-singular values. The damping parameter was determined by performing a series of one-iteration inversions using different damping values and the specific source-receiver configuration and model parameterization, and choosing that value which minimized the data variance while only moderately increasing model variance. Velocity is parameterized by a discrete set of velocities on a 3D, arbitrarily-spaced grid, and velocity at any position is defined as a weighted average of the surrounding eight nodes.

331

Finally, we solved for Vp/Vs within the velocity grid volume by including the (S–P) travel times in the inversion. These data were given a 0.75 weighting relative to the P wave times because of the greater uncertainty in picking S waves. A smaller Vp/Vs damping value of 3 s2/km was used for the inversion since the S velocity field is in general much more uniform and slowly varying than the P velocity field. The S wave structure is then calculated from the P and Vp/Vs structures determined in the inversion. The target volume encompassed a 25×40-km area in and around Aso Caldera. Figure 8 shows the tomography model parameterization for the preferred inversion for P-wave velocity structure. The plane-layered velocity model, parameterized as a 3D volume in Cartesian coordinates, has nodes 5 km apart in X and Y, and 3 km apart in Z. The starting velocity model used to locate the earthquakes is shown in Fig. 7. In determining the best model parameterization, we examined grid spacings from 3 to 8 km in the horizontal, and 2 to 3 km in the vertical, direction. The spacings were based on the distribution of sources and receivers and a priori knowledge of the expected scale of heterogeneity. A trade-off between the scale of resolved heterogeneity and parameter resolution is a well-known effect, and our preferred model, 5 km in X and Y and 3 km in Z, attempts to extract the most structural information at the finest scale over the greatest volume given the a priori information. Fig. 8 Velocity model parameterization for inversions. Nodes (squares) were placed 5 km apart in X and Y, and every 3 km in Z (depth). Triangles show the locations of the AVL network

In the initial inversions, earthquake sources were held fixed throughout to obtain a first-order estimate of the magnitude and dimensions of velocity heterogeneity in the Aso region, and to determine the optimal velocity inversion parameters (Kong and Sudo 1993). After a preliminary velocity structure was determined, we performed inversions solving for both earthquake location and P-wave structure simultaneously, while holding the Vp/Vs ratio constant throughout the inversion. A Vp/Vs ratio of 1.704 determined from a modified Wadati diagram was used. The Vp damping factor was chosen after running a series of inversions with values ranging from 0.5 to 200 s2/km. Our tests indicated that the use of small values (D <1 s2/km), approximating the ratio on the data to model variance, results in unjustified oscillations in velocity which approximate the grid spacing, whereas large values (D >50 s2/km) overcompensate by allowing little or no perturbation to the initial model. For the final solution, a Vp damping value of 15 s2/km was used. For other inversions that tested model robustness, damping ranged from 6 to 20 s2/km. The final velocity model (inversion 1 in Table 2) corresponds to the inversion using the best estimates of the true hypocenters, and results in a 68% reduction in the travel-time residual variance between the initial and final velocity models.

Solution quality To investigate the quality of the determined model, other inversions were carried out in which the starting model, the numbers of earthquakes, and the nodal grid spacing were varied (Table 3), as well as inversions which included the S wave times (using instead the equivalent P-wave times from the Vp/Vs ratio). Our tests indicate that the model solution is robust, as the patterns and magnitudes of velocity heterogeneity, and especially that of the low-velocity zone beneath the active crater with Aso Caldera, are reproduced in all the models. We examined the effect of using a smaller, higherquality, data set to determine the 3D structure. Inversion 7 used 405 earthquakes, each with at least 14 arrival times (seven P- and seven S-wave times) and hypocenters similarly distributed as in inversion 5 (800 earthquakes). No significant differences in the magnitude or patterns of velocity heterogeneity were observed, although slightly higher RMS travel-time residuals and variances were observed because the smaller data set provided less ray coverage than the larger data set. The resolution of a given node is taken as the size of the diagonal component of its resolution matrix. It has been shown that information on the velocity structure can be recovered from nodes where the diagonal resolution is low if the node is compact (Stauber et al. 1988; Toomey and Foulger 1989), or if the averaging vectors or resolving kernels (rows of R) are non-zero only near

332 Table 3 Summary of tomographic inversion runs Run

nx, ny, nz dx, dy, dz (km)

#Eqs #Vel (km)

Traceb P, S

%LIPc P, S

Solution typed Data type

#Iter Dampe

rmssf rmsf g

vars varf

%rms Redg %var Redg

Inversion 1

6, 9, 5 5, 5, 3 5, 6, 5 7, 7, 3 10, 13, 5 3, 3, 3 6, 8, 7 5, 5, 2 6, 8, 5 5, 5, 3 6, 8, 5 5, 5, 3 6, 8, 5 5, 5, 3 6, 8, 5 5, 5, 3

800 270 800 150 800 650 800 336 800 240 800 240 405 240 867 240

40 86 28 64 90 129 60 102 64

15 32 19 43 14 20 18 30 27

64

27

42

18

45

19

P, Vp/Vs P, S–P P, Vp/Vs P, S–P P, Vp/Vs P, S–P P, Vp/Vs P, S–P P P, S–P P P, S–P P P, S–P P P

6 15 9 20 9 5 9 8 3 6 5 6 3 10 3 15

0.1369 0.0774 0.1369 0.0781 0.1369 0.0723 0.1369 0.0744 0.1369 0.0908 0.1249 0.0772 0.1374 0.0919 0.1378 0.0911

0.0188 0.0060 0.0190 0.0061 0.0188 0.0052 0.019 0.0055 0.0188 0.0082 0.0156 0.0060 0.0189 0.0084 0.0190 0.0083

43 68 43 67 47 72 47 72 34 56 38 62 33 55 34 56

Inversion 2 Inversion 3 Inversion 4 Inversion 5h Inversion 6h Inversion 7h, i Inversion 8

a No. of P- or S-velocity parameters (m) b Trace, K c Percentage of linearly independent parameters d Inversion solves for P, or P and Vp/Vs e Damping factor used in velocity inversion

(k/m)

f Weighted by data contribution to solution g Percent reduction h Vp/Vs nodal values fixed at 1.704 i Earthquakes with at least 13 weighted readings

Fig. 9 Complete resolution matrices for selected P and Vp/Vs nodes. The complete matrix was examined to evaluate the compactness of each velocity node. Nodes with nodal resolution values as low as 0.20 are compact since the magnitude of resolution at adjacent nodes is negligibly small (e.g., node 58), indicating that velocity information can be reliably extracted from these nodes

the principal diagonal (Wiggins 1972). For a perfect data set, the resolution matrix, R, corresponds to the identity matrix. In reality, however, the resolving kernels never take on a delta-like appearance, and we must examine the complete resolution matrix to see whether nodes with low resolution also provide important velocity information (Fig. 9). Resolution is affected by the data kernel, and ray density can be quantified by calculating the derivative weight sum (DWS; Toomey and Foulger 1989; Arnott and Foulger 1994). The DWS gives a more realistic assessment of ray density than a simple count of the number of rays passing near to a node, because each ray is weighted by its distance to the node. The size of the DWS is dependent on the incremental arc length so that smaller step lengths yield larger values of DWS. Wellconstrained regions are those where the DWS is smooth

and slowly varying. When many rays sample the node in similar directions and at similar depths, however, a smearing in the velocity variations can result because the node is not compact despite a large DWS. To assess compactness, we examined the complete resolution matrix to determine what threshold of resolution and DWS values would result in nodes with negligibly small, neighboring resolution values (Figs. 9, 10). In general, the side lobes of nodes with resolution greater than 0.2 were small. Based on this analysis, we interpret regions where the DWS is greater than 1500 km (resolution >0.20) as regions where velocity can be constrained. We compared tomography results from runs in which the grid spacing was varied from 3 to 8 km in the X and Y directions, and 2 to 4 km in Z direction. Final RMS residuals increased with increasing grid spacing

333

Because all of the earthquake sources have both P- and S-wave readings, we were able to invert for the S-wave structure as well. Our tests comparing P- and both P- and S-arrival times to solve for P or both P and Vp/Vs showed that the inclusion of the P- and S-wave data resulted in lower final RMS values and greater reductions in variance (inversions 1, 7, and 8; Table 3), and significantly improved the resolution of the P-wave structure, thus permitting a greater volume of the study area to be constrained. Inversions using only P-wave data resolved a smaller volume, but the patterns and magnitude of velocity heterogeneity were similar to the results using both P- and S-wave data.

Fig. 10 Resolution (value of diagonal element) for inversion 1; stippled regions indicate where the DWS ≥1500 km. As explained in the text, we use the DWS to assess solution robustness. A comparison of the resolution and DWS shows similarity in the shapes in the contours, and a rough correspondence of the 0.2–0.4 resolution contour with the 1500-value DWS contour

(0.0723–0.0781 s from 3 to 7 km spacing; Table 3). Our tests also show that the trace (k) decreases with increasing nodal spacing (decreasing number of velocity parameters), whereas the percentage of linearly independent velocity parameters (%LIP=k/m) increases. The trace is calculated as the sum of the diagonal elements of the resolution matrix, and for a finite number of m velocity parameters, provides an estimate of the number of degrees of freedom of the system (Wiggins 1972). For our simulations, k decreased from 90 to 28 for P (129 to 64 for S), as node spacing increased from 3 to 7 km in the horizontal directions (Z spacing held constant). At the same time, the %LIP increased from 14 to 19% (20–43% for S). Because k remains fixed while m increases, a decrease in nodal spacing and thus a decrease in %LIP also effectively decreases the resolution values of the velocity parameters. In addition, for an underdetermined inverse problem, the resolution matrix will in general have nonzero off-diagonal elements, which need to be examined for compactness. We tested whether the use of an initial 3D model resulted in significantly better results. The use of a 3D model reduced the RMS residual by less than 5% (e.g., inversions 5 and 6; Table 2), and resulted in greater heterogeneity, but included large, spurious perturbations in poorly resolved regions. Because of this, we opted to use a homogeneous plane-layered model as the initial starting velocity model for our final inversion. The simple, 1D model has the advantage of not introducing a priori biases to the structure, which could be incorrectly biased. Furthermore, in the final assessment, the good correlation of the structure with known geologic and geophysical anomalies gave us confidence that the results are reliable and that the resolved portions of the inversion solution provide an accurate image of the true velocity structure.

Results The results of 3D velocity structure of the Aso Caldera region are shown in series of vertical and horizontal sections from Figs. 11, 12, 13, and 14. An examination of the velocity perturbations for well-resolved regions clearly reveals velocity anomalies that are associated with the active volcano. These results were obtained using 10298 (5453 P and 4845 S) travel times from 800 earthquakes and ten shots to determine estimates of 270 velocity parameters; of these, 40% of the parameters had good resolution. Velocity was well resolved over much of the study volume from 3 to 9 km depth, with regions proximal to the stations resolved at the 0-km depth velocity nodes. The average travel time residual, weighted by the importance of each travel time observation, was reduced by 43% to 0.777 s and variance reduced by 68% to 0.006 s2 after simultaneous inversion for hypocenter and velocity structure (inversion 1; Table 3). The final variance approximates the estimated maximum data error (0.005). The average of the individual time residuals was reduced to near zero, and initial travel-time residuals biases (as indicated by the ~0.1 s standard deviation in P times) caused by near-surface heterogeneity were removed after simultaneous inversion (Fig. 15). The mean formal model uncertainty, derived from the covariance matrix as a standard error, was 0.05 km/s and 0.04 for resolvable P-wave velocity and Vp/Vs parameters, respectively. The absolute mean velocity error is thought to be approximately two times the mean standard error (Thurber 1981), or for this study, 0.10 km/s for P-wave velocity and 0.07 for Vp/Vs or S-wave velocity. Velocity heterogeneity up to 26% slower and 18% faster in P velocity (range from –1.6 to 0.8 km/s), and up to 31% slower and 22% faster in S velocity (range from –1.1 to 0.6 km/s), than the 1D model were imaged. An average perturbation of 0.4±11% (median 5%) for P and 1.2±13% (median 7%) for S was found, indicating that our 1D starting model constitutes a good representation of the average velocity field in the study region. Vp/Vs perturbations ranged from –8 to 9%, with a mean of –0.3±5%. P- and S-velocity anomalies were well correlated, with R value of 0.92. Few events occurred at depths greater than 12 km (Fig. 3), resulting in poor res-

334 Fig. 11 Horizontal P-wave sections for the study region. Velocity resolution over the entire study area is very good for depth slices at 3, 6, and 9 km. Triangles show the locations of the seismic stations used to locate the earthquakes (station names shown on Z=0 km and Z=12-km sections). Heavy solid line delineates that outline of the caldera. Stippled areas show the regions where velocity is well resolved (DWS ≥1500). Velocity and DWS are contoured at 0.5-km/s and 1000-km intervals, respectively. Perturbations are contoured at 5% intervals with solid lines indicating positive values, or velocities faster than the average model, and dashed lines indicating negative values, or velocities slower than the average model. Main lower-velocity anomalous regions are indicated as LA, LB, LC, LD, and LE, and higher regions as HA, HB, HC, HD, HE, HF, and HG

olution of 12-km depth nodes except in a small region west of the caldera where no significant anomaly was observed. When interpreting the velocity structure, we utilized the solutions using finer grid spacings and different node locations to estimate the boundaries and dimensions of the constrained anomalies. Use of a finer grid spacing can help to constrain model details, because the travel time residual is ascribed to a smaller averaging volume. We assume a velocity anomaly to be significant if it is characterized by a DWS greater than 1500, has a

magnitude greater than 10%, and is observed in over more than one node in the P- and/or S-wave structure. Table 4 lists the velocity anomalous regions. The largest anomaly in the study region is seen over the upper 11 km over the central and northern half beneath the central cones of Aso Caldera. Horizontal P- and S-wave velocity sections at 3, 6, and 9 km depth (Figs. 11, 13), and the vertical section at X=3 km (Fig. 12, 14), clearly image two regions of low velocity. The first region encompasses the upper 3 km of crust and

335 Fig. 11 (continued)

is centered 2 km northwest of station avl at (X, Y)= (–1, 0), near the base of the collapsed caldera wall at Tateno Valley (named LA). The 10×15-km region covers the central north portion of the caldera (at Z=0 km, X=–7 to 3 km, Y=–4 to 11 km) and is characterized by P velocities up to 22% slower (27% for S) than the 1D model. The second low-velocity region is associated with the active central cones (LD). The lowest velocities were centered at 6 km depth at (3, 4), 1 km south of station ksm (Figs. 11, 12). Velocities as low as 4.3 km/s (up to 26%) in P and 2 km/s (31% slower) in S were found. The anomalous low-velocity volume covers a 7×6×7-km region and is roughly spherical in shape extending downwards from Z=2 to 11 km, and extending from 0 to 7 km and 1 to 7 km in the X and Y directions, respectively (X=3 km; Fig. 12). The upper 8 km of the low-velocity zone (LVZ) is well resolved, but while that from 8 to 11 km is not as well constrained, the patterns of velocity heterogeneity suggest that it extends to 11-km depth and its shape becomes more elliptical in the northwest direction and flattens at depth by approximately 5 km. Velocity anomalies associated with heterogeneity beneath stations hdk, mkn, nbr, and tat, in addition to stations avl and ksm within the caldera (described above), were found (Figs. 11, 13). Velocities greater than the average model characterize the upper several kilometers be-

neath stations hdk, mkn, and tat, whereas slower velocities characterize the region beneath stations nbr, avl, and ksm. The large anomaly beneath station tat (up to 18% in P, 22% in S) to the northwest of Aso Caldera extends to 6 km depth and trends downward to the east ending outside the caldera wall (HA). P velocities up to 17% (14% in S) slower than the 1D model are seen in both the 0and 3-km depth slices beneath station nbr north of Aso Caldera (LB). The relative magnitude and sign of the anomalies beneath the stations are consistent with the character of heterogeneity estimated from the determination of station delays for the 1D model (Fig. 6). At 3 km depth several kilometers northeast of station mak, a small pocket (≤5 km in width, ≤3 km in thickness) of higher velocities (up to 14% in P, 12% in S) is imaged (HB). Within the caldera, one pocket of higher velocities (up to 14% in P, 15% in S) is centered at (X, Y, Z)=(3, –2, 3), approximately 3 km south of station avl in the southwest portion of the caldera, and deepens to 6 km depth beneath the caldera wall 10 km to the southwest (HC). The second pocket is located at the southeastern edge of the caldera (and our study region) at least 5 km east and south of station ksm (HD). Velocities up to 11% faster are observed over a 5- to 10-km region in the Y direction and 8–10 km in the X direction from 2 to 6 km depth. At 9-km depth, higher velocities

336

Fig. 12 Vertical P-velocity sections for the study region. See Fig. 11 for map symbol explanations. Station names shown on X=8-km section

are found centered at (–2, –10, 9), 7 km southwest of the breached caldera wall in Tateno Valley (HF). The anomaly is characterized by velocities up to 11% faster in P, and up to 12% faster in S, and extends from the mouth of Tateno Valley where it meets the caldera wall over a 10×5-km southwest-trending region.

337 Fig. 13 Horizontal S-velocity sections for the study region. See Fig. 11 for map symbol explanations

Smaller-magnitude velocity anomalies were imaged in which the S-wave velocity perturbations were greater than 10%, but P perturbations were approximately 50% smaller than the S anomalies. At 9 km depth, an elliptical region of velocities up to 14% faster in S and 7% faster in P extends from inside the northern part of the caldera (–7, 9, 9), northward to approximately 7 km outside the wall (HG). Small anomalies (≤5 km in width, ≤3 km in thickness) of faster than average velocities (up to 11% in S, 5% in P) are centered at (–12, –6, 3), 9 km west of the caldera wall and 5 km west of station mak (HE), and slower-than-average velocities are found at

(–2, 14, 3), in the northern part of the caldera at the base of the wall (up to 11% in S, 7% in P; LC), and at (–2, –6, 6), in Tateno Valley 3 km outside of the wall (up to 11% in S, 4% in P; LE).

Hypocentral relocation In our study, the simultaneous inversion of the earthquake data set for hypocenter and 3D velocity structure resulted in relatively small changes in hypocenter from the starting locations estimated from a 1D model with

338 Fig. 13 (continued)

Table 4 List of velocity anomalous regions Low-velocity region Name Location

Center

Velocity

LAa LB LCb LDc LE

(–1, 0, 0) (–19, 19, 0) (–2, 14, 3) (3, 4, 6) (–2, –6, 6)

P –22%; S –27% 10×15×3 P –17%; S –14% 7×5×3 P –7%; S –11% 3×3×3 P –26%; S –31% 7×6×7 P –4%; S –11% 3×3×2

Name Location

Center

Velocity

Dimension (km)

HA HB HCb HDb HE HF HG

(–19, 4, 0) (–7, -2, 3) (3, -2, 3) (10, 8, 3) (–12, -6, 3) (–2, -10, 9) (–7, 9, 9)

P +18%; S +22% P +14%; S +12% P +14%; S +15% P +11%; S +11% P +5%; S +11% P +11%; S +12% P +7%; S +14%

5×6×5 3×3×3 5×3×3 3×5×4 2×2×2 3×3×2 3×3×2

2 km NW of avl near nbr NNE of caldera 1 km S of ksm Tateno Valley

Dimension (km)

High-velocity region

near tat 3 km E of mak 3 km S of avl SE of caldera 5 km W of mak 6 km N of hdk No. of caldera

a Corresponds b Corresponds c Corresponds

to the low-gravity region Oita-Kumamoto Tectonic Line the magma chamber

station delays (Fig. 16). The average differences in epicenter, depth, and origin time were 0.510±.34 km, 1.06±0.64 km, and –0.08±0.08 s, respectively (Fig. 17). For our 800-earthquake data set, epicenter, depth, and origin-time differences ranged up to 2.03 km, 3.97 km, and 0.10 s, respectively. The average depth of the relocated data set was 7.84±2.38 km, with source depths from 0 to 15.4 km. Because of the slower shallow velocities within the much of the caldera, there was a systematic shallowing in earthquake depth upon relocation. The translation implies that the relative locations are better known than their absolute locations, and is a manifestation of the non-uniqueness of earthquake location algorithms, especially in the presence of 3D heterogeneity. However, because we are seeking to extract velocity variations at wavelengths much larger than the differences in earthquake hypocenter, the absolute hypocentral estimates are not critical and the observed shifts in epicenter and depth do not affect the patterns of heterogeneity.

Discussion Low-velocity region LD In previous work to image the upper crustal structure beneath Aso Caldera, Sudo (1991) found the presence of

339

Fig. 14 Vertical S-velocity sections for the study region. See Fig. 11 for map symbol explanations

attenuating material in the central region beneath the caldera that extended to the north of the central cones at depths from approximately 6 to 9 km. A Qp value of the material was approximately 100.

In this study we confirm the presence of a lowvelocity region of anomalous material beneath the central cones, but do not see evidence for low velocities north of the central cones. The tomographic images show that the low velocities (up to 31%) cover a 7-km wide region extending from approximately 3 to 11 km depth (Z=3–9 km, Y=4 km, X=3 km sections; LD;

340 Fig. 15 Travel-time residuals, relative to the final velocity model, after simultaneous hypocentral relocation and velocity structural inversion. The 3D model results in smaller average for the residuals and a decrease in standard deviation

Fig. 17 Summary of the differences in epicenter (DX), depth (DZ), and origin time (DOT) for 800 earthquakes after simultaneous inversion for hypocenter and velocity structure. Comparisons are made between tomography inversion 1 (nodal spacing 5 km in X and Y, and 3 km in Z) and hypocenters located with the 1D model with station corrections

Fig. 16 Relocated hypocenters of 800 earthquakes using the 3D velocity model determined in this study. Circles represent the relocated location, and a line connects the original location using the 1D model with station corrections (see Fig. 3 and Table 1) with the relocated location

Figs. 11, 12, 13, 14) and are centered 1 km south of station ksm. This low-velocity region (LD) appears to be roughly spherical in shape, centered at 6 km depth, and not directly beneath the active crater of Mt. Naka, the present focus of magmatic activity, but rather between the cen-

341

tral cones of Mt. Kishima, Mt. Eboshi, and Mt. Naka. Although not as well resolved, the velocity contours suggest that this region may flatten at 10 km depth toward the northwest. These features appear to be robust and insensitive to the grid spacing in the X, Y, and Z directions we use for our tomographic inversions. This region is also marked by little or no recorded seismicity, implying that the region cannot support brittle failure. Because of the aseismicity, we are unable to further constrain details of the region with only the earthquake sources; however, this lack of seismicity also provides further evidence that the region is the center of active volcanism. Numerous measurements have been made of the ground deformation and gravity near the LD region. In the Aso Volcano area, leveling surveys have been repeated since 1937 along the route (17 km in length) from the north foot of the central cones to the active crater of Mt. Naka (Yoshikawa 1954; Kikuchi and Ono 1985). The change in the vicinity of the Mt. Naka crater is not significant (ca. 0.5 cm), but a greater change is seen to the south of Mt. Kishima which may correspond to the LD region (maximum ca. 4.5 cm). Precise gravity measurements have also been made repeatedly along the above-mentioned leveling route since 1964 (Kubotera and Satomura 1985; Sudo and Yoshikawa 1995). The gravity change of each leveling point from 1971 to the present is seen to be small; however, large changes were seen at the points near Mt. Kishima. The change was approximately 80 µgal between 1979 (for this year and 1989–1990, the volcano was very active) and 1994 (this year was quiet). The position where the gravity change is seen to be largest is not near the active crater of Mt. Naka but the region approximately 3 km to the west of Mt. Naka (LD). At Aso Volcano the volcanic tremor observations were made with broad-band seismometers. Even during the volcanic quiet period, many isolated very long-period tremors are recorded (Kaneshima et al. 1996; Yamamoto et al. 1999a, 1999b). The period of tremor is ca. 15 s as a fundamental mode. The source of this long-period tremor is not beneath the active crater of Mt. Naka but to the west, 1 km from the crater at a depth of ca. 1 km. This fact suggests that the volcanic gas, which plays an important role in the occurrence of the tremor, is supplied from the west of the crater (LD). Moreover, at Aso Volcano, sub-surface temperatures were measured at a depth of 1 km below sea level in drill holes inside and outside of the western part of caldera, at the Yunotani geothermal area near the LD region and the Futagawa Fault area (NEDO 1996). The results show temperatures in excess of 190 and 78°C, respectively. The contour maps of underground temperatures at 50, 100, 150, and 200°C indicate to be the concentric circles centered at the south point of Mt. Kishima near the LD region. The magnetotelluric electrical resistivity measurements in the same part of caldera also reveal pockets of high resistivity that are indicative of hot rock associated

with the geothermal activity at the Yunotani area (near the LD region; NEDO 1996). Therefore, from the above-mentioned results, we interpret this low-velocity region LD to be associated with the magma chamber at Aso Volcano, as would be expected to be a partially molten magma system. The velocity of P wave is decreased by 10–40% when a material with 10% partial melting, whereas the S-wave velocity decreases by more than 20% (Mavko 1980; Sato et al. 1989). In the present case for the LD region, the pressure situation is different from the cases of Mavko (1980) and Sato et al. (1989), but it is possible that the decrease of velocity in the LD region is due to more than 10% of partial melting and the LD region is a magma chamber. Figure 18 indicates the horizontal and vertical sections of Vp/Vs ratio perturbations at the depths of Z=0, 3, 6, and 9 km and at the lines of X=–3, 2, and 7 km, respectively. In the LD region, Vp/Vs ratio perturbations are only slightly larger than normal. Numerous studies of the changes of P- and S-wave velocities or Vp/Vs ratios (Aster and Meyer 1989; Chiarabba et al. 1995; Wu and Lees 1999) have concluded that both P- and S-wave velocities decrease with increasing temperature for a uniform composition, but at different rates. Vs generally decreases faster, so the Vp/Vs ratio also increases with increasing temperature. Our results of lower velocities and high attenuation (Sudo 1991) in the LD region therefore indicate the presence of high-temperature material. Low-velocity-region LA A region of low velocities characterizes much of the shallowmost crust in the central and northern portions (Z=0–3 km; LA; Figs. 11, 13) coincident with local negative Bouguer gravity anomalies associated with the northern central cones (Komazawa 1995), Nango Dani (the southern caldera floor), and Aso Dani (the northern caldera floor), of Aso Caldera, and also consistent with magnetic modeling which reveals a magnetic source in the upper 3 km of the crust just west of Mt. Eboshi within Aso Caldera (Okubo and Shibuya 1993). Low velocities often characterize tectonically disrupted areas, such as fault zones, where the crust becomes broken up by the fracturing allowing water to percolate downward and become heated by the magma source. The low velocities are probably also due to the presence of several hundred meters of low-velocity, low-density, lake sediments which filled crevasses after collapse episodes in Aso Dani and the northern central cone region (Kubotera and Otsuka 1971; Watanabe et al. 1986). The largest anomalies (up to 27%; LA) are centered near station avl, and extend several kilometers outside the western caldera wall just north of Tateno Valley. Detailed gravity data in and around the Aso volcanic area show that Aso Caldera is a piston/cylinder-type structure rather than a funnel-shaped structure with a single low anomaly (Komazawa 1995), as was previously believed (e.g., Kubotera et al. 1969). A piston/cylinder-

342

Other velocity-anomalous regions LC, HC, HD Tateno Valley, lying along the active SWW/NEE-trending, northward-dipping Futagawa Fault of Oita-Kumamoto Tectonic Line, is the only portion of the Aso caldera wall that has collapsed. Gravity anomaly data delineate a normal fault through the gap and magnetic data have been modeled by a highly magnetized prism source dipping north beneath the fault (Okubo and Shibuya 1993), consistent with normal faulting mechanisms found for earthquakes occurring beneath the valley (Sudo 1993). Figures 11 and 13 (at the depth of Z=3 km, Y=3 km; LC and HC, HD) show a marked contrast in velocities, with the southwest region showing high velocities and the northeast region low velocities. This sharp velocity contrast appears to delineate the different abutting structures associated with Oita-Kumamoto Tectonic Line. However, because this region is not well resolved, more detailed structural studies are needed to confirm this observation.

Conclusion

Fig. 18 Horizontal and vertical Vp/Vs ratio perturbation sections for the study region. See Fig. 11 for map symbol explanations

type caldera can be formed by collapse along ring fractures after a large pyroclastic eruption (Williams 1941; Smith and Bailey 1968; Ono 1974; Komazawa 1995), with ensuing volcanic activity or geothermal hot springs peripheral to the collapse along narrow circular zones marginal to the gravity anomaly low.

In order to investigate the velocity heterogeneity associated with this active Aso Volcano, we carried out tomographic inversions for P- and S-velocity structure. An examination of the velocity perturbations for wellresolved regions clearly reveals velocity anomalies associated with the active crater at Aso Volcano. These results were obtained using 10298 (5453 P and 4845 S) travel times from 800 earthquakes and ten shots to determine estimates of 270 velocity parameters. Large-velocity anomalies are seen over the upper 11 km over the central and northern half beneath the central cones of Aso Caldera. The shallowest low-velocity region (LA) encompasses the upper 3 km of crust and is centered near the base of the collapsed caldera wall at Tateno Valley. The 10×15-km region covers the central north portion of the caldera and is characterized by P velocities up to 22% slower (27% for S). The second lowvelocity region (LD) is associated with the active central cones. The lowest velocities were centered at 6-km depth south of Mt. Kishima. Velocities as low as 4.3 km/s (up to 26%) in P and 2 km/s (31% lower) in S were found. The anomalous low-velocity volume covers a 7×6×7-km region and is roughly spherical in shape extending downward to 11 km, and extending to northeast. Its shape becomes more elliptical and flattens at depth by approximately 5 km. At 9-km depth, higher velocities are found centered in the southwestern part of the breached caldera wall in Tateno Valley (HF). The anomaly is characterized by velocities up to 11% faster in P, and up to 12% faster in S, and extends from the mouth of Tateno Valley. Small anomalies (5 km in width, 3 km in thickness) of faster-than-average velocities (up to 11% in S, 5% in P) are centered in the western portion of the caldera (HE) and slower-than-average velocities are found in the northern part of the caldera at the base of the wall (LC) and in Tateno Valley (LE).

343

We interpret the main low-velocity anomaly region (LD) imaged in this study to be associated with the active magma conduit system at Aso Volcano. The magma chamber that appears to be roughly spherical in shape, centered at 6 km depth, and is not located directly beneath Mt. Naka, the present location of magmatic activity, but rather between the central cones of Mt. Kishima, Mt. Eboshi, and Mt. Naka. The velocity contours suggest that the magma chamber may flatten at 10 km depth. The marked contrast in velocities, with the southwest and southeast regions showing high velocities (HC and HD) and the northeast region low velocities (LC and LD), may be delineating different abutting structures associated with Oita-Kumamoto Tectonic Line. Acknowledgements The data used in this study are part of the long-term, continuous volcano monitoring efforts of the staffs of Aso Volcanological Laboratory of Kyoto University. Thus, the authors are indebted to all staff members of AVL, especially H. Masuda, M. Sako, and T. Hoka for their assistance during the observations. Also thanks to A. Amato, A. Hurst, J.-F. Lénat, and an anonymous reviewer whose comments and suggestions helped to improve this paper. This study was funded under a NSF-Japan Postdoctoral Fellowship awarded to L.S.L.K., with further support provided by the USGS/Hawaiian Volcano Observatory.

References Arnott SK, Foulger GR (1994) The Krafla spreading segment, Iceland. 1. Three-dimensional crustal structure and the spatial and temporal distribution of local earthquakes. J Geophys Res 99:23801–23825 Aster RC, Meyer RP (1989) Determination of shear- and compressional-wave velocity variations and hypocenter locations in a rapidly inflating caldera: the Campi Flegrei. Phys Earth Planet Interior 55:313–325 Chatterjee SN, Pitt AM, Iyer HM (1985) Vp/Vs ratios in the Yellowstone National Park region, Wyoming. J Volcanol Geothermal Res 26:213–230 Chiarabba C, Amato A, Evans JR (1995) Variations on the NeHT high-resolution tomography method: a test of technique and results for Medicine Lake Volcano, northern California. J Geophys Res 100:4035–4052 Eberhart-Phillips D (1986) Three-dimensional velocity structure in northern California Coast Ranges for inversion of local earthquakes arrival times. Bull Seismol Soc Am 76:1025– 1052 Eberhart-Phillips D (1993) Local tomography: earthquake source regions. In: Iyer HM, Hirahara K (eds) Seismic tomography: theory and practice. Chapman and Hall, London, pp 613–643 Ehara S (1984) Seismic wave attenuation beneath geothermal areas in central Kyushu Japan. GRC Trans 8:489–492 Evans JR, Eberhart-Phillips D, Thurber CH (1994) User’s manual for SIMULPS12 for imaging Vp and Vp/Vs; a derivative of the “Thurber” tomographic inversion SIMUL3 for local earthquakes and explosions. US Geol Surv Open-File Rep 94: 1–101 Foulger GR, Miller AD, Julian BR, Evans JR (1995) Threedimensional vp and vp/vs structure of the Hengill Triple Junction and geothermal area Iceland and the repeatability of tomographic inversion. Geophys Res Lett 22:1309–1312 Francis TJG (1976) The ratio of compressional to shear velocity and rock porosity on the axis of the Mid-Atlantic Ridge. J Geophys Res 81:4361–4364 Kaneshima S, Kawakatsu H, Matsubayashi H, Sudo Y, Tsutsui T, Ohminato T, Ito H, Uhira K, Yamasato H, Oikawa J, Takeo M, Iidaka T (1996) Mechanism of phreatic eruptions at Aso Vol-

cano inferred from Near-Field Broadband Seismic Observations. Science 273:642–645 Kikuchi S, Ono H (1985) Observation of crustal deformation at Volcano Aso. Comparative study for geophysical background of volcanic activity and its relation to prediction of eruption disasters, pp 55–63 (in Japanese) Kissling E (1988) Geotomography with local earthquake data. Rev Geophys 26:659–698 Klein FW (1978) Hypocenter location program HYPOINVERSE 1 User’s guide to versions 1, 2, 3, 4. US Geol Surv Open-File Rep 78:1–103 Klein FW (1989) User’s guide to HYPOINVERSE a program for VAX computers to solve for earthquake locations and magnitudes. US Geol Surv Open-File Rep 89:1–58 Komazawa M (1995) Gravimetric analysis of Aso Volcano and its interpretation. J Geodetic Soc Japan 41:17–45 Kong LSL, Sudo Y (1993) Three-dimensional seismic velocity structure beneath Aso volcano Kyushu Japan. EOS Trans AGU 74:421 Kubotera A, Otsuka M (1971) Nature of microtremors observed on the thick deposit. Butsuri-Tanko (Geophys Explor) 24: 126–134 (in Japanese) Kubotera A, Mitsunami T (1980) An earthquake swarm in the northern part of the Aso caldera and migration of their foci. Tectonophysics 70:223–236 Kubotera A, Satomura M (1985) Repeated gravity measurements on the Aso Volcano. Comparative study for geophysical background of volcanic activity and its relation to prediction of eruption disasters, pp 64–68 (in Japanese) Kubotera A, Tajima H, Sumitomo N, Doi H, Izutuya S (1969) Gravity surveys on Aso Kuju volcanic region Kyushu district Japan. Bull Earthq Res Inst Univ Tokyo 47:215–255 Mavko GM (1980) Velocity and attenuation in partially molten rocks. J Geophys Res 85:5173–5189 New Energy Development Organization (NEDO) (1996) Geothermal development promotion survey in the western region of Aso Volcano, 1508 pp (in Japanese) Okubo Y, Shibuya A (1993) Thermal and crustal structure of the Aso Volcano and surrounding regions constrained by gravity and magnetic data, Japan. J Volcanol Geotherm Res 55:337– 350 Ono K (1974) On collapse-calderas. Association for the Geological Collaboration in Japan, Monogr 18, 4:55–61 (in Japanese) Ono K, Watanabe K (1985) Geological map of Aso Volcano. Geological map of volcanoes, no. 4. Geological Survey of Japan Pavlis GL, Booker JR (1980) The mixed discrete continuous inverse problem: application to the simultaneous determination of earthquake hypocenters and velocity structure. J Geophys Res 85:4801–4810 Sato H, Sacks IS, Murase T (1989) The use of laboratory velocity data for estimating temperature and partial melt fraction in the low-velocity zone: comparison with heat flow and electrical conductivity studies. J Geophys Res 94:5689–5704 Smith RL, Bailey RA (1968) Resurgent cauldrons. Mem Geol Soc Am 116:83–104 Stauber DA, Green SM, Iyer HM (1988) Three-dimensional P velocity structure of the crust below Newberry Volcano, Oregon. J Geophys Res 93:95–107 Sudo Y (1975a) Seismic activity in the northern region of Aso Volcano. Report about Aso Earthquake Swarm, January, 1975, pp 33–39 (in Japanese) Sudo Y (1975b) Seismometric study on earthquakes at the western region of the Aso Caldera: focal mechanism and tectonic line. Kazan Bull Volcanol Soc Japan 20:1–12 (in Japanese with English abstract) Sudo Y (1978) Activity of volcanic micro-tremors (1976–1977). Joint Geophysical and Geochemical Observations of Aso Volcano (August to December 1977), pp 1–3 (in Japanese) Sudo Y (1981) Seismic activities at the western region of the Aso caldera. Kazan Bull Volcanol Soc Japan 23:263–279 (in Japanese with English abstract)

344 Sudo Y (1984) Activity of volcanic micro-tremors (1978–1982). Joint Geophysical and Geochemical Observations of Aso Volcano (August to December 1981), pp 1–6 (in Japanese) Sudo Y (1988) Upper crustal structure of the Aso caldera. Kazan Bull Volcanol Soc Japan 33:130–134 (in Japanese) Sudo Y (1990) Characteristics of seismic activity around Aso Caldera. Chikyu Monthly 12:379–382 (in Japanese) Sudo Y (1991) An attenuating structure beneath the Aso caldera determined from the propagation of seismic waves. Bull Volcanol 55:99–111 Sudo Y (1993) Seismic activity and tectonic significance near the volcanic areas in central Kyushu, Japan. Mem Geol Soc Japan 41:19–34 (in Japanese with English abstract) Sudo Y, Yoshikawa S (1995) Gravity measurements at Aso Volcano. Chikyu Monthly 11:75–79 (in Japanese) Sudo Y, Yamada T, Masuda H (1984) Seismic activity and focal mechanism around the Aso caldera. Ann Dis Prev Res Kyoto Univ 27B-1:35–44 (in Japanese with English abstract) Thurber CH (1981) Earth structure and earthquake locations in the Coyote Lake area central California. PhD thesis, Mass Inst of Technol, Cambridge, Mass., 332 pp Thurber CH (1983) Earthquake locations and three-dimensional crustal structure in the Coyote Lake area, central California. J Geophys Res 88:8226–8236 Thurber CH (1987) Seismic structure and tectonics of Kilauea Volcano volcanism in Hawaii. US Geol Surv Prof Pap 1350: 919–934 Thurber CH (1993) Local earthquake tomography: velocities and vp/vs theory. In: Iyer HM, Hirahara K (eds) Seismic tomography: theory and practice. Chapman and Hall, London, pp 563–583 Toomey DR, Foulger GR (1989) Tomographic inversion of local earthquake data from the Hengill-Grensdalur central volcano complex Iceland. J Geophys Res 94:17497–17510

Tsutsui T, Sudo Y, Yoshikawa S, Ikawa T, Kuroda T (1997) A seismic exploration of Nango Valley Aso Caldera. Kazan Bull Volcanol Soc Japan 42:257–268 (in Japanese with English abstract) Watanabe K (1984) Active faults and their significance in the region to the west of Aso Caldera, Kumamoto Prefecture. Mem Faculty Educ Kumamoto Univ Nat Sci 33:35–47 (in Japanese with English abstract) Watanabe K, Momikura Y, Tsuruta K (1979) Active faults and parasitic eruption centers on the west flank of Aso caldera. Jpn Quaternary Res 18:89–101 (in Japanese with English abstract) Watanabe S, Kano N, Yamaguchi K, Yokokura T (1986) A subterranean structure of Aso-Dani by reflection seismic exploration Aso Caldera. Program Abstr Seis Soc Japan 1:73 (in Japanese) Wiggins RA (1972) The general linear inverse problem: implication of surface waves and free oscillations for earth structure. Rev Geophys Space Phys 10:251–285 Williams H (1941) Caldera and their origin. Bull Dept Geol Sci Univ Calif Publ 25:239–346 Wu H, Lees JM (1999) Three-dimensional P and S wave velocity structures of the Coso Geothermal Area, California, from microseismic travel time data. J Gephys Res 104:13217– 13233 Yamamoto M, Kawakatsu H, Kaneshima S, Mori T, Tsutsui T, Sudo Y, Morita Y (1999a) Detection of a crack-like conduit beneath the active crater at Aso volcano, Japan. Geophys Res Lett 26:3677–3680 Yamamoto M, Kawakatsu H, Kaneshima S, Iidaka T, Oikawa J, Watada S, Morita Y, Mori T, Tsutsui T, Sudo Y, Yoshikawa M, Hashimoto T, Nakabo M (1999b) ASOBOI97: Aso Seismic Observation with Broadband Instruments in 1997. Bull Earthq Res Inst 74:267–285 Yoshikawa K (1954) On the crustal movement of Volcano Aso (part 1). Zishin J Seismol Soc Japan 7:151–154 (in Japanese with English abstract)