## 1. Introduction

Subsurface structures reveal the current form of the static basement and crust and may also indicate the result of past crustal deformation. Consequently, subsurface structures are essentially fossil evidence of crustal movement, such as faulting, dyke emplacement, and/or volcanic activity, and the restoration of these subsurface structures provides important information on the tectonic and/or volcanic history of a region.

In volcanic areas, large eruptions produce widespread volcanic ejecta at the Earth’s surface, such as ash and lava flows, as well as volcanic depressions in the Earth’s crust. And, there are fossil magma chambers, that is, plutons. These are frequently recorded as gravity and magnetic anomalies, and especially collapse calderas formed by partial or total collapse of a magma chamber roof (e.g., [1–3]) have low gravity anomalies and it is a roughly concentric circular shape. Large eruption forming calderas occur with a very low frequency, however, since much smaller calderas such as the Miyake-jima caldera (e.g., [4]) form more frequently (e.g., [5]), understanding the mechanisms behind caldera formation is very important not only for advancing scientific knowledge but also for social purposes such as the construction of hazard maps.

Thus, studies on caldera formation have been conducted not only by geological surveys in the field but with analogue experiments, numerical simulations, and theoretical approaches (e.g., [6–21]). Although results obtained by these experiments have been compared with surface topography and ring-fault distributions, they have not been compared with subsurface structures or gravity anomalies reflecting subsurface structures. This is a result of the difficulty inherent in estimating caldera subsurface structures and in transforming analogue experiment results into gravity anomalies. Nevertheless, many researchers have recognized the importance of caldera subsurface structures and observed significant relationships between the dip of the caldera wall, the radius of the magma chamber, processes of caldera formation, and the type of caldera (e.g., [14, 19–21]).

Geoelectromagnetic and gravity surveys are popular geophysical techniques used in volcanic areas and have been employed frequently to estimate shallow to medium depth subsurface structures. In recent years, gravity gradiometry has been introduced. This measures the gravity gradient tensor generated by a source body, which consists of six components of three-dimensional (3D) gravity gradients. Gravity gradiometry survey has higher sensitivity than gravity surveys. Various analysis techniques for gravity gradient tensors have been developed and have given excellent results in subsurface structure estimation and edge detection, for example (e.g., [22–25]).

In this chapter, we first describe the gravity gradient tensor and its characteristics. Generally, the gravity gradient tensor is obtained by gravity gradiometry. However, surveys of this type have been made in only a few areas, so that tensor data are rarely available. Therefore, in Section 2.2, we present a method based on the work of Mickus and Hinojosa [26] that is used to obtain the tensor from the gravity anomaly. We then review the semiautomatic interpretation methods used to extract information on subsurface structures without additional geological and geophysical information, and apply some of these techniques to the volcanic zone of central Kyushu in Japan.

## 2. Gravity gradient tensor

### 2.1. Characteristics of gravity gradient tensor

The gravity gradient tensor Γ is defined by the differential coefficients of the gravitational potential *W* (e.g., [27, 28]), as follows:

Defining *g*_{x}, *g*_{y}, and *g*_{z}, as the first derivative of *W* along the *x, y*, and *z* directions, we can rewrite Eq. (1) as follows:

Here, *g*_{z} is the well-known gravity anomaly and *g*_{zz} is the vertical gradient of the gravity anomaly. The gravity gradient tensor is symmetric (e.g., [27]) and the sum of its diagonal components is zero since the gravitational potential satisfies Laplace’s equation (Δ*W* = 0).

Here, we set a point mass model under the surface (**Figure 1**). In this case, the gravitational potential *W* is:

where

*G* is the gravitational constant and *M* is the mass. Taking the first derivative along the *x*-, *y*-, and *z*-directions gives

and

These are *xyz* components of the gravitation vector, *F*, due to the mass anomaly *M*.

Taking the second derivatives of the potential *W* along *x, y*, and *z* directions, six of the nine second derivatives are given by

and

These are six independent components of the gravity gradient tensor. The other three second derivatives can be obtained using the symmetry characteristics of the second derivative. And, since the sum of Eqs. (8), (11), and (13) is zero, namely

diagonal components satisfy Laplace's equation (Δ*W* = 0).

We show a numerical example for the point mass model in **Figure 2**. In the numerical example, we set the point mass model with a depth of 4 km and a mass anomaly of –5.7 × 10^{4} kg. Here, negative mass anomaly means lack of mass caused by negative density contrast such as hole in the crust, and the mass anomaly of –5.7 × 10^{4} kg is equivalent to a sphere having radius of 3 km and density contrast of –500 kg/m^{3}. The vertical component, *g*_{z}, of the gravitational vector is negative because of a negative mass anomaly (**Figure 2(i)**). The horizontal components of the gravitational vector, *g*_{x} (**Figure 2(a)**) and *g*_{y} (**Figure 2(e)**), are positive in the area where the terms (*x* – *x*') or (*y* – *y*') in Eqs. (5) and (6) are positive, and negative in the area where the terms are negative, because material surrounding the point mass is denser and gravitation vectors point to dense area, that is, the outside of the point mass.

In the gravity gradient tensor, Γ, we first find the triplet pattern, negative-positive-negative, in the *g*_{xx} component (**Figure 2(b)**). This pattern can be understood by examining the slopes of the anomalies of *g*_{x} as we proceed from left to right along the *x*-axis. We notice that the steepest slope is at the center of the map of *g*_{x} and it is positive; the zero slopes are at the trough and peak of *g*_{x}, and the gentle negative slopes are to the left and right of the trough and peak, respectively. In a similar manner, we can explain the triplet pattern of the *g*_{yy} component (**Figure 2(g)**). In the *g*_{xy} component (**Figure 2(c)**), we find the quadruplet pattern, negative-positive-negative-positive. This pattern can be understood by examining the slopes of the anomalies of *g*_{y} as we proceed from left to right along the *x*-axis, or by examining the slopes of the anomalies of *g*_{x} as we proceed from left to right along the *y*-axis. The components *g*_{xz} (**Figure 2(d)**), *g*_{yz} (**Figure 2(h)**), and *g*_{zz} (**Figure 2(l)**) are differential coefficients of *g*_{x}, *g*_{y}, and *g*_{z} in the *z* direction and emphasize the high frequencies of each gravitational component without any changes in the location or shape of the anomalies.

### 2.2. Derivation of gravity gradient tensor from gravity anomaly

In general, the gravity gradient tensor is measured with gravity gradiometry (e.g., [29–31]). However, we can obtain the tensor from gravity anomaly data using calculations shown in [26], even if the gravity gradiometry surveys have not yet been undertaken.

As a technique for calculating the gravity gradient tensor from the gravity anomaly data, Mickus and Hinojosa [26] used the following procedure: (1) application of a Fourier transformation to the gravity anomaly, (2) estimation of the gravitational potential by integration of the gravity anomaly in the Fourier domain, (3) calculation of the gravity gradient components from second-order derivatives of the potential in each direction, and (4) application of a Fourier inverse transformation to finally obtain all components of the tensor in the spatial domain.

When the Fourier transformation of a function *f*(*x, y*) is denoted by the shorthand notation [*f*], that is,

The Fourier inverse transformation of the function *F*(*k*_{x}, *k*_{y}) in the Fourier domain is denoted by the shorthand notation ^{−1}[*F*] (e.g., [32]), that is,

where *k*_{x} and *k*_{y} are wave numbers and are inversely related to wavelengths *λ*_{x} and *λ*_{y} in the *x* and *y* directions, respectively: *k*_{x} = 2π/*λ*_{x} and *k*_{y} = 2π/*λ*_{y}. If *G*_{z} is the Fourier transform of the gravity anomaly *g*_{z}, the relationships between *g*_{z} and *G*_{z} are as follows:

and we use the notation *g*_{z} ↔ *G*_{z} for this relationship. Differentiation of a function in the space domain is equivalent to multiplication by a power of the wavenumber in the Fourier domain (e.g., [32]). If *f*(*x, y*) ↔ *F*(*k*_{x}, *k*_{y}), then,

Integration of the function in the Fourier domain is given by division of the wavenumber. From these and relationships between gravity and gravitational potential, we derive the following equations on gravitational vector components:

where *U* is the Fourier transform of the gravitational potential, *W*, and |**k**| = (*k*_{x}^{2} + *k*_{y}^{2})^{1/2}. From Eq. (21), we can obtain the gravitational potential *U* from the gravity anomaly by

and Eqs. (19) and (20) are rewritten as follows:

Therefore, we can obtain *g*_{x} and *g*_{y} from *g*_{z} by the Fourier inverse transformation of *G*_{x} and *G*_{y}.

In a similar manner, we derive equations to obtain each component of the gravity gradient tensor:

and

From these equations, as summarized in [26], we can obtain all components of the gravity gradient tensor from the gravity anomaly as follows:

## 3. Semiautomatic interpretation techniques

### 3.1. Edge emphasis

Edge emphasis techniques are extraction techniques used to find locations (namely, edge) where the gravity anomaly changes abruptly due to density variations, and these techniques play an important role in interpretation of potential field data (e.g., [33–35]).

In these techniques, the horizontal and vertical gravity gradient methods have frequently been employed to find structural boundaries such as faults or contacts between different materials (e.g., [36–39]). The vertical gravity gradient (*VG*) is the *g*_{zz} component of the gravity gradient tensor, and the horizontal gravity gradient (*HG*) is given using the components of the gravity gradient tensor shown in Eq. (2) as follows:

In addition to these gravity gradients, the second vertical derivative, ∂^{2}*g*_{z}/∂*z*^{2}, is well known and is a classical edge emphasis technique for edge detection (e.g., [40–43]).

Miller and Singh [44] proposed the *TDR* (tilt derivative), defined by the arctangent of the ratio of the vertical gravity gradient to the *HG*:

This emphasis technique does not extract extremely large amplitude signals in the short wavelength range and is known as a balanced method. However, it is noted that the *TDR* extracts a phantom boundary at the zero line if positive and negative anomalies exist side by side (e.g., [45, 46]).

Wijns et al. [47] proposed the *THETA* map that could be used to obtain the structural boundaries,

and Cooper and Cowan [44] proposed the *TDX* defined by the following equation:

Li et al. [48] conducted a numerical test of these emphasis techniques for structural boundary extraction. As a result, they found that (1) the *HG* method can determine structural boundaries of the causative body but its ability decreases with depth, (2) the *TDR* and *THETA* methods are also able to retrieve the structural boundaries but their extracted shapes are not clear and, as mentioned above, these methods introduce a phantom boundary at the zero line in case of side by side positive and negative anomalies, and (3) the *TDX* has the same phantom boundary limitation but the method can give clearer boundary location than both *TDR* and *THETA*.

Based on these discussions, Li et al. [48] suggested the *CLP* method, which extracts the outlines of subsurface structures without the phantom boundary. The *CLP* is defined as follows:

where *p* is a previously determined constant, typically 1/10 of the maximum value of the horizontal gravity gradient. And *k* is defined as follows:

In recent years, specialized edge emphasis techniques have been developed to find the edges of potential field data related to multi subsurface structures. For example, Ma [46] suggested the Improved Local Phase method (*ILP*):

And, Ferreira et al. [49] have suggested the *TAHG*:

Using numerical tests, Zhang et al. [50] showed that the *ILP* and the *TAHG* successfully extract the edges of potential field data relating to multi subsurface structures, and they also suggested a more sensitive edge detection method, *THVH*. The *THVH* is defined as follows:

Kusumoto [51] pointed out that because the *CLP, ILP, TAHG*, and *THVH* are all very sensitive high pass filters, which respond to very small signals (perhaps also including noise), other geological or geophysical data of the region would be required in order to fully understand the results obtained by these techniques.

### 3.2. Curvature of the potential field

Curvatures of potential field data vary in response to density changes in the subsurface structure, and these curvatures are described by the *g*_{xx}, *g*_{yy}, and *g*_{xy} components (e.g., [27, 28]).

Pedersen and Rasmussen [52] defined the invariant ratio, *I*, of the gravity gradient tensor as follows:

Here, *I*_{1} and *I*_{2} are invariants of the tensor. Each invariant is given by three eigenvalues (*λ*_{1}, *λ*_{2}, *λ*_{3}) of the gravity gradient tensors, as follows:

The invariant ratio *I* varies between 0 and +1. If the body causing a gravity anomaly is a 2D structure, *I* is 0. If the causative body is 3D, *I* is +1. Beiki and Pedersen [53] named this ratio the dimensionality index and suggested *I* = 0.5 as the threshold between a 2D and 3D causative structure. In addition, for the case of a 2D causative structure such as a dike, the maximum eigenvector of the gravity gradient tensor points to the causative body and the minimum eigenvector is parallel to the strike direction of the structure.

Cevallos et al. [54] and Cevallos [55] pointed out that the shape index (e.g., [56, 57]) is useful for determining the characteristics of potential fields and is defined by

The shape index, *Si*, varies from –1 to +1, and it is suggested that values of –1, –0.5, 0, +0.5, and +1 correspond to bowl, valley, flat, ridge, and dome structures, respectively. Cevallos [55] found that a map of the shape index at depth was consistent with deep structures integrated by geophysical and geological data in the King Sound area of the Canning Basin, Western Australia.

Studies on the curvature of the potential field have recently been developed as new interpretation techniques because they use a lot of the characteristics of the gravity gradient tensor (e.g., [58–62]). However, some problems relating to its practical usage remain unsolved. In fact, Li [63] conducted numerical tests for 13 well-known semiautomatic interpretation methods based on the curvature of the potential field and concluded that the shape of the potential does not always correspond to the subsurface structure, and therefore care should be taken when interpreting these results.

### 3.3. Euler deconvolution

Euler deconvolution is a semiautomatic interpretation method based on Euler’s homogeneity equation and is often employed to estimate locations and/or outlines of causative bodies. Because the Euler deconvolution technique can provide rapid interpretations of any potential field data in terms of depth and geological structures, it has been used by several researchers for analyzing both magnetic anomalies (e.g., [64–66]) and gravity anomalies (e.g., [67]).

The Euler deconvolution based on the three orthogonal gradient components of the potential field is simply called Euler deconvolution or conventional Euler deconvolution. For the gravity anomaly *g*_{z}, it uses the following equation:

In Eq. (44), *x*_{0},* y*_{0}, and *z*_{0} are unknown location parameters of the causative body center or edge that have to be estimated, and *x, y*, and *z* are known location parameters of the observation point of the gravity and its gradients. *N* is the structural index and *B*_{z} is the regional component of the gravity anomaly that has to be estimated.

Rewriting Eq. (44), we obtain the following equation in which we can separate unknown parameters from known parameters:

There are four unknown parameters in Eq. (45). If we have enough data, namely *n* ≥ 4, to solve for these unknowns within a selected window, we can estimate these four parameters using the least-squares method.

Zhang et al. [68] suggested the Tensor Euler deconvolution, designed to consider the full gravity gradient tensor (Eq. (2)) and all components (*g*_{x}, *g*_{y}, *g*_{z}) of the gravitational vector. The Tensor Euler deconvolution is defined by the following simultaneous equations:

(46) |

Here, *B*_{x} and *B*_{y} are the regional components of the gravity anomaly that have to be estimated. Thus, there are six unknown parameters in this equation, and if we have enough data to solve for these six unknowns within a selected window, namely *n* ≥ 6, these parameters can be estimated using the least squares method. Using a numerical study and applying it to field data, Zhang et al. [68] showed that the Tensor Euler deconvolution produces a significantly improved resolution compared to the conventional Euler deconvolution.

### 3.4. Dip estimation of a fault or structural boundary

In recent years, a technique estimating the dip of structures such as dikes, faults, or other geological structure boundaries has been developed and produced good results as a semiautomatic interpretation method using the gravity gradient tensor.

Beiki [69] suggested that the dip, *α*, of the causative body could be estimated from the three components (*x, y, z*) of the maximum eigenvector (*v*_{1}) of the gradient tensor of the potential field (**Figure 3**) as follows:

where *v*_{1x}, *v*_{1y}, and *v*_{1z} are the *x, y*, and *z* components of the maximum eigenvector, *v*_{1}. Beiki [69] applied this method to an aeromagnetic data set of the Åsele area in Sweden, and provided useful information on the dip of dike swarms. Kusumoto [70] applied this technique to estimate the dip of the Median Tectonic Line located in southwest Japan, and obtained a dip distribution compatible with results from other geophysical surveys, including a reflection survey. Since gravity gradient tensor consists of differential coefficients of gravitation caused by causative body the tensor has higher accuracy than gravity anomaly. In addition, since analysis using eigenvalues and eigenvectors of the tensor considers all components of the tensor, this technique would give good results. In this technique, we do not need to know or assume density contrast between layers or structures.

## 4. Application of semiautomatic interpretation method

Here, we apply the semiautomatic interpretation methods shown in the previous section to a volcanic area located in central Kyushu, Japan. This area consists of the Aso caldera and the Hohi Volcanic Zone containing a buried caldera, the Shishimuta caldera, where many previous geological and geophysical studies have been conducted. In addition, a database of the Bouguer gravity anomaly is available for use (e.g., [71]), although a gravity gradiometory survey covering this area has not yet been carried out. Since this gravity anomaly database employed in this study gives users 1 km × 1 km mesh data, we discussed structures larger than several kilometers in this study. If we would like to discuss fine structures, it would be possible by employing denser mesh data.

### 4.1. Background of Hohi Volcanic Zone and Aso caldera

#### 4.1.1. Hohi Volcanic Zone and Shishimuta caldera

The Hohi Volcanic Zone (e.g., [73]) is located in the eastern part of central Kyushu (**Figure 4**) and is characterized by a wedge-shaped low gravity anomaly area, which becomes narrow toward the east (**Figure 5**). In the Hohi Volcanic Zone, there are four conspicuous low anomaly areas toward the west-southwest from Beppu Bay, which correspond to the Beppu Bay, Shonai Basin, Shishimuta caldera, and Kuju basin, respectively. The Beppu Bay has the lowest gravity anomaly, which reaches –50 mGal.

Subsurface structures in this region have been estimated by gravity anomaly data, drill core data, reflection survey, and surface geology data (e.g., [74, 75]). Kusumoto et al. [75] showed that, although the average depth of the basement was estimated to be approximately 1.2 km in the Hohi Volcanic Zone, the small basins (the Shonai basin and the Kuju basin) along the Oita-Kumamoto Tectonic Line (the tectonic line suggested by a steep gravity gradient (e.g., [76])) are deeper than 2 km. In addition, it was shown that the basement depth in Beppu Bay reaches to approximately 4 km. Detailed geological surveys have shown that this volcanic zone consists of half-grabens detailing volcanic activity that began 6 Ma (e.g., [73]) and pull-apart basins that formed after the half-graben formation (e.g., [77]).

During formation of the Hohi Volcanic Zone, it is known that eruption of more than 5000 km^{3} of material occurred in the last six million years, and that eruption rate gradually decreased. The volcanic zone erupted 2900 km^{3}/Ma in its initial stage of activity from 5 to 4 Ma. After that, the zone erupted approximately 1300, 900, 400, and 200 km^{3} between 4–3, 3–2, 2–1, and 1–0 Ma, respectively [73]. The formation of the pull-apart basin began 1.5 Ma; this tectonic change occurred due to a counter-clockwise shift of the subducting Philippine Sea Plate (e.g., [77]).

Kusumoto et al. [78] attempted to restore the tectonic structures by numerical simulations based on the dislocation theory (e.g., [79]). They showed that tectonic structures such as half-grabens and pull-apart basins can be restored by obeying the tectonics suggested by Itoh et al. [77]. Kusumoto et al. [80] assumed a simple two-layer structure model, which consisted of the basement and the sediment, and estimated the gravity anomaly field from the vertical displacement field of the basement obtained by tectonic modeling [78]. As a result, it was shown that a gravity anomaly of volcanic origin cannot be restored. In this model, the unrestored gravity anomaly was that relating to the Shishimuta caldera.

The Shishimuta caldera is a buried caldera and is the origin of the Yabakei and Imaichi pyroclastic flows that are widely distributed throughout central Kyushu (e.g., [81, 82]). The volumes of these flows have been estimated as 40 km^{3} (e.g., [81]) and 90 km^{3} (e.g., [82]), respectively, based on detailed surface geological surveys. The red circle shown in **Figure 4** indicates the caldera wall at 1 km depth with a diameter of approximately 8–10 km, estimated from drilling core data and gravity anomalies (e.g., [81]). Density contrast between basement and sedimentary layers is estimated to be in the range of 300–700 kg/m^{3} by Kusumoto et al. [75]. As shown in **Figure 5**, this Shishimuta caldera does not have the concentrically circular low gravity anomaly found in many calderas; rather, this caldera has a triangular low gravity anomaly.

The subsurface structure of the Shishimuta caldera was estimated from a geological perspective, and the structure shallower than 3 km depth was revealed by drill core data and gravity anomalies. On the basis of drill core data, gravity anomalies, active fault distribution, microearthquake activity, and geological surveys in and near the Shishimuta caldera, Kamata [81] estimated that the shape of this caldera is a funnel type. In addition, Kamata [81] estimated a magma chamber depth of the Shishimuta caldera would be 7–12 km depth. This was estimated by assuming that the original depth of the Yabakei pyroclastic flow is equal to the source depth of the Ito pyroclastic flow of Aira caldera, Kagoshima (e.g., [83]).

#### 4.1.2. Aso caldera

The Aso caldera, like the Shishimuta caldera, is located in central Kyushu (**Figure 4**). However, it extends approximately 25 km in the south-north direction and about 18 km in the east-west direction and is thus elliptical in plan view, and is one of the largest calderas in the world.

It is known that the Aso caldera was formed by four eruptions with large-scale pyroclastic flows. The first large pyroclastic flow is called Aso-1 and flowed out from the present Aso caldera prior to 0.27 Ma. The second is called Aso-2 and was erupted before 0.14 Ma. The Aso-2 was the smallest pyroclastic flow and its volume is estimated at approximately 25 km^{3}. The third pyroclastic flow is called Aso-3 and is dated at 0.12 Ma. The fourth pyroclastic flow is called Aso-4 and is the largest pyroclastic flow with a volume of more than 80 km^{3}. From detailed geological surveys, it has been found that Aso-4 would have flown over the sea (the Seto inland sea and Tachibana Bay) and reached Yamaguchi Prefecture and Shimabara Peninsula. In the Aso caldera, there is also a central group of volcanic cones, which are currently active (e.g., [84]).

Gravity surveys of this region have been initiated by Tsuboi et al. [76] and Kubotera et al. [85]. Yokoyama [86] estimated that the subterranean caldera structure is similar to a funnel-shape. Komazawa [87] accumulated gravity data and conducted high accuracy, high-resolution analyses for the compiled gravity data. As a result, they found that the Aso caldera has a piston-cylinder type structure rather than a funnel-shaped structure with a single low anomaly. In addition, five minor local gravity lows exist in the caldera, which make up the major low gravity zone of the Aso caldera. The major gravity low has a steep gradient inside the caldera rim, and the central area of the minor gravity anomalies has a relatively flat bottom.

### 4.2. Gravity gradient tensor in central Kyushu

In **Figure 6**, we show all components of the gravitation vector and gravity gradient tensor. Since the gravity gradiometry survey covering this region has not been conducted, the gravitation vector and gravity gradient tensor were calculated from the Bouguer gravity anomaly data shown in **Figure 5**.

Although all components of the gravity gradient tensor were theoretically estimated by the method of Mickus and Hinojosa [26], we also obtained components of the gravity gradient tensor by numerical differentiation of *g*_{x}, *g*_{y}, and *g*_{z} in the space domain, since short wavelength noise due to differentiation in the Fourier domain can sometimes be generated. The gravitation vector was, however, calculated using the method of Mickus and Hinojosa [26].

## 5. Results and discussion

### 5.1. Edge emphasis

In this section, we applied the 10 edge emphasis techniques shown in Section 3.1 to the field data. The results obtained by each method are shown in **Figures 7**–**16**.

The maps of the horizontal gravity gradient (*HG*: Eq. (32)), the vertical gravity gradient (*VG: g*_{zz}), and the second vertical derivative (∂^{2}*g*_{z}/∂*z*^{2}, or *g*_{zzz}) are shown in **Figures 7**–**9**, respectively. The horizontal gravity gradient map extracts the outlines of the Hohi Volcanic Zone, the Shishimuta caldera, and the Aso caldera. In addition to these caldera or volcanic depressions, tectonic lines such as the Median Tectonic Line and Oita-Kumamoto Tectonic Line, which are important structures in this region, are also extracted (**Figure 7**). These lineaments were extracted as a high gravity gradient belt with an NE-SW trend crossing the Aso caldera. The land area of the belt corresponds to the Oita-Kumamoto Tectonic Line, and the sea area of the belt corresponds to the Median Tectonic Line. The Median Tectonic Line is the largest tectonic line in the southwest Japan, and it is known that the tectonic line has moved as a right lateral fault with normal faulting in the Quaternary. Furthermore, we also detect a large circular belt with a high gravity gradient outside of the high gradient belt that represents the Shishimuta caldera.

The vertical gravity gradient extracts the boundaries of the low gravity and high gravity areas, and emphasizes areas of low density (**Figure 8**). In particular, the caldera rim of the Aso caldera and boundaries of the tectonic lines and sedimentary basins are clearly extracted. In addition, areas of low-density material such as ash, pyroclastic flows, and sediments are accurately predicted and emphasized, and correspond to the inner areas of the Aso caldera, Shishimuta caldera, Kuju basin, Shonai basin, and Beppu Bay. The *g*_{zzz} map (**Figure 9**) shows each rim of the Aso caldera and Shishimuta caldera more clearly and emphasizes shorter wavelength signals. However, it seems that this technique does not highlight the density structures as effectively as the vertical gravity gradient, *g*_{zz}.

The maps of the *TDR* (Eq. (33)), the *THETA* (Eq. (34)), the *TDX* (Eq. (35)), and the *CLP* (Eq. (36)) are shown in **Figures 10**–**13**, respectively. The *TDR* shows the original structural boundaries at zero lines [44]. If we obeyed this judgment criterion, the boundaries of the green and yellow areas would be the structural boundaries such as lineaments, faults, and caldera rims. And, if we did not persist in the criterion, from this *TDR* we can derive information on both the structural boundaries and the density structures extracted by *HG* and *VG*. In fact, the *TDR* shows local low gravity areas more clearly or in greater detail than the vertical gravity gradient map.

The *THETA* map, the *TDX* map, and the *CLP* map extract structural boundaries by lines. In particular, it appears that the *CLP* technique is an excellent boundary extraction technique as it extracts very distinctly the caldera rims of the Aso caldera and the Shishimuta caldera, as well as locations of the Median Tectonic Line and the Oita-Kumamoto Tectonic Line (**Figure 13**). In addition, we also observe the larger circular lines outside the yellow lines indicating the Shishimuta caldera rim and the Aso caldera rim. The larger circular line outside of the Shishimuta caldera was also shown in the horizontal gravity gradient (*HG*) map (**Figure 7**).

The maps of the *ILP* (Eq. (38)), *TAHG* (Eq. (39)), and *THVH* (Eq. (40)) are shown in **Figures 14**–**16**, respectively. These techniques have been developed for finding the edges of the potential field related to multi subsurface structures, and generally they emphasize very small signals caused by the deep structures. The *ILP* map extracts the Aso caldera rim, but it is difficult to interpret what the other strong short wavelength lines represent. The *TAHG* and *THVH* maps are similar to the *ILP* map and have similar interpretation difficulties.

### 5.2. Curvature analysis

**Figures 17** and **18** are the dimensionality index map and the shape index map. The dimensionality index map indicates whether the subsurface structure in the region is 2D, such as a vertical or subvertical dike or fault, or 3D, such as a dome structure. If the subsurface structure is 2D, the dimensionality index, *I* < 0.5. Since areas of *I* ≤ 0.5 are widely distributed in **Figure 17**, we can conclude that the subsurface structures in the study area predominantly consist of 2D structures. However, in the inner area of the Shishimuta caldera rim, the dimensionality index is more than 0.7, and 3D structures are therefore predicted (**Figure 17**). In addition, the dimensionality index is also more than 0.7 in part of the inner area of the Aso caldera rim. This area corresponds to the volcanic cone group in the Aso caldera.

The shape index in the Shishimuta caldera rim reaches –1 (**Figure 18**). This indicates that the potential shape of the low-density material is either bowl-like, a downward convex structure, or both. In fact, since Kamata [81] predicted that the Shishimuta caldera would be a funnel type caldera, the obtained results are considered suitable. Furthermore, these results are compatible with the result shown in **Figure 17**; the inner structure of the Shishimuta caldera would be 3D. In addition, locally low shape index areas are distributed around the Shishimuta caldera.

In the Aso caldera, the shape index is in the range of –1 to +1. This indicates that the subsurface structures would be complex since the Aso caldera has experienced four caldera formations. Moreover, areas of locally low shape index less than –0.75 correspond to the five local minor gravity lows found by Komazawa [87].

### 5.3. Euler deconvolution

We employed here the Tensor Euler deconvolution [68] to obtain information on the subsurface structures. In the calculation, we set a 24 km width window (the required amount of data to solve Eq. (9) is 576) with a structure index of 0.001. Generally, a structure index of 0 is employed for sills, dikes, and faults (e.g., [88]), but calculations with this value did not give solutions in our study field. The value of 0.001 was obtained through trail and error.

In **Figure 19**, we show the solution clouds given by the Tensor Euler deconvolution. We identify a deep and flat structure like a sedimentary basin at the east of the Aso caldera. There are large structural gaps between this structure and the Hohi Volcanic Zone, and the gaps appear to resemble a cliff (**Figure 19b**). In the Aso caldera, structural boundaries within the caldera are extracted as a void area of the solutions. However, these inner boundaries were not extracted in the Shishimuta caldera. In addition, we detect a large circular wall structure surrounding the Shishimuta caldera. Although partially observable in map view (**Figure 19a**), the structure is more clearly distinguished in bird’s-eye view (**Figure 19b**).

### 5.4. Dip estimation

Dips of faults and/or structural boundaries were estimated in the areas that satisfied the conditions *I* (dimensionality index) < 0.5 and *HG* (horizontal gravity gradient) ≥ 20 E (**Figure 20**). From **Figure 20**, we found that dip of the Median Tectonic Line and Oita-Kumamoto Tectonic Line are in the range of 40–70°, and are interpreted as high angle faults. In part of the Median Tectonic Line, the dip is very high, exceeding 70°. Since it is known that the Median Tectonic Line moved as a right lateral fault with normal faulting in the Quaternary, it is reasonable that this tectonic line has a high dip reaching 70° or 80°. The dip of the Median Tectonic Line becomes gradually lower to the north, toward Beppu Bay area. Seismic reflection surveys (e.g., [89]) confirm that dip of normal faults distributed in the Beppu Bay are <45°. Consequently, the obtained dip distribution is consistent with the observed data, suggesting that this technique is a useful semiautomatic interpretation method. In our results, the dip of the Oita-Kumamoto Tectonic Line is estimated to be approximately 60°, but its actual dip is unknown. Since, for many numerical simulations evaluating deformation fields by fault motions, the dip of the fault plays an important role, this technique will give us useful information on fault shape even if the result is a first approximation solution.

In the Aso caldera, dip of its caldera wall are in the range of 50–70° (inward dipping), as shown in **Figure 20**. Dips become gradually lower toward the caldera center, which is almost flat. These subsurface structural characteristics in the inner area of the Aso caldera have already been observed by detailed gravity analysis [87], and our results and Komazawa's results support to each other. However, this is the first time that caldera wall dips have been explicitly determined. In the Shishimuta caldera, because the condition of the dimensionality index was not satisfied, the distribution of caldera wall dips could not be estimated. Since the dip of the caldera wall leads to a discussion of the radius of the magma chamber, the process of caldera formation, and the type of caldera (e.g., [9, 14, 19–20]), the dip of the wall is thought to be as important as the fault dip.

The dips of the large circular wall structure surrounding the Shishimuta caldera are in the range of approximately 45–60° (inward dipping). This large structure appears in *HG* (**Figure 7**), *VG* (**Figure 8**), *TDR* (**Figure 10**), *CLP* (**Figure 13**), shape index (**Figure 18**), and Tensor Euler deconvolution (**Figure 19**) data, and it has a diameter of approximately 35 km. Since the Hohi Volcanic Zone was formed by eruption of more than 5000 km^{3} or material during six million years (e.g., [73]), it is suggested that this structure represents a depression or group of depressions resulting from these eruptions. However, because the obvious sedimentary basins such as the Kuju basin and the Shonai basin are included within this area, more detailed surveys and studies would be required to discuss this large structure surrounding the Shishimuta caldera.

## 6. Conclusions

By obtaining the gravity gradient tensor from the gravity anomaly, we can use a variety of techniques to extract subsurface structures. In this study, we reviewed semiautomatic interpretation methods using the gravity gradient tensor, and applied some of the techniques to the volcanic zone of central Kyushu, Japan. This area consists of the Aso caldera and the Hohi Volcanic Zone containing a buried caldera, the Shishimuta caldera, and has large tectonic lines such as the Median Tectonic Line and the Oita-Kumamoto Tectonic Line, and tectonic sedimentary basins such as the Beppu Bay, the Kuju basin, and the Shonai basin.

Most edge detection methods extracted the outlines of these structures, and some of them indicated density structures. In spite of classical methods, the horizontal gravity gradient method and the vertical gravity gradient method were excellent edge detection methods. As for recent methods, the *CLP* method clearly extracted structural boundaries such as the caldera rim and tectonic lines as lines.

In the curvature analysis, we obtained useful information from the shape index and the dimensionality index that indicated caldera shape. In the estimation of the dip of structural boundaries using the eigenvector of the gravity gradient tensor, we obtained the fault dip of the Median Tectonic Line, which is consistent with seismic reflection surveys, and estimated the caldera wall dip of the Aso caldera, which corresponds to that obtained by detailed gravity analysis.

Through these data analyses, we distinguished a large circular structure with a diameter of 35 km surrounding the Shishimuta caldera. This structure appeared also in solution clouds obtained by the Tensor Euler deconvolution. However, we cannot confidently judge what this structure represents using only the gravity gradient tensor. More detailed surveys and studies are required to further discuss this issue.