Structural Analysis of Calderas by Semiautomatic Interpretation of the Gravity Gradient Tensor: A Case Study in Central Kyushu, Japan

Abstract

In order to understand the formation history and mechanism of volcanoes and their related structures, including calderas, their subsurface structures play an important role. In recent years, gravity gradiometry survey has been introduced, and new analyses techniques for gravity gradient tensors obtained by the survey have been developed. In this chapter, we first describe the gravity gradient tensor and its characteristics, and the method for obtaining the tensor from the gravity anomaly. Next, we review the semiautomatic interpretation methods for extracting information on subsurface structures, and apply some of the techniques to the volcanic zone of central Kyushu, Japan. The results showed that the horizontal and vertical gravity gradient methods, and the CLP method were useful for extracting outlines of important volcanic and tectonic structures in this region. Using the maximum eigenvector of the gravity gradient tensor, the caldera wall dip of the Aso caldera was successfully estimated to be in the range of 50–70°, and the dip of the Median Tectonic Line which was the largest tectonic line in the southwest Japan was consistent with seismic reflection surveys. In addition, a large circular structure surrounding the Shishimuta caldera with a diameter of 35 km was distinguished in some analyses.

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:

$\Gamma =\left[\begin{array}{ccc}\frac{{\partial }^{2}W}{\partial x^2}& \frac{{\partial }^{2}W}{\partial x\partial y}& \frac{{\partial }^{2}W}{\partial x\partial z}\\ \frac{{\partial }^{2}W}{\partial y\partial x}& \frac{{\partial }^{2}W}{\partial y^2}& \frac{{\partial }^{2}W}{\partial y\partial z}\\ \frac{{\partial }^{2}W}{\partial z\partial x}& \frac{{\partial }^{2}W}{\partial z\partial y}& \frac{{\partial }^{2}W}{\partial z^2}\end{array}\right]$E1

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:

$\Gamma =\left[\begin{array}{ccc}\frac{\partial g_x}{\partial x}& \frac{\partial g_x}{\partial y}& \frac{\partial g_x}{\partial z}\\ \frac{\partial g_y}{\partial x}& \frac{\partial g_y}{\partial y}& \frac{\partial g_y}{\partial z}\\ \frac{\partial g_z}{\partial x}& \frac{\partial g_z}{\partial y}& \frac{\partial g_z}{\partial z}\end{array}\right]=\left[\begin{array}{ccc}{g}_{xx}& g_{xy}& g_{xz}\\ g_{yx}& g_{yy}& g_{yz}\\ g_{zx}& g_{zy}& g_{zz}\end{array}\right]$E2

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:

$W=-\frac{GM}{r}$E3

where

$r=\sqrt{{\left(x-x^{\prime }\right)}^2+{\left(y-y^{\prime }\right)}^2+{\left(z-z^{\prime }\right)}^2}$E4

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

$g_x=W_x=\frac{\partial W}{\partial x}=GM\frac{x-x^{\prime }}{r^3}$E5

$g_y=W_y=\frac{\partial W}{\partial y}=GM\frac{y-y^{\prime }}{r^3}$E6

and

$g_z=W_z=\frac{\partial W}{\partial z}=GM\frac{z-z^{\prime }}{r^3}$E7

Figure 1.

Point mass model. Schematic illustration of the gravitational vector due to a positive mass anomaly.

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

$g_{xx}=W_{xx}=\frac{{\partial }^{2}W}{\partial x^2}=GM\frac{-3{\left(x-x^{\prime }\right)}^2+r^2}{r^5}$E8

$g_{xy}=W_{xy}=\frac{{\partial }^{2}W}{\partial x\partial y}=GM\frac{-3\left(x-x^{\prime }\right)\left(y-y^{\prime }\right)}{r^5}$E9

$g_{xz}=W_{xz}=\frac{{\partial }^{2}W}{\partial x\partial z}=GM\frac{-3\left(x-x^{\prime }\right)\left(z-z^{\prime }\right)}{r^5}$E10

$g_{yy}=W_{yy}=\frac{{\partial }^{2}W}{\partial y^2}=GM\frac{-3{\left(y-y^{\prime }\right)}^2+r^2}{r^5}$E11

$g_{yz}=W_{yz}=\frac{{\partial }^{2}W}{\partial y\partial z}=GM\frac{-3\left(y-y^{\prime }\right)\left(z-z^{\prime }\right)}{r^5}$E12

and

$g_{zz}=W_{zz}=\frac{{\partial }^{2}W}{\partial z^2}=GM\frac{-3{\left(z-z^{\prime }\right)}^2+r^2}{r^5}$E13

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

$\begin{array}{lll}{g}_{xx}+g_{yy}+g_{zz}\hfill & =\hfill & GM\left[\frac{-3{\left(x-x^{\prime }\right)}^2+r^2}{r^5}+\frac{-3{\left(y-y\text{'}\right)}^2+r^2}{r^5}+\frac{-3{\left(z-z\text{'}\right)}^2+r^2}{r^5}\right]\hfill \\ \hfill & =\hfill & GM\left\{\frac{-3\left[{\left(x-x^{\prime }\right)}^2+{\left(y-y\text{'}\right)}^2+{\left(z-z\text{'}\right)}^2\right]+3r^2}{r^5}\right\}=0\hfill \end{array}$E140

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⁴ 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⁴ kg is equivalent to a sphere having radius of 3 km and density contrast of –500 kg/m³. 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.

Figure 2.

Components of the gravitational vector and gravity gradient tensor caused by a point mass with a depth of 4 km and a mass anomaly of –5.7 × 10⁴ kg. Negative mass anomaly means lack of mass caused by negative density contrast such as hole in the crust. Negative areas are colored gray. The gravitational vector unit is given in mGal, and mGal = 10^–5 m/s². The unit of the gravity gradient tensor is given in E (Eötvös), and 1 E = 0.1 mGal/km. (a) g_x component of gravitational vector, (b) g_xx component of gravity gradient tensor, (c) g_xy component, (d) g_xz component, (e) g_y component, (f) g_yx component, which is equal to g_xy, (g) g_yy component, (h) g_yz component, (i) g_z component; gravity anomaly, (j) g_zx component, which is equal to g_xz, (k) g_zy component, which is equal to g_yz, and (l) g_zz component.

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,

$F\left[f\right]=F\left(k_x,k_y\right)={\displaystyle {\int }_{-\infty }^{\infty }{\displaystyle {\int }_{-\infty }^{\infty }f\left(x,y\right)e^{-i\left(k_{x}x+k_{y}y\right)}dxdy}}$E14

The Fourier inverse transformation of the function F(k_x, k_y) in the Fourier domain is denoted by the shorthand notation ⁻¹[F] (e.g., [32]), that is,

$F^{-1}\left[F\right]=f\left(x,y\right)=\frac{1}{4{\pi }^2}{\displaystyle {\int }_{-\infty }^{\infty }{\displaystyle {\int }_{-\infty }^{\infty }F\left(k_x,k_y\right)e^{i\left(k_{x}x+k_{y}y\right)}d{k}_{x}d{k}_y}}$E15

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:

$G_z=F\left[g_z\right]$E16

$g_z=F^{-1}\left[G_z\right]$E17

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,

$\frac{{\partial }^m}{\partial x^m}\frac{{\partial }^n}{\partial y^n}f\left(x,y\right)\leftrightarrow {\left(i{k}_x\right)}^m{\left(i{k}_y\right)}^{n}F\left(k_x,k_y\right)$E18

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:

$G_x=i{k}_{x}U$E19

$G_y=i{k}_{y}U$E20

$G_z=\left|k\right|U$E21

where U is the Fourier transform of the gravitational potential, W, and |k| = (k_x² + k_y²)^1/2. From Eq. (21), we can obtain the gravitational potential U from the gravity anomaly by

$U=\frac{G_z}{\left|k\right|}$E22

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

$G_x=\frac{i{k}_x}{\left|k\right|}{G}_z$E23

$G_y=\frac{i{k}_y}{\left|k\right|}{G}_z$E24

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:

$G_{xx}=-k_x^{2}U=\frac{-k_x^2}{\left|k\right|}{G}_z$E25

$G_{yy}=-k_y^{2}U=\frac{-k_y^2}{\left|k\right|}{G}_z$E26

$G_{zz}=\left|k\right|G_z$E27

$G_{xy}=i{k}_{x}i{k}_{y}U=\frac{-k_x{k}_y}{\left|k\right|}{G}_z$E28

$G_{zx}=i{k}_x{G}_z$E29

and

$G_{zy}=i{k}_y{G}_z$E30

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

$\left[\begin{array}{ccc}{g}_{xx}& g_{xy}& g_{xz}\\ g_{yx}& g_{yy}& g_{yz}\\ g_{zx}& g_{zy}& g_{zz}\end{array}\right]=F^{-1}\left\{\left[\begin{array}{ccc}\frac{-k_x^2}{\left|k\right|}& \frac{-k_x{k}_y}{\left|k\right|}& i{k}_x\\ \frac{-k_x{k}_y}{\left|k\right|}& \frac{-k_y^2}{\left|k\right|}& i{k}_y\\ i{k}_x& i{k}_y& \left|k\right|\end{array}\right]G_z\left(k\right)\right\}$E31

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:

$HG=\sqrt{g_{zx}^2+g_{zy}^2}$E32

In addition to these gravity gradients, the second vertical derivative, ∂²g_z/∂z², 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:

$TDR=\arctan \left(\frac{g_{zz}}{\sqrt{g_{zx}^2+g_{zy}^2}}\right)$E33

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,

$THETA=\arccos \left(\frac{\sqrt{g_{zx}^2+g_{zy}^2}}{\sqrt{g_{zx}^2+g_{zy}^2+g_{zz}^2}}\right)$E34

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

$TDX=\arctan \left(\frac{\sqrt{g_{zx}^2+g_{zy}^2}}{\left|g_{zz}\right|}\right)$E35

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:

$CLP=\arctan \left[\frac{\sqrt{{\left(\frac{\partial g_z}{\partial x}\right)}^2+{\left(\frac{\partial g_z}{\partial y}\right)}^2}}{p+k\times \left|\frac{\partial g_z}{\partial z}\right|}\right]$E36

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:

$k=\frac{mean\left|\frac{\partial g_z}{\partial z}\right|}{mean\left|\frac{{\partial }^2{g}_z}{\partial z^2}\right|}$E37

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):

$ILP=\arcsin \left[\frac{\sqrt{{\left(\frac{\partial g_z}{\partial x}\right)}^2+{\left(\frac{\partial g_z}{\partial y}\right)}^2}}{\sqrt{{\left(\frac{\partial g_z}{\partial x}\right)}^2+{\left(\frac{\partial g_z}{\partial y}\right)}^2+{\left(\frac{{\partial }^2{g}_z}{\partial x^2}+\frac{{\partial }^2{g}_z}{\partial y^2}\right)}^2}}\right]$E38

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

$TAHG=\arctan \left[\frac{\frac{\partial HG}{\partial z}}{\sqrt{{\left(\frac{\partial HG}{\partial x}\right)}^2+{\left(\frac{\partial HG}{\partial y}\right)}^2}}\right]$E39

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:

$THVH=\arctan \left[\frac{\frac{{\partial }^{2}HG}{\partial z^2}}{\sqrt{{\left(\frac{\partial }{\partial x}\frac{\partial HG}{\partial z}\right)}^2+{\left(\frac{\partial }{\partial y}\frac{\partial HG}{\partial z}\right)}^2}}\right]$E40

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:

$I=\frac{-27I_2^2}{4I_1^3}$E41

Here, I₁ and I₂ are invariants of the tensor. Each invariant is given by three eigenvalues (λ₁, λ₂, λ₃) of the gravity gradient tensors, as follows:

$\begin{array}{l}{I}_1={\lambda }_1{\lambda }_2+{\lambda }_2{\lambda }_3+{\lambda }_1{\lambda }_3\\ I_2={\lambda }_1{\lambda }_2{\lambda }_3\end{array}$E42

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

$Si=\frac{2}{\pi }\arctan \left[\frac{g_{zz}}{\sqrt{{\left(g_{xx}-g_{yy}\right)}^2+4g_{xy}^2}}\right]$E43

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:

$(x-x_0)g_{zx}+(y-y_0)g_{zy}+(z-z_0)g_{zz}=N(B_z-g_z)$E44

In Eq. (44), x₀, y₀, and z₀ 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:

$x_0{g}_{zx}+y_0{g}_{zy}+z_0{g}_{zz}+N{B}_z=x{g}_{zx}+y{g}_{zy}+z{g}_{zz}+N{g}_z$E45

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:

$\begin{array}{c}\left(x-x_0\right)g_{xx}+\left(y-y_0\right)g_{xy}+\left(z-z_0\right)g_{xz}=N\left(B_x-g_x\right)\\ \left(x-x_0\right)g_{yx}+\left(y-y_0\right)g_{yy}+\left(z-z_0\right)g_{yz}=N\left(B_y-g_y\right)\\ \left(x-x_0\right)g_{zx}+\left(y-y_0\right)g_{zy}+\left(z-z_0\right)g_{zz}=N\left(B_z-g_z\right)\end{array}$E46

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₁) of the gradient tensor of the potential field (Figure 3) as follows:

$\alpha =\arctan \left(\frac{v_{1z}}{\sqrt{{\left(v_{1x}\right)}^2+{\left(v_{1y}\right)}^2}}\right)$E47

where v_1x, v_1y, and v_1z are the x, y, and z components of the maximum eigenvector, v₁. 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.

Figure 3.

Schematic illustration of the eigenvectors for a 2D structure such as dike. In this figure, v₁, v₂, and v₃ are the maximum, intermediate, and minimum eigenvectors of the gravity gradient tensor, respectively. The maximum eigenvector points to the causative body and the minimum eigenvector is parallel to the strike direction of the body. The angle α between the surface and maximum eigenvector is the angle of the structural boundary and/or fault.

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.

Figure 4.

Location map of study area. Blue dashed line indicates the outline of the Hohi volcanic zone. Red lines indicate the caldera. Green lines are active faults (e.g., [72]).

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³ of material occurred in the last six million years, and that eruption rate gradually decreased. The volcanic zone erupted 2900 km³/Ma in its initial stage of activity from 5 to 4 Ma. After that, the zone erupted approximately 1300, 900, 400, and 200 km³ 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.

Figure 5.

Residual Bouguer gravity anomaly map of the study area and surroundings. The contour interval is 4 mGal. In the residual Bouguer gravity anomaly map, the regional linear trend of the Bouguer gravity anomalies, such as the effect of the subducting plate, is estimated and removed by the least squares method from the Bouguer gravity anomalies. The Bouguer gravity anomaly used here is based on the gravity anomaly database (1 km × 1 km mesh data) by Komazawa [71], and a Bouguer density of 2670 kg/m³ is employed.

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³ (e.g., [81]) and 90 km³ (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³ 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³. 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³. 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.

Figure 6.

Components of gravitational vector and gravity gradient tensor calculated from g_z data shown in Figure 5, using the method of Mickus and Hinojosa [24]. The gravitational vector unit is given in mGal, and mGal = 10⁻⁵ m/s². The unit of the gravity gradient tensor is given in E (Eötvös), and 1 E = 0.1 mGal/km. (a) g_x component of gravitational vector, (b) g_xx component of gravity gradient tensor, (c) g_xy component, (d) g_xz component, (e) g_y component, (f) g_yx component, which is equal to g_xy, (g) g_yy component, (h) g_yz component, (i) g_z component; gravity anomaly, (j) g_zx component, which is equal to g_xz, (k) g_zy component, which is equal to g_yz, and (l) g_zz component.

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 (∂²g_z/∂z², 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.

Figure 7.

Horizontal gravity gradient map for values over 20 E. Units are given in E (Eötvös). Colored areas indicate locations of structural boundaries such as faults or contacts between different materials.

Figure 8.

Vertical gravity gradient (g_zz) map. Units are given in E (Eötvös). Vertical gravity gradient evaluates vertical gradient of the Bouguer gravity anomaly, and estimates geological or structure boundary by line g_zz = 0.

Figure 9.

Vertical gravity gradient (second differential coefficient: ∂²g_z/∂z², g_zzz) map. Units are given in E/km. E (Eötvös) is unit of the gravity gradient tensor, and 1 E = 0.1 mGal/km. This map indicates geological or structure boundary more sharply than g_zz.

Figure 10.

TDR map. Units are radians. This map indicates geological or structure boundary by line TDR = 0.

Figure 11.

THETA map. Units are radians. This map indicates geological or structure boundary by line THETA = 0.

Figure 12.

TDX map. Units are radians. This map indicates geological or structure boundary by line TDX = max (here, TDX = 1.6).

Figure 13.

CLP map. Units are radians. This map indicates geological or structure boundary by line CLP = max (here, CLP = 1.6).

Figure 14.

ILP map. Units are radians. This map extracts the edge of the gravity anomaly related to multi subsurface structures.

Figure 15.

TAHG map. Units are radians. This map extracts the edge of the gravity anomaly related to multi subsurface structures.

Figure 16.

THVH map. Units are radians. This map extracts the edge of the gravity anomaly related to multi subsurface structures.

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.

Figure 17.

Map of dimensionality index, I. 2D structures are predicted in areas where I < 0.5. Locations of the caldera rim are shown by black lines. The northern small caldera is the Shishimuta caldera and its diameter is 8–10 km. The southern large caldera is the Aso caldera and its diameter is approximately 25 km in the south-north direction and about 18 km in the east-west direction.

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.

Figure 18.

Shape index map. Values of –1, –0.5, 0, +0.5, and +1 indicate bowl, valley, flat, ridge, and dome structures, respectively. Locations of the caldera rim are shown by black lines. The northern small caldera is the Shishimuta caldera and its diameter is 8–10 km. The southern large caldera is the Aso caldera and its diameter is approximately 25 km in the south-north direction and about 18 km in the east-west direction.

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.

Figure 19.

Solution clouds given by the Tensor Euler deconvolution. In the calculation, a 14km width window and a structure index of 0.5 are assumed. (a) Distribution of the solutions in map view. Units are km. Locations of the caldera rim are shown by black lines. The northern small caldera is the Shishimuta caldera and its diameter is 8–10 km. The southern large caldera is the Aso caldera and its diameter is approximately 25 km in the south-north direction and about 18 km in the east-west direction. (b) Distribution of the solutions in bird’s-eye view. Units are km.

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.

Figure 20.

Distribution of estimated dip of faults and/or structural boundaries. Locations of the caldera rim are shown by black lines. The northern small caldera is the Shishimuta caldera and its diameter is 8–10 km. The southern large caldera is the Aso caldera and its diameter is approximately 25 km in the south-north direction and about 18 km in the east-west direction. Dip estimation was conducted for the areas that satisfied the conditions I < 0.5 and HG (horizontal gravity gradient) ≥ 20 E.

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³ 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.

Acknowledgments

This work was supported partially by JSPS (Japan Society for the Promotion of Science) KAKENHI Grant Number 15K14274. The author is grateful to the agency. The author is also most grateful to Agust Gudmundsson for very useful comments that have considerably improved the manuscript. In addition, the author is most grateful to Karoly Nemeth and Dajana Pemac for their editorial advices and cooperation.